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 3432 for branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 – NEMO

Ignore:
Timestamp:
2012-07-11T13:22:58+02:00 (12 years ago)
Author:
trackstand2
Message:

Merge branch 'ksection_partition'

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r3211 r3432  
    3636   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90 
    3737 
     38!FTRANS zdiu zdju zdiv zdjv zdj1u zdj1v zfuw zfvw :I :z 
    3839   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)  
    3940   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      - 
     
    114115      !!---------------------------------------------------------------------- 
    115116      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    116       USE wrk_nemo, ONLY:   ziut  => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt => wrk_2d_3    ! 2D workspace 
    117       USE wrk_nemo, ONLY:   zivf  => wrk_2d_4 , zdku  => wrk_2d_5 , zdkv => wrk_2d_6    ! 2D workspace 
     117#if ! defined key_z_first 
     118      ! Then these workspace arrays can be left as 2D 
     119      USE wrk_nemo, ONLY:   zjvt  => wrk_2d_3    ! 2D workspace 
     120      USE wrk_nemo, ONLY:   zivf  => wrk_2d_4    ! 2D workspace 
     121      USE wrk_nemo, ONLY:   ziut  => wrk_2d_1 , zjuf  => wrk_2d_2 
     122      USE wrk_nemo, ONLY:   zdku  => wrk_2d_5 , zdkv => wrk_2d_6 
    118123      USE wrk_nemo, ONLY:   zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8 
     124#endif 
     125      USE timing,   ONLY:   timing_start, timing_stop 
    119126      ! 
    120127      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     
    125132      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      - 
    126133      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
     134      REAL(wp) ::   zcof, zrecip 
     135#if defined key_z_first 
     136      REAL(wp) ::   zdku, zdk1u, zdki1u, zdk1i1u, zdkim1u, zdk1im1u 
     137      REAL(wp) ::         zdkj1u, zdk1j1u, zdkjm1u, zdk1jm1u 
     138      REAL(wp) ::   zdkv, zdk1v, zdki1v, zdk1i1v, zdkim1v, zdk1im1v 
     139      REAL(wp) ::         zdkj1v, zdk1j1v, zdkjm1v, zdk1jm1v 
     140      REAL(wp) ::   aziut, aziuti1, azjuf, azjufjm1 
     141      REAL(wp) ::   azivf, azivfim1, azjvtj1, azjvt 
     142#endif 
    127143      !!---------------------------------------------------------------------- 
     144 
     145! ziut zjvt zjuf zivf :I :I :z 
     146!!$#if defined key_z_first 
     147!!$      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ziut, zjuf, zivf, zjvt 
     148!!$      ALLOCATE( ziut(jpi,jpj,jpk), zjuf(jpi,jpj,jpk),  & 
     149!!$                zivf(jpi,jpj,jpk), zjvt(jpi,jpj,jpk) ) 
     150!!$#endif 
     151      CALL timing_start('dyn_ldf_iso') 
    128152 
    129153      IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 
     
    170194      ENDIF 
    171195 
     196      CALL timing_start('dyn_ldf_iso_hslab') 
     197 
     198#if defined key_z_first 
     199 
     200      ! Vertical u- and v-shears at level jk and jk+1 
     201      ! --------------------------------------------- 
     202      ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 
     203      !                             zdkv(jk=1)=zdkv(jk=2) 
     204!!$      DO jj = 1, jpj, 1 
     205!!$         DO ji = 1, jpi, 1 
     206!!$ 
     207!!$            ! jk=1 special case 
     208!!$            !zdk1u(ji,jj,1) = ( ub(ji,jj,1) -ub(ji,jj,2) ) * umask(ji,jj,2) 
     209!!$            zdk1v(ji,jj,1) = ( vb(ji,jj,1) -vb(ji,jj,2) ) * vmask(ji,jj,2) 
     210!!$            !zdku(ji,jj,1) = zdk1u(ji,jj,1) 
     211!!$            zdkv(ji,jj,1) = zdk1v(ji,jj,1) 
     212!!$ 
     213!!$            DO jk = 2, jpkm1, 1 
     214!!$               !zdk1u(ji,jj,jk) = ( ub(ji,jj,jk) -ub(ji,jj,jk+1) ) * umask(ji,jj,jk+1) 
     215!!$               zdk1v(ji,jj,jk) = ( vb(ji,jj,jk) -vb(ji,jj,jk+1) ) * vmask(ji,jj,jk+1) 
     216!!$ 
     217!!$               !zdku(ji,jj,jk) = ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) * umask(ji,jj,jk) 
     218!!$               zdkv(ji,jj,jk) = ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) * vmask(ji,jj,jk) 
     219!!$            END DO !jk  
     220!!$         END DO 
     221!!$      END DO 
     222 
     223!!$      !                               -----f----- 
     224!!$      ! Horizontal fluxes on U             |   
     225!!$      ! --------------------===        t   u   t 
     226!!$      !                                    |   
     227!!$      ! i-flux at t-point             -----f----- 
     228!!$ 
     229!!$      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
     230!!$         DO jj = 2, jpjm1 
     231!!$            DO ji = 2, jpi 
     232!!$               DO jk = 1, jpkm1, 1 
     233!!$                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 
     234!!$ 
     235!!$                  zmskt = 1._wp/MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     236!!$                     &              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1._wp ) 
     237!!$ 
     238!!$                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5_wp * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     239!!$    
     240!!$                  ziut(ji,jj,jk) = (  zabe1 * ( ub(ji,jj,jk)    - ub(ji-1,jj,jk) )   & 
     241!!$                     &              + zcof1 * ( zdku  + zdk1im1u     & 
     242!!$                     &                         +zdk1u + zdkim1u )  ) * tmask(ji,jj,jk) 
     243!!$               END DO 
     244!!$            END DO 
     245!!$         END DO 
     246!!$      ELSE                   ! other coordinate system (zco or sco) : e3t 
     247!!$         DO jj = 2, jpjm1 
     248!!$            DO ji = 2, jpi    
     249!!$               DO jk = 1, jpkm1, 1 
     250!!$                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
     251!!$ 
     252!!$                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     253!!$                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
     254!!$ 
     255!!$                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     256!!$ 
     257!!$                  ziut(ji,jj,jk) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
     258!!$                     &              + zcof1 * ( zdku(ji,jj,jk) + zdk1u(ji-1,jj,jk)     & 
     259!!$                     &                         +zdk1u(ji,jj,jk) + zdku (ji-1,jj,jk) )  ) * tmask(ji,jj,jk) 
     260!!$ 
     261!!$               END DO 
     262!!$            END DO 
     263!!$         END DO 
     264!!$      ENDIF 
     265 
     266      ! j-flux at f-point 
     267! BLOCKABLE(ji,jj,jk) 
     268! BLOCKING SIZE (0) 
     269!!$      DO jj = 1, jpjm1, 1 
     270!!$! BLOCKING SIZE (0) 
     271!!$         DO ji = 1, jpim1, 1 
     272!!$! BLOCKING SIZE (4) 
     273!!$            DO jk = 1, jpkm1, 1 
     274!!$               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
     275!!$ 
     276!!$               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
     277!!$                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
     278!!$ 
     279!!$               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     280!!$ 
     281!!$               zjuf(ji,jj,jk) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
     282!!$                  &              + zcof2 * ( zdku (ji,jj+1,jk) + zdk1u(ji,jj,jk)     & 
     283!!$                  &                         +zdk1u(ji,jj+1,jk) + zdku (ji,jj,jk) )  ) * fmask(ji,jj,jk) 
     284!!$            END DO 
     285!!$         END DO 
     286!!$      END DO 
     287 
     288      !                                |   t   | 
     289      ! Horizontal fluxes on V         |       | 
     290      ! --------------------===        f---v---f 
     291      !                                |       | 
     292      !                                |   t   | 
     293       
     294!!$      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
     295!!$         DO jj = 2, jpj 
     296!!$            DO ji = 1, jpim1 
     297!!$               DO jk = 1, jpkm1, 1 
     298!!$ 
     299!!$                  ! i-flux at f-point              |   t   | 
     300!!$                  zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
     301!!$ 
     302!!$                  zmskf = 1._wp/MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
     303!!$                                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) 
     304!!$ 
     305!!$                  zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     306!!$ 
     307!!$                  zivf(ji,jj,jk) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        & 
     308!!$                                    + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk)     & 
     309!!$                                              +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) )  ) * fmask(ji,jj,jk) 
     310!!$ 
     311!!$ 
     312!!$                  ! j-flux at t-point 
     313!!$                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 
     314!!$ 
     315!!$                  zmskt = 1._wp/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     316!!$                                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) 
     317!!$ 
     318!!$                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     319!!$ 
     320!!$                  zjvt(ji,jj,jk) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     321!!$                                    + zcof2 * ( zdkv(ji,jj-1,jk) + zdk1v(ji,jj,jk)     & 
     322!!$                                               +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) )  ) * tmask(ji,jj,jk) 
     323!!$               END DO 
     324!!$            END DO 
     325!!$         END DO 
     326!!$      ELSE                   ! other coordinate system (zco or sco) : e3t 
     327!!$         DO jj = 2, jpj 
     328!!$            DO ji = 1, jpim1 
     329!!$               DO jk = 1, jpkm1 
     330!!$ 
     331!!$                  ! i-flux at f-point  
     332!!$                  zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
     333!!$ 
     334!!$                  zmskf = 1._wp/MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
     335!!$                                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) 
     336!!$ 
     337!!$                  zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     338!!$ 
     339!!$                  zivf(ji,jj,jk) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        & 
     340!!$                                    + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk)     & 
     341!!$                                               +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) )  ) * fmask(ji,jj,jk) 
     342!!$                  ! j-flux at t-point 
     343!!$                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
     344!!$ 
     345!!$                  zmskt = 1._wp/MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     346!!$                                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) 
     347!!$ 
     348!!$                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     349!!$ 
     350!!$                  zjvt(ji,jj,jk) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     351!!$                                    + zcof2 * ( zdkv (ji,jj-1,jk) + zdk1v(ji,jj,jk)     & 
     352!!$                                               +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) )  ) * tmask(ji,jj,jk) 
     353!!$               END DO 
     354!!$            END DO 
     355!!$         END DO 
     356!!$ 
     357!!$      ENDIF 
     358 
     359 
     360 
     361      IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u) 
     362 
     363         DO jj = 2, jpjm1 
     364            DO ji = 2, jpim1 
     365!DIR$ SHORTLOOP 
     366               ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 
     367               ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 
     368               ! when jk is 1 but that doesn't matter. 
     369!DIR$ SAFE_ADDRESS 
     370               DO jk = 1, jpkm1, 1 
     371 
     372                  ! Vertical u- and v-shears at level jk and jk+1 
     373                  ! --------------------------------------------- 
     374                  ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 
     375                  !                             zdkv(jk=1)=zdkv(jk=2) 
     376                  zdk1u    = ( ub(ji  ,jj  ,jk) -ub(ji  ,jj  ,jk+1) ) * umask(ji  ,jj,jk+1) 
     377                  zdk1i1u  = ( ub(ji+1,jj  ,jk) -ub(ji+1,jj  ,jk+1) ) * umask(ji+1,jj,jk+1) 
     378                  zdk1j1u  = ( ub(ji  ,jj+1,jk) -ub(ji  ,jj+1,jk+1) ) * umask(ji  ,jj+1,jk+1) 
     379                  zdk1im1u = ( ub(ji-1,jj  ,jk) -ub(ji-1,jj  ,jk+1) ) * umask(ji-1,jj,jk+1) 
     380                  zdk1jm1u = ( ub(ji  ,jj-1,jk) -ub(ji  ,jj-1,jk+1) ) * umask(ji  ,jj-1,jk+1) 
     381                  IF(jk > 1)THEN 
     382                     zdku   = ( ub(ji  ,jj  ,jk-1) - ub(ji  ,jj  ,jk) ) * umask(ji,jj,jk) 
     383                     zdki1u = ( ub(ji+1,jj  ,jk-1) - ub(ji+1,jj  ,jk) ) * umask(ji+1,jj,jk) 
     384                     zdkim1u= ( ub(ji-1,jj  ,jk-1) - ub(ji-1,jj  ,jk) ) * umask(ji-1,jj,jk) 
     385                     zdkj1u = ( ub(ji  ,jj+1,jk-1) - ub(ji  ,jj+1,jk) ) * umask(ji,jj+1,jk) 
     386                     zdkjm1u= ( ub(ji  ,jj-1,jk-1) - ub(ji  ,jj-1,jk) ) * umask(ji,jj-1,jk) 
     387                  ELSE 
     388                     zdku   = zdk1u 
     389                     zdki1u = zdk1i1u 
     390                     zdkim1u= zdk1im1u 
     391                     zdkj1u = zdk1j1u 
     392                     zdkjm1u= zdk1jm1u 
     393                  END IF 
     394 
     395                  zdku    = zdku + zdk1u 
     396                  zdki1u  = zdki1u + zdk1i1u 
     397                  zdkim1u = zdkim1u + zdk1im1u 
     398                  zdkj1u  = zdkj1u +zdk1j1u 
     399                  zdkjm1u = zdk1jm1u + zdkjm1u 
     400                   
     401                  ! volume elements 
     402                  zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     403 
     404                  ! horizontal component of isopycnal momentum diffusive trends 
     405 
     406                  !                               -----f----- 
     407                  ! Horizontal fluxes on U             |   
     408                  ! --------------------===        t   u   t 
     409                  !                                    |   
     410                  ! i-flux at t-point             -----f----- 
     411 
     412                  ! z-coordinate - partial steps : min(e3u) 
     413                  aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji+1,jj)) * & 
     414                       ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * &  
     415                       (1._wp/MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   & 
     416                               + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1._wp )) * 0.5  * & 
     417                               ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & 
     418                               ( zdku + zdki1u )  ) * tmask(ji+1,jj,jk) 
     419 
     420                  aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & 
     421                    ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * &  
     422                    (1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     423                            + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * & 
     424                    ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & 
     425                    ( zdku + zdkim1u )  ) * tmask(ji,jj,jk) 
     426 
     427                  ! j-flux at f-point - identical to non-ln_zps 
     428 
     429                  azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & 
     430                    ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & 
     431                    (1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
     432                            + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * & 
     433                    ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) *  & 
     434                    ( zdku + zdkj1u )  ) * fmask(ji,jj,jk) 
     435 
     436                  azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & 
     437                    ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & 
     438                    (1./MAX(  umask(ji,jj,jk  )+umask(ji,jj-1,jk+1)   & 
     439                            + umask(ji,jj,jk+1)+umask(ji,jj-1,jk  ), 1. )) * 0.5  * & 
     440                    ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) *  & 
     441                    ( zdku  + zdkjm1u  )  ) * fmask(ji,jj-1,jk) 
     442 
     443                  ! Second derivative (divergence) and add to the general trend 
     444                  ! ----------------------------------------------------------- 
     445 
     446                  zuah =( & 
     447!                       ziut (ji+1,jj,jk) - & 
     448                       aziuti1 - & 
     449!                       ziut (ji  ,jj,jk) + & 
     450                       aziut + & 
     451!                       zjuf (ji  ,jj,jk) - & 
     452                       azjuf - & 
     453!                       zjuf (ji,jj-1,jk)   & 
     454                       azjufjm1 & 
     455                       ) / zbu 
     456 
     457                  ! add the trends to the general trends 
     458                  ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 
     459 
     460               END DO 
     461 
     462!DIR$ SHORTLOOP 
     463! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 
     464! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 
     465! when jk is 1 but that doesn't matter. 
     466!DIR$ SAFE_ADDRESS 
     467               DO jk = 1, jpkm1, 1 
     468 
     469                  zdk1v    = ( vb(ji  ,jj  ,jk) -vb(ji  ,jj  ,jk+1) ) * vmask(ji  ,jj,jk+1) 
     470                  zdk1i1v  = ( vb(ji+1,jj  ,jk) -vb(ji+1,jj  ,jk+1) ) * vmask(ji+1,jj,jk+1) 
     471                  zdk1im1v = ( vb(ji-1,jj  ,jk) -vb(ji-1,jj  ,jk+1) ) * vmask(ji-1,jj,jk+1) 
     472                  zdk1j1v  = ( vb(ji  ,jj+1,jk) -vb(ji  ,jj+1,jk+1) ) * vmask(ji  ,jj+1,jk+1) 
     473                  zdk1jm1v = ( vb(ji  ,jj-1,jk) -vb(ji  ,jj-1,jk+1) ) * vmask(ji  ,jj-1,jk+1) 
     474                  IF(jk > 1)THEN 
     475                     zdkv   = ( vb(ji  ,jj  ,jk-1) - vb(ji  ,jj  ,jk) ) * vmask(ji,jj,jk) 
     476                     zdki1v = ( vb(ji+1,jj  ,jk-1) - vb(ji+1,jj  ,jk) ) * vmask(ji+1,jj,jk) 
     477                     zdkim1v= ( vb(ji-1,jj  ,jk-1) - vb(ji-1,jj  ,jk) ) * vmask(ji-1,jj,jk) 
     478                     zdkj1v = ( vb(ji  ,jj+1,jk-1) - vb(ji  ,jj+1,jk) ) * vmask(ji,jj+1,jk) 
     479                     zdkjm1v= ( vb(ji  ,jj-1,jk-1) - vb(ji  ,jj-1,jk) ) * vmask(ji,jj-1,jk) 
     480                  ELSE 
     481                     zdkv   = zdk1v 
     482                     zdki1v = zdk1i1v 
     483                     zdkim1v= zdk1im1v 
     484                     zdkj1v = zdk1j1v 
     485                     zdkjm1v= zdk1jm1v 
     486                  END IF 
     487                   
     488                  zdkv = zdkv + zdk1v 
     489                  zdki1v = zdk1i1v + zdki1v 
     490                  zdkim1v = zdkim1v +zdk1im1v 
     491                  zdkj1v = zdk1j1v + zdkj1v 
     492                  zdkjm1v = zdkjm1v +zdk1jm1v 
     493 
     494                  !                                |   t   | 
     495                  ! Horizontal fluxes on V         |       | 
     496                  ! --------------------===        f---v---f 
     497                  !                                |       | 
     498                  !                                |   t   | 
     499 
     500                  ! i-flux at f-point - identical to non-ln_zps case 
     501                  azivf = (  (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & 
     502                    e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        & 
     503                    + (- aht0 * e2f(ji,jj) / & 
     504                    MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
     505                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * & 
     506                    0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & 
     507                    ( zdkv + zdki1v )  ) * fmask(ji,jj,jk) 
     508 
     509                  azivfim1 = (  (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & 
     510                    fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & 
     511                    + (- aht0 * e2f(ji-1,jj) / & 
     512                    MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   & 
     513                    + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1._wp ) * 0.5_wp * & 
     514                    ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & 
     515                    ( zdkim1v + zdkv )  ) * fmask(ji-1,jj,jk) 
     516 
     517                  ! j-flux at t-point - min(e3u) instead of e3t 
     518                  azjvtj1 = (  & 
     519                    ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * MIN( fse3v(ji,jj+1,jk), fse3v(ji,jj,jk) ) / & 
     520                    e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) )  & 
     521                    + (- aht0 * e1t(ji,jj+1) / & 
     522                    MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   & 
     523                    + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1._wp ) * & 
     524                    0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & 
     525                    ( zdkv + zdkj1v )  & 
     526                         ) * tmask(ji,jj+1,jk) 
     527 
     528                  azjvt = (  & 
     529                    ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / & 
     530                    e2t(ji,jj)) * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     531                    + (- aht0 * e1t(ji,jj) /MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     532                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * 0.5_wp * & 
     533                    ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 
     534                    ( zdkjm1v + zdkv )  ) * tmask(ji,jj,jk) 
     535 
     536                  ! volume elements 
     537                  zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     538 
     539 
     540                  ! Second derivative (divergence) and add to the general trend 
     541                  ! ----------------------------------------------------------- 
     542                  zvah =( & 
     543!                       zivf (ji  ,jj  ,jk) - & 
     544                       azivf - & 
     545!                       zivf (ji-1,jj  ,jk) + & 
     546                       azivfim1 + & 
     547!                       zjvt (ji  ,jj+1,jk) - & 
     548                       azjvtj1 - & 
     549!                       zjvt (ji  ,jj  ,jk)   & 
     550                       azjvt & 
     551                       ) / zbv 
     552 
     553                  ! add the trend to the general trend 
     554                  va (ji,jj,jk) = va (ji,jj,jk) + zvah 
     555               END DO 
     556            END DO 
     557         END DO 
     558 
     559      ELSE                   ! other coordinate system (zco or sco) : e3t 
     560 
     561         DO jj = 2, jpjm1 
     562            DO ji = 2, jpim1 
     563!DIR$ SHORTLOOP 
     564               ! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 
     565               ! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 
     566               ! when jk is 1 but that doesn't matter. 
     567!DIR$ SAFE_ADDRESS 
     568               DO jk = 1, jpkm1, 1 
     569 
     570                  ! Vertical u- and v-shears at level jk and jk+1 
     571                  ! --------------------------------------------- 
     572                  ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 
     573                  !                             zdkv(jk=1)=zdkv(jk=2) 
     574                  zdk1u    = ( ub(ji  ,jj  ,jk) -ub(ji  ,jj  ,jk+1) ) * umask(ji  ,jj,jk+1) 
     575                  zdk1i1u  = ( ub(ji+1,jj  ,jk) -ub(ji+1,jj  ,jk+1) ) * umask(ji+1,jj,jk+1) 
     576                  zdk1j1u  = ( ub(ji  ,jj+1,jk) -ub(ji  ,jj+1,jk+1) ) * umask(ji  ,jj+1,jk+1) 
     577                  zdk1im1u = ( ub(ji-1,jj  ,jk) -ub(ji-1,jj  ,jk+1) ) * umask(ji-1,jj,jk+1) 
     578                  zdk1jm1u = ( ub(ji  ,jj-1,jk) -ub(ji  ,jj-1,jk+1) ) * umask(ji  ,jj-1,jk+1) 
     579                  IF(jk > 1)THEN 
     580                     zdku   = ( ub(ji  ,jj  ,jk-1) - ub(ji  ,jj  ,jk) ) * umask(ji,jj,jk) 
     581                     zdki1u = ( ub(ji+1,jj  ,jk-1) - ub(ji+1,jj  ,jk) ) * umask(ji+1,jj,jk) 
     582                     zdkim1u= ( ub(ji-1,jj  ,jk-1) - ub(ji-1,jj  ,jk) ) * umask(ji-1,jj,jk) 
     583                     zdkj1u = ( ub(ji  ,jj+1,jk-1) - ub(ji  ,jj+1,jk) ) * umask(ji,jj+1,jk) 
     584                     zdkjm1u= ( ub(ji  ,jj-1,jk-1) - ub(ji  ,jj-1,jk) ) * umask(ji,jj-1,jk) 
     585                  ELSE 
     586                     zdku   = zdk1u 
     587                     zdki1u = zdk1i1u 
     588                     zdkim1u= zdk1im1u 
     589                     zdkj1u = zdk1j1u 
     590                     zdkjm1u= zdk1jm1u 
     591                  END IF 
     592 
     593                  zdku    = zdku + zdk1u 
     594                  zdki1u  = zdki1u + zdk1i1u 
     595                  zdkim1u = zdkim1u + zdk1im1u 
     596                  zdkj1u  = zdkj1u +zdk1j1u 
     597                  zdkjm1u = zdk1jm1u + zdkjm1u 
     598                   
     599                  ! volume element 
     600                  zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     601 
     602                  ! horizontal component of isopycnal momentum diffusive trends 
     603 
     604                  !                               -----f----- 
     605                  ! Horizontal fluxes on U             |   
     606                  ! --------------------===        t   u   t 
     607                  !                                    |   
     608                  ! i-flux at t-point             -----f----- 
     609 
     610                  ! other coordinate system (zco or sco) : e3t 
     611                  aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * & 
     612                       ( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * &  
     613                       (1./MAX(  umask(ji,jj,jk  )+umask(ji+1,jj,jk+1)   & 
     614                               + umask(ji,jj,jk+1)+umask(ji+1,jj,jk  ), 1. )) * 0.5  * & 
     615                               ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * & 
     616                               ( zdku + zdki1u )  ) * tmask(ji+1,jj,jk) 
     617 
     618                  aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * & 
     619                    ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * &  
     620                    (1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
     621                            + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * & 
     622                    ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * & 
     623                    ( zdku + zdkim1u )  ) * tmask(ji,jj,jk) 
     624 
     625                  azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * & 
     626                    ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * & 
     627                    (1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
     628                            + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )) * 0.5  * & 
     629                    ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) *  & 
     630                    ( zdku + zdkj1u )  ) * fmask(ji,jj,jk) 
     631 
     632                  azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * & 
     633                    ( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * & 
     634                    (1./MAX(  umask(ji,jj,jk  )+umask(ji,jj-1,jk+1)   & 
     635                            + umask(ji,jj,jk+1)+umask(ji,jj-1,jk  ), 1. )) * 0.5  * & 
     636                    ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) *  & 
     637                    ( zdku  + zdkjm1u  )  ) * fmask(ji,jj-1,jk) 
     638 
     639 
     640                  ! Second derivative (divergence) and add to the general trend 
     641                  ! ----------------------------------------------------------- 
     642                  zuah =( & 
     643!                       ziut (ji+1,jj,jk) - & 
     644                       aziuti1 - & 
     645!                       ziut (ji  ,jj,jk) + & 
     646                       aziut + & 
     647!                       zjuf (ji  ,jj,jk) - & 
     648                       azjuf - & 
     649!                       zjuf (ji,jj-1,jk)   & 
     650                       azjufjm1 & 
     651                       ) / zbu 
     652 
     653                  ! add the trend to the general trend 
     654                  ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 
     655 
     656               END DO 
     657 
     658!DIR$ SHORTLOOP 
     659! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the 
     660! If condition on jk>1. This does mean that there will be out-of-bounds *reads* 
     661! when jk is 1 but that doesn't matter. 
     662!DIR$ SAFE_ADDRESS 
     663               DO jk = 1, jpkm1, 1 
     664 
     665                  zdk1v    = ( vb(ji  ,jj  ,jk) -vb(ji  ,jj  ,jk+1) ) * vmask(ji  ,jj,jk+1) 
     666                  zdk1i1v  = ( vb(ji+1,jj  ,jk) -vb(ji+1,jj  ,jk+1) ) * vmask(ji+1,jj,jk+1) 
     667                  zdk1im1v = ( vb(ji-1,jj  ,jk) -vb(ji-1,jj  ,jk+1) ) * vmask(ji-1,jj,jk+1) 
     668                  zdk1j1v  = ( vb(ji  ,jj+1,jk) -vb(ji  ,jj+1,jk+1) ) * vmask(ji  ,jj+1,jk+1) 
     669                  zdk1jm1v = ( vb(ji  ,jj-1,jk) -vb(ji  ,jj-1,jk+1) ) * vmask(ji  ,jj-1,jk+1) 
     670                  IF(jk > 1)THEN 
     671                     zdkv   = ( vb(ji  ,jj  ,jk-1) - vb(ji  ,jj  ,jk) ) * vmask(ji,jj,jk) 
     672                     zdki1v = ( vb(ji+1,jj  ,jk-1) - vb(ji+1,jj  ,jk) ) * vmask(ji+1,jj,jk) 
     673                     zdkim1v= ( vb(ji-1,jj  ,jk-1) - vb(ji-1,jj  ,jk) ) * vmask(ji-1,jj,jk) 
     674                     zdkj1v = ( vb(ji  ,jj+1,jk-1) - vb(ji  ,jj+1,jk) ) * vmask(ji,jj+1,jk) 
     675                     zdkjm1v= ( vb(ji  ,jj-1,jk-1) - vb(ji  ,jj-1,jk) ) * vmask(ji,jj-1,jk) 
     676                  ELSE 
     677                     zdkv   = zdk1v 
     678                     zdki1v = zdk1i1v 
     679                     zdkim1v= zdk1im1v 
     680                     zdkj1v = zdk1j1v 
     681                     zdkjm1v= zdk1jm1v 
     682                  END IF 
     683                   
     684                  zdkv = zdkv + zdk1v 
     685                  zdki1v = zdk1i1v + zdki1v 
     686                  zdkim1v = zdkim1v +zdk1im1v 
     687                  zdkj1v = zdk1j1v + zdkj1v 
     688                  zdkjm1v = zdkjm1v +zdk1jm1v 
     689 
     690                  azivf = (  (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / & 
     691                    e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )        & 
     692                    + (- aht0 * e2f(ji,jj) / & 
     693                    MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
     694                    + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * & 
     695                    0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * & 
     696                    ( zdkv + zdki1v )  ) * fmask(ji,jj,jk) 
     697 
     698                  azivfim1 = (  (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * & 
     699                    fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) & 
     700                    + (- aht0 * e2f(ji-1,jj) / & 
     701                    MAX(  vmask(ji,jj,jk  )+vmask(ji-1,jj,jk+1)   & 
     702                    + vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk  ), 1._wp ) * 0.5_wp * & 
     703                    ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * & 
     704                    ( zdkim1v + zdkv )  ) * fmask(ji-1,jj,jk) 
     705 
     706                  azjvtj1 = (  & 
     707                    ((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / & 
     708                    e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) )  & 
     709                    + (- aht0 * e1t(ji,jj+1) / & 
     710                    MAX(  vmask(ji,jj,jk  )+vmask(ji,jj+1,jk+1)   & 
     711                    + vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk  ), 1._wp ) * & 
     712                    0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * & 
     713                    ( zdkv + zdkj1v )  & 
     714                         ) * tmask(ji,jj+1,jk) 
     715 
     716                  azjvt = (  ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * & 
     717                    ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     718                    + (- aht0 * e1t(ji,jj) /MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     719                    + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1._wp ) * 0.5_wp * & 
     720                    ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * & 
     721                    ( zdkjm1v + zdkv )  ) * tmask(ji,jj,jk) 
     722 
     723                  ! volume elements 
     724                  zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     725 
     726                  ! Second derivative (divergence) and add to the general trend 
     727                  ! ----------------------------------------------------------- 
     728                  zvah =( & 
     729!                       zivf (ji  ,jj  ,jk) - & 
     730                       azivf - & 
     731!                       zivf (ji-1,jj  ,jk) + & 
     732                       azivfim1 + & 
     733!                       zjvt (ji  ,jj+1,jk) - & 
     734                       azjvtj1 - & 
     735!                       zjvt (ji  ,jj  ,jk)   & 
     736                       azjvt & 
     737                       ) / zbv 
     738 
     739                  ! add the trend to the general trend 
     740                  va (ji,jj,jk) = va (ji,jj,jk) + zvah 
     741               END DO 
     742            END DO 
     743         END DO 
     744 
     745      END IF ! ln_zps 
     746#else 
    172747      !                                                ! =============== 
    173748      DO jk = 1, jpkm1                                 ! Horizontal slab 
     
    320895      END DO                                           !   End of slab 
    321896      !                                                ! =============== 
     897#endif 
     898      CALL timing_stop('dyn_ldf_iso_hslab','section') 
    322899 
    323900      ! print sum trends (used for debugging) 
     
    325902         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    326903 
    327  
     904      CALL timing_start('dyn_ldf_iso_vslab') 
     905 
     906#if defined key_z_first 
     907      !                                                ! =============== 
     908      DO jj = 2, jpjm1                                 !  Vertical slab 
     909         !                                             ! =============== 
     910 
     911         ! I. vertical trends associated with the lateral mixing 
     912         ! ===================================================== 
     913         !  (excluding the vertical flux proportional to dk[t] 
     914 
     915 
     916         ! I.1 horizontal momentum gradient 
     917         ! -------------------------------- 
     918 
     919         DO ji = 2, jpi 
     920            DO jk = 1, jpk 
     921               ! i-gradient of u at jj 
     922               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) ) 
     923               ! j-gradient of u and v at jj 
     924               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) ) 
     925               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) ) 
     926               ! j-gradient of u and v at jj+1 
     927               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) ) 
     928               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) ) 
     929            END DO 
     930         END DO 
     931 
     932         DO ji = 1, jpim1 
     933            DO jk = 1, jpk 
     934               ! i-gradient of v at jj 
     935               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) ) 
     936            END DO 
     937         END DO 
     938 
     939 
     940         ! I.2 Vertical fluxes 
     941         ! ------------------- 
     942 
     943         ! Surface and bottom vertical fluxes set to zero 
     944!!$         DO ji = 1, jpi 
     945!!$            zfuw(ji, 1 ) = 0.e0 
     946!!$            zfvw(ji, 1 ) = 0.e0 
     947!!$            zfuw(ji,jpk) = 0.e0 
     948!!$            zfvw(ji,jpk) = 0.e0 
     949!!$         END DO 
     950 
     951         ! interior (2=<jk=<jpk-1) on U field 
     952         DO ji = 2, jpim1 
     953            DO jk = 2, jpkm1 
     954               zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 
     955 
     956               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
     957               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 
     958 
     959               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
     960                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
     961               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   & 
     962                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. ) 
     963 
     964               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 
     965               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj 
     966               ! vertical flux on u field 
     967               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     & 
     968                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   & 
     969                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     & 
     970                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) 
     971               ! update avmu (add isopycnal vertical coefficient to avmu) 
     972               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
     973               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 
     974            END DO 
     975         END DO 
     976 
     977         ! interior (2=<jk=<jpk-1) on V field 
     978         DO ji = 2, jpim1 
     979            DO jk = 2, jpkm1 
     980               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 
     981 
     982               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
     983               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 
     984 
     985               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   & 
     986                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. ) 
     987               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   & 
     988                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. ) 
     989 
     990               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi 
     991               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj 
     992               ! vertical flux on v field 
     993               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     & 
     994                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
     995                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
     996                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
     997               ! update avmv (add isopycnal vertical coefficient to avmv) 
     998               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
     999               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 
     1000            END DO 
     1001!         END DO 
     1002 
     1003 
     1004         ! I.3 Divergence of vertical fluxes added to the general tracer trend 
     1005         ! ------------------------------------------------------------------- 
     1006!         DO ji = 2, jpim1 
     1007            zfuw(ji,1)   = 0.0_wp 
     1008            zfuw(ji,jpk) = 0.0_wp 
     1009            zfvw(ji,1)   = 0.0_wp 
     1010            zfvw(ji,jpk) = 0.0_wp 
     1011 
     1012            DO jk = 1, jpkm1 
     1013               ! volume elements 
     1014               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     1015               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     1016               ! part of the k-component of isopycnal momentum diffusive trends 
     1017               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 
     1018               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 
     1019               ! add the trends to the general trends 
     1020               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 
     1021               va(ji,jj,jk) = va(ji,jj,jk) + zvav 
     1022            END DO 
     1023         END DO 
     1024         !                                             ! =============== 
     1025      END DO                                           !   End of vertical slab 
     1026      !                                                ! =============== 
     1027#else 
    3281028      !                                                ! =============== 
    3291029      DO jj = 2, jpjm1                                 !  Vertical slab 
     
    4411141      END DO                                           !   End of slab 
    4421142      !                                                ! =============== 
     1143#endif 
     1144      CALL timing_stop('dyn_ldf_iso_vslab','section') 
    4431145 
    4441146      IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn_ldf_iso: failed to release workspace arrays') 
     1147      ! 
     1148      CALL timing_stop('dyn_ldf_iso','section') 
    4451149      ! 
    4461150   END SUBROUTINE dyn_ldf_iso 
Note: See TracChangeset for help on using the changeset viewer.