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 12377 for NEMO/trunk/src/OCE/LDF/ldftra.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk

    • Property svn:externals
      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
  • NEMO/trunk/src/OCE/LDF/ldftra.F90

    r12296 r12377  
    9494 
    9595   !! * Substitutions 
    96 #  include "vectopt_loop_substitute.h90" 
     96#  include "do_loop_substitute.h90" 
    9797   !!---------------------------------------------------------------------- 
    9898   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    152152      ! ================================= 
    153153      ! 
    154       REWIND( numnam_ref ) 
    155154      READ  ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901) 
    156155901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_ldf in reference namelist' ) 
    157  
    158       REWIND( numnam_cfg ) 
    159156      READ  ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 ) 
    160157902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' ) 
     
    383380 
    384381 
    385    SUBROUTINE ldf_tra( kt ) 
     382   SUBROUTINE ldf_tra( kt, Kbb, Kmm ) 
    386383      !!---------------------------------------------------------------------- 
    387384      !!                  ***  ROUTINE ldf_tra  *** 
     
    404401      !!---------------------------------------------------------------------- 
    405402      INTEGER, INTENT(in) ::   kt   ! time step 
     403      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
    406404      ! 
    407405      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    412410         !                                ! =F(growth rate of baroclinic instability) 
    413411         !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S 
    414          CALL ldf_eiv( kt, aei0, aeiu, aeiv ) 
     412         CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm ) 
    415413      ENDIF 
    416414      ! 
     
    425423            ahtv(:,:,1) = aeiv(:,:,1) 
    426424         ELSE                                            ! compute aht.  
    427             CALL ldf_eiv( kt, aht0, ahtu, ahtv ) 
     425            CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm ) 
    428426         ENDIF 
    429427         ! 
     
    431429         zaht_min = 0.2_wp * aht0                                       ! minimum value for aht 
    432430         zDaht    = aht0 - zaht_min                                       
    433          DO jj = 1, jpj 
    434             DO ji = 1, jpi 
    435                !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 
    436                !!     ==>>>   The Coriolis value is identical for t- & u_points, and for v- and f-points 
    437                zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 
    438                zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 
    439                ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  )     ! min value zaht_min 
    440                ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  )     ! increase within 20S-20N 
    441             END DO 
    442          END DO 
     431         DO_2D_11_11 
     432            !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 
     433            !!     ==>>>   The Coriolis value is identical for t- & u_points, and for v- and f-points 
     434            zaht = ( 1._wp -  MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 
     435            zahf = ( 1._wp -  MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 
     436            ahtu(ji,jj,1) = (  MAX( zaht_min, ahtu(ji,jj,1) ) + zaht  )     ! min value zaht_min 
     437            ahtv(ji,jj,1) = (  MAX( zaht_min, ahtv(ji,jj,1) ) + zahf  )     ! increase within 20S-20N 
     438         END_2D 
    443439         DO jk = 1, jpkm1                             ! deeper value = surface value + mask for all levels 
    444440            ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 
     
    449445         IF( ln_traldf_lap     ) THEN          !   laplacian operator |u| e /12 
    450446            DO jk = 1, jpkm1 
    451                ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12   ! n.b. ub,vb are masked 
    452                ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 
     447               ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12   ! n.b. uu,vv are masked 
     448               ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12 
    453449            END DO 
    454450         ELSEIF( ln_traldf_blp ) THEN      ! bilaplacian operator      sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e 
    455451            DO jk = 1, jpkm1 
    456                ahtu(:,:,jk) = SQRT(  ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12  ) * e1u(:,:) 
    457                ahtv(:,:,jk) = SQRT(  ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12  ) * e2v(:,:) 
     452               ahtu(:,:,jk) = SQRT(  ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12  ) * e1u(:,:) 
     453               ahtv(:,:,jk) = SQRT(  ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12  ) * e2v(:,:) 
    458454            END DO 
    459455         ENDIF 
     
    511507      ENDIF 
    512508      ! 
    513       REWIND( numnam_ref ) 
    514509      READ  ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 
    515510901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' ) 
    516511      ! 
    517       REWIND( numnam_cfg ) 
    518512      READ  ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 
    519513902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' ) 
     
    626620 
    627621 
    628    SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv ) 
     622   SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm ) 
    629623      !!---------------------------------------------------------------------- 
    630624      !!                  ***  ROUTINE ldf_eiv  *** 
     
    638632      !!---------------------------------------------------------------------- 
    639633      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index 
     634      INTEGER                         , INTENT(in   ) ::   Kmm            ! ocean time level indices 
    640635      REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s] 
    641636      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s] 
     
    652647      !                       ! Compute lateral diffusive coefficient at T-point 
    653648      IF( ln_traldf_triad ) THEN 
    654          DO jk = 1, jpk 
    655             DO jj = 2, jpjm1 
    656                DO ji = 2, jpim1 
    657                   ! Take the max of N^2 and zero then take the vertical sum  
    658                   ! of the square root of the resulting N^2 ( required to compute  
    659                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    660                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
    661                   zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 
    662                   ! Compute elements required for the inverse time scale of baroclinic 
    663                   ! eddies using the isopycnal slopes calculated in ldfslp.F :  
    664                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    665                   ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    666                   zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 
    667                   zhw(ji,jj) = zhw(ji,jj) + ze3w 
    668                END DO 
    669             END DO 
    670          END DO 
     649         DO_3D_00_00( 1, jpk ) 
     650            ! Take the max of N^2 and zero then take the vertical sum  
     651            ! of the square root of the resulting N^2 ( required to compute  
     652            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
     653            zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
     654            zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 
     655            ! Compute elements required for the inverse time scale of baroclinic 
     656            ! eddies using the isopycnal slopes calculated in ldfslp.F :  
     657            ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     658            ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     659            zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 
     660            zhw(ji,jj) = zhw(ji,jj) + ze3w 
     661         END_3D 
    671662      ELSE 
    672          DO jk = 1, jpk 
    673             DO jj = 2, jpjm1 
    674                DO ji = 2, jpim1 
    675                   ! Take the max of N^2 and zero then take the vertical sum  
    676                   ! of the square root of the resulting N^2 ( required to compute  
    677                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    678                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
    679                   zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 
    680                   ! Compute elements required for the inverse time scale of baroclinic 
    681                   ! eddies using the isopycnal slopes calculated in ldfslp.F :  
    682                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    683                   ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    684                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    685                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 
    686                   zhw(ji,jj) = zhw(ji,jj) + ze3w 
    687                END DO 
    688             END DO 
    689          END DO 
    690       ENDIF 
    691  
    692       DO jj = 2, jpjm1 
    693          DO ji = fs_2, fs_jpim1   ! vector opt. 
    694             zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 
    695             ! Rossby radius at w-point taken betwenn 2 km and  40km 
    696             zRo(ji,jj) = MAX(  2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 )  ) 
    697             ! Compute aeiw by multiplying Ro^2 and T^-1 
    698             zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
    699          END DO 
    700       END DO 
     663         DO_3D_00_00( 1, jpk ) 
     664            ! Take the max of N^2 and zero then take the vertical sum  
     665            ! of the square root of the resulting N^2 ( required to compute  
     666            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
     667            zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
     668            zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 
     669            ! Compute elements required for the inverse time scale of baroclinic 
     670            ! eddies using the isopycnal slopes calculated in ldfslp.F :  
     671            ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
     672            ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     673            zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     674               &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 
     675            zhw(ji,jj) = zhw(ji,jj) + ze3w 
     676         END_3D 
     677      ENDIF 
     678 
     679      DO_2D_00_00 
     680         zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 
     681         ! Rossby radius at w-point taken betwenn 2 km and  40km 
     682         zRo(ji,jj) = MAX(  2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 )  ) 
     683         ! Compute aeiw by multiplying Ro^2 and T^-1 
     684         zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
     685      END_2D 
    701686 
    702687      !                                         !==  Bound on eiv coeff.  ==! 
    703688      z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  ) 
    704       DO jj = 2, jpjm1 
    705          DO ji = fs_2, fs_jpim1   ! vector opt. 
    706             zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease 
    707             zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0 
    708          END DO 
    709       END DO 
     689      DO_2D_00_00 
     690         zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease 
     691         zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0 
     692      END_2D 
    710693      CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. )       ! lateral boundary condition 
    711694      !                
    712       DO jj = 2, jpjm1                          !== aei at u- and v-points  ==! 
    713          DO ji = fs_2, fs_jpim1   ! vector opt. 
    714             paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1) 
    715             paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1) 
    716          END DO  
    717       END DO  
     695      DO_2D_00_00 
     696         paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj  ) ) * umask(ji,jj,1) 
     697         paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji  ,jj+1) ) * vmask(ji,jj,1) 
     698      END_2D 
    718699      CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. )      ! lateral boundary condition 
    719700 
     
    726707 
    727708 
    728    SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype ) 
     709   SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm, Krhs ) 
    729710      !!---------------------------------------------------------------------- 
    730711      !!                  ***  ROUTINE ldf_eiv_trp  *** 
     
    742723      !!                                    velocity and heat transport (call ldf_eiv_dia) 
    743724      !! 
    744       !! ** Action  : pun, pvn increased by the eiv transport 
    745       !!---------------------------------------------------------------------- 
    746       INTEGER                         , INTENT(in   ) ::   kt       ! ocean time-step index 
    747       INTEGER                         , INTENT(in   ) ::   kit000   ! first time step index 
    748       CHARACTER(len=3)                , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    749       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean transport components   [m3/s] 
    750       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean transport components   [m3/s] 
    751       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv                [m3/s] 
     725      !! ** Action  : pu, pv increased by the eiv transport 
     726      !!---------------------------------------------------------------------- 
     727      INTEGER                         , INTENT(in   ) ::   kt        ! ocean time-step index 
     728      INTEGER                         , INTENT(in   ) ::   kit000    ! first time step index 
     729      INTEGER                         , INTENT(in   ) ::   Kmm, Krhs ! ocean time level indices 
     730      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype    ! =TRA or TRC (tracer indicator) 
     731      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu      ! in : 3 ocean transport components   [m3/s] 
     732      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv      ! out: 3 ocean transport components   [m3/s] 
     733      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw      ! increased by the eiv                [m3/s] 
    752734      !! 
    753735      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     
    767749      zpsi_uw(:,:,jpk) = 0._wp   ;   zpsi_vw(:,:,jpk) = 0._wp 
    768750      ! 
    769       DO jk = 2, jpkm1 
    770          DO jj = 1, jpjm1 
    771             DO ji = 1, fs_jpim1   ! vector opt. 
    772                zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   & 
    773                   &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * wumask(ji,jj,jk) 
    774                zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   & 
    775                   &                                    * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * wvmask(ji,jj,jk) 
    776             END DO 
    777          END DO 
    778       END DO 
    779       ! 
    780       DO jk = 1, jpkm1 
    781          DO jj = 1, jpjm1 
    782             DO ji = 1, fs_jpim1   ! vector opt.                
    783                pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
    784                pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
    785             END DO 
    786          END DO 
    787       END DO 
    788       DO jk = 1, jpkm1 
    789          DO jj = 2, jpjm1 
    790             DO ji = fs_2, fs_jpim1   ! vector opt. 
    791                pwn(ji,jj,jk) = pwn(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   & 
    792                   &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) ) 
    793             END DO 
    794          END DO 
    795       END DO 
     751      DO_3D_10_10( 2, jpkm1 ) 
     752         zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk) )   & 
     753            &                                    * ( aeiu (ji,jj,jk-1) + aeiu (ji  ,jj,jk) ) * wumask(ji,jj,jk) 
     754         zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk) )   & 
     755            &                                    * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj  ,jk) ) * wvmask(ji,jj,jk) 
     756      END_3D 
     757      ! 
     758      DO_3D_10_10( 1, jpkm1 ) 
     759         pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
     760         pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
     761      END_3D 
     762      DO_3D_00_00( 1, jpkm1 ) 
     763         pw(ji,jj,jk) = pw(ji,jj,jk) + (  zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj  ,jk)   & 
     764            &                             + zpsi_vw(ji,jj,jk) - zpsi_vw(ji  ,jj-1,jk) ) 
     765      END_3D 
    796766      ! 
    797767      !                              ! diagnose the eddy induced velocity and associated heat transport 
    798       IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw ) 
     768      IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 
    799769      ! 
    800770    END SUBROUTINE ldf_eiv_trp 
    801771 
    802772 
    803    SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw ) 
     773   SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm ) 
    804774      !!---------------------------------------------------------------------- 
    805775      !!                  ***  ROUTINE ldf_eiv_dia  *** 
     
    812782      !!---------------------------------------------------------------------- 
    813783      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   psi_uw, psi_vw   ! streamfunction   [m3/s] 
     784      INTEGER                         , INTENT(in   ) ::   Kmm   ! ocean time level indices 
    814785      ! 
    815786      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     
    832803      ! 
    833804      DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw] 
    834          zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) ) 
     805         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) ) 
    835806      END DO 
    836807      CALL iom_put( "uoce_eiv", zw3d ) 
    837808      ! 
    838809      DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw] 
    839          zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) ) 
     810         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) ) 
    840811      END DO 
    841812      CALL iom_put( "voce_eiv", zw3d ) 
    842813      ! 
    843       DO jk = 1, jpkm1                                         ! e1 e2 w_eiv = dk[psix] + dk[psix] 
    844          DO jj = 2, jpjm1 
    845             DO ji = fs_2, fs_jpim1  ! vector opt. 
    846                zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    & 
    847                   &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj) 
    848             END DO 
    849          END DO 
    850       END DO 
     814      DO_3D_00_00( 1, jpkm1 ) 
     815         zw3d(ji,jj,jk) = (  psi_vw(ji,jj,jk) - psi_vw(ji  ,jj-1,jk)    & 
     816            &              + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj  ,jk)  ) / e1e2t(ji,jj) 
     817      END_3D 
    851818      CALL lbc_lnk( 'ldftra', zw3d, 'T', 1. )      ! lateral boundary condition 
    852819      CALL iom_put( "woce_eiv", zw3d ) 
     
    872839        zw2d(:,:)   = 0._wp  
    873840        zw3d(:,:,:) = 0._wp  
    874         DO jk = 1, jpkm1 
    875            DO jj = 2, jpjm1 
    876               DO ji = fs_2, fs_jpim1   ! vector opt. 
    877                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   & 
    878                     &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji+1,jj,jk,jp_tem) )  
    879                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
    880               END DO 
    881            END DO 
    882         END DO 
     841        DO_3D_00_00( 1, jpkm1 ) 
     842           zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   & 
     843              &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji+1,jj,jk,jp_tem,Kmm) )  
     844           zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     845        END_3D 
    883846        CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 
    884847        CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 
     
    897860      zw2d(:,:)   = 0._wp  
    898861      zw3d(:,:,:) = 0._wp  
    899       DO jk = 1, jpkm1 
    900          DO jj = 2, jpjm1 
    901             DO ji = fs_2, fs_jpim1   ! vector opt. 
    902                zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   & 
    903                   &                            * ( tsn   (ji,jj,jk,jp_tem) + tsn   (ji,jj+1,jk,jp_tem) )  
    904                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
    905             END DO 
    906          END DO 
    907       END DO 
     862      DO_3D_00_00( 1, jpkm1 ) 
     863         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   & 
     864            &                            * ( ts    (ji,jj,jk,jp_tem,Kmm) + ts    (ji,jj+1,jk,jp_tem,Kmm) )  
     865         zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     866      END_3D 
    908867      CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 
    909868      CALL iom_put( "veiv_heattr", zztmp * zw2d )                  !  heat transport in j-direction 
    910869      CALL iom_put( "veiv_heattr", zztmp * zw3d )                  !  heat transport in j-direction 
    911870      ! 
    912       IF( ln_diaptr )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 
     871      IF( iom_use( 'sophteiv' ) )   CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 
    913872      ! 
    914873      zztmp = 0.5_wp * 0.5 
     
    916875        zw2d(:,:) = 0._wp  
    917876        zw3d(:,:,:) = 0._wp  
    918         DO jk = 1, jpkm1 
    919            DO jj = 2, jpjm1 
    920               DO ji = fs_2, fs_jpim1   ! vector opt. 
    921                  zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)      - psi_uw(ji,jj,jk)          )   & 
    922                     &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji+1,jj,jk,jp_sal) )  
    923                  zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
    924               END DO 
    925            END DO 
    926         END DO 
     877        DO_3D_00_00( 1, jpkm1 ) 
     878           zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1)          - psi_uw(ji  ,jj,jk)            )   & 
     879              &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji+1,jj,jk,jp_sal,Kmm) )  
     880           zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     881        END_3D 
    927882        CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 
    928883        CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 
    929884        CALL iom_put( "ueiv_salttr", zztmp * zw2d )                  ! salt transport in i-direction 
    930         CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                  ! salt transport in i-direction 
     885        CALL iom_put( "ueiv_salttr3d", zztmp * zw3d )                ! salt transport in i-direction 
    931886      ENDIF 
    932887      zw2d(:,:) = 0._wp  
    933888      zw3d(:,:,:) = 0._wp  
    934       DO jk = 1, jpkm1 
    935          DO jj = 2, jpjm1 
    936             DO ji = fs_2, fs_jpim1   ! vector opt. 
    937                zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)      - psi_vw(ji,jj,jk)          )   & 
    938                   &                            * ( tsn   (ji,jj,jk,jp_sal) + tsn   (ji,jj+1,jk,jp_sal) )  
    939                zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
    940             END DO 
    941          END DO 
    942       END DO 
     889      DO_3D_00_00( 1, jpkm1 ) 
     890         zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1)          - psi_vw(ji,jj  ,jk)            )   & 
     891            &                            * ( ts    (ji,jj,jk,jp_sal,Kmm) + ts    (ji,jj+1,jk,jp_sal,Kmm) )  
     892         zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 
     893      END_3D 
    943894      CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 
    944895      CALL iom_put( "veiv_salttr", zztmp * zw2d )                  !  salt transport in j-direction 
    945896      CALL iom_put( "veiv_salttr", zztmp * zw3d )                  !  salt transport in j-direction 
    946897      ! 
    947       IF( ln_diaptr ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 
     898      IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 
    948899      ! 
    949900      ! 
Note: See TracChangeset for help on using the changeset viewer.