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/TRA/eosbn2.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/TRA/eosbn2.F90

    r11993 r12377  
    2929   !!   eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass 
    3030   !!   eos_insitu_2d : Compute the in situ density for 2d fields 
    31    !!   bn2           : Compute the Brunt-Vaisala frequency 
    3231   !!   bn2           : compute the Brunt-Vaisala frequency 
    3332   !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature 
    … …  
    180179 
    181180   !! * Substitutions 
    182 #  include "vectopt_loop_substitute.h90" 
     181#  include "do_loop_substitute.h90" 
    183182   !!---------------------------------------------------------------------- 
    184183   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    238237      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    239238         ! 
    240          DO jk = 1, jpkm1 
    241             DO jj = 1, jpj 
    242                DO ji = 1, jpi 
    243                   ! 
    244                   zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    245                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    246                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    247                   ztm = tmask(ji,jj,jk)                                         ! tmask 
     239         DO_3D_11_11( 1, jpkm1 ) 
     240            ! 
     241            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     242            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     243            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     244            ztm = tmask(ji,jj,jk)                                         ! tmask 
     245            ! 
     246            zn3 = EOS013*zt   & 
     247               &   + EOS103*zs+EOS003 
     248               ! 
     249            zn2 = (EOS022*zt   & 
     250               &   + EOS112*zs+EOS012)*zt   & 
     251               &   + (EOS202*zs+EOS102)*zs+EOS002 
     252               ! 
     253            zn1 = (((EOS041*zt   & 
     254               &   + EOS131*zs+EOS031)*zt   & 
     255               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     256               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     257               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     258               ! 
     259            zn0 = (((((EOS060*zt   & 
     260               &   + EOS150*zs+EOS050)*zt   & 
     261               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     262               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     263               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     264               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     265               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     266               ! 
     267            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     268            ! 
     269            prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     270            ! 
     271         END_3D 
     272         ! 
     273      CASE( np_seos )                !==  simplified EOS  ==! 
     274         ! 
     275         DO_3D_11_11( 1, jpkm1 ) 
     276            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     277            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     278            zh  = pdep (ji,jj,jk) 
     279            ztm = tmask(ji,jj,jk) 
     280            ! 
     281            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     282               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     283               &  - rn_nu * zt * zs 
     284               !                                  
     285            prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked) 
     286         END_3D 
     287         ! 
     288      END SELECT 
     289      ! 
     290      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
     291      ! 
     292      IF( ln_timing )   CALL timing_stop('eos-insitu') 
     293      ! 
     294   END SUBROUTINE eos_insitu 
     295 
     296 
     297   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
     298      !!---------------------------------------------------------------------- 
     299      !!                  ***  ROUTINE eos_insitu_pot  *** 
     300      !! 
     301      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
     302      !!      potential volumic mass (Kg/m3) from potential temperature and 
     303      !!      salinity fields using an equation of state selected in the 
     304      !!     namelist. 
     305      !! 
     306      !! ** Action  : - prd  , the in situ density (no units) 
     307      !!              - prhop, the potential volumic mass (Kg/m3) 
     308      !! 
     309      !!---------------------------------------------------------------------- 
     310      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     311      !                                                                ! 2 : salinity               [psu] 
     312      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
     313      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     314      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
     315      ! 
     316      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     317      INTEGER  ::   jdof 
     318      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     319      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     320      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
     321      !!---------------------------------------------------------------------- 
     322      ! 
     323      IF( ln_timing )   CALL timing_start('eos-pot') 
     324      ! 
     325      SELECT CASE ( neos ) 
     326      ! 
     327      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     328         ! 
     329         ! Stochastic equation of state 
     330         IF ( ln_sto_eos ) THEN 
     331            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     332            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     333            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     334            DO jsmp = 1, 2*nn_sto_eos, 2 
     335              zsign(jsmp)   = 1._wp 
     336              zsign(jsmp+1) = -1._wp 
     337            END DO 
     338            ! 
     339            DO_3D_11_11( 1, jpkm1 ) 
     340               ! 
     341               ! compute density (2*nn_sto_eos) times: 
     342               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     343               ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     344               DO jsmp = 1, nn_sto_eos*2 
     345                  jdof   = (jsmp + 1) / 2 
     346                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     347                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     348                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     349                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     350                  ztm    = tmask(ji,jj,jk)                                         ! tmask 
    248351                  ! 
    249352                  zn3 = EOS013*zt   & 
    … …  
    260363                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    261364                     ! 
    262                   zn0 = (((((EOS060*zt   & 
     365                  zn0_sto(jsmp) = (((((EOS060*zt   & 
    263366                     &   + EOS150*zs+EOS050)*zt   & 
    264367                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    … …  
    268371                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    269372                     ! 
    270                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     373                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     374               END DO 
     375               ! 
     376               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     377               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     378               DO jsmp = 1, nn_sto_eos*2 
     379                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    271380                  ! 
    272                   prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked) 
    273                   ! 
     381                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
    274382               END DO 
    275             END DO 
    276          END DO 
    277          ! 
    278       CASE( np_seos )                !==  simplified EOS  ==! 
    279          ! 
    280          DO jk = 1, jpkm1 
    281             DO jj = 1, jpj 
    282                DO ji = 1, jpi 
    283                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    284                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    285                   zh  = pdep (ji,jj,jk) 
    286                   ztm = tmask(ji,jj,jk) 
    287                   ! 
    288                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    289                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    290                      &  - rn_nu * zt * zs 
    291                      !                                  
    292                   prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked) 
    293                END DO 
    294             END DO 
    295          END DO 
    296          ! 
    297       END SELECT 
    298       ! 
    299       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk ) 
    300       ! 
    301       IF( ln_timing )   CALL timing_stop('eos-insitu') 
    302       ! 
    303    END SUBROUTINE eos_insitu 
    304  
    305  
    306    SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
    307       !!---------------------------------------------------------------------- 
    308       !!                  ***  ROUTINE eos_insitu_pot  *** 
    309       !! 
    310       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
    311       !!      potential volumic mass (Kg/m3) from potential temperature and 
    312       !!      salinity fields using an equation of state selected in the 
    313       !!     namelist. 
    314       !! 
    315       !! ** Action  : - prd  , the in situ density (no units) 
    316       !!              - prhop, the potential volumic mass (Kg/m3) 
    317       !! 
    318       !!---------------------------------------------------------------------- 
    319       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
    320       !                                                                ! 2 : salinity               [psu] 
    321       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-] 
    322       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
    323       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    324       ! 
    325       INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
    326       INTEGER  ::   jdof 
    327       REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
    328       REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
    329       REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    330       !!---------------------------------------------------------------------- 
    331       ! 
    332       IF( ln_timing )   CALL timing_start('eos-pot') 
    333       ! 
    334       SELECT CASE ( neos ) 
    335       ! 
    336       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    337          ! 
    338          ! Stochastic equation of state 
    339          IF ( ln_sto_eos ) THEN 
    340             ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
    341             ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
    342             ALLOCATE(zsign(1:2*nn_sto_eos)) 
    343             DO jsmp = 1, 2*nn_sto_eos, 2 
    344               zsign(jsmp)   = 1._wp 
    345               zsign(jsmp+1) = -1._wp 
    346             END DO 
    347             ! 
    348             DO jk = 1, jpkm1 
    349                DO jj = 1, jpj 
    350                   DO ji = 1, jpi 
    351                      ! 
    352                      ! compute density (2*nn_sto_eos) times: 
    353                      ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
    354                      ! (2) for t-dt, s-ds (with the opposite fluctuation) 
    355                      DO jsmp = 1, nn_sto_eos*2 
    356                         jdof   = (jsmp + 1) / 2 
    357                         zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    358                         zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
    359                         zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
    360                         zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
    361                         ztm    = tmask(ji,jj,jk)                                         ! tmask 
    362                         ! 
    363                         zn3 = EOS013*zt   & 
    364                            &   + EOS103*zs+EOS003 
    365                            ! 
    366                         zn2 = (EOS022*zt   & 
    367                            &   + EOS112*zs+EOS012)*zt   & 
    368                            &   + (EOS202*zs+EOS102)*zs+EOS002 
    369                            ! 
    370                         zn1 = (((EOS041*zt   & 
    371                            &   + EOS131*zs+EOS031)*zt   & 
    372                            &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    373                            &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    374                            &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    375                            ! 
    376                         zn0_sto(jsmp) = (((((EOS060*zt   & 
    377                            &   + EOS150*zs+EOS050)*zt   & 
    378                            &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    379                            &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    380                            &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    381                            &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    382                            &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    383                            ! 
    384                         zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
    385                      END DO 
    386                      ! 
    387                      ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
    388                      prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
    389                      DO jsmp = 1, nn_sto_eos*2 
    390                         prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    391                         ! 
    392                         prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
    393                      END DO 
    394                      prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
    395                      prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
    396                   END DO 
    397                END DO 
    398             END DO 
     383               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     384               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     385            END_3D 
    399386            DEALLOCATE(zn0_sto,zn_sto,zsign) 
    400387         ! Non-stochastic equation of state 
    401388         ELSE 
    402             DO jk = 1, jpkm1 
    403                DO jj = 1, jpj 
    404                   DO ji = 1, jpi 
    405                      ! 
    406                      zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    407                      zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    408                      zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    409                      ztm = tmask(ji,jj,jk)                                         ! tmask 
    410                      ! 
    411                      zn3 = EOS013*zt   & 
    412                         &   + EOS103*zs+EOS003 
    413                         ! 
    414                      zn2 = (EOS022*zt   & 
    415                         &   + EOS112*zs+EOS012)*zt   & 
    416                         &   + (EOS202*zs+EOS102)*zs+EOS002 
    417                         ! 
    418                      zn1 = (((EOS041*zt   & 
    419                         &   + EOS131*zs+EOS031)*zt   & 
    420                         &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    421                         &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    422                         &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    423                         ! 
    424                      zn0 = (((((EOS060*zt   & 
    425                         &   + EOS150*zs+EOS050)*zt   & 
    426                         &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    427                         &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    428                         &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    429                         &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    430                         &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    431                         ! 
    432                      zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    433                      ! 
    434                      prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    435                      ! 
    436                      prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
    437                   END DO 
    438                END DO 
    439             END DO 
    440          ENDIF 
    441           
    442       CASE( np_seos )                !==  simplified EOS  ==! 
    443          ! 
    444          DO jk = 1, jpkm1 
    445             DO jj = 1, jpj 
    446                DO ji = 1, jpi 
    447                   zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
    448                   zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
    449                   zh  = pdep (ji,jj,jk) 
    450                   ztm = tmask(ji,jj,jk) 
    451                   !                                                     ! potential density referenced at the surface 
    452                   zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
    453                      &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
    454                      &  - rn_nu * zt * zs 
    455                   prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 
    456                   !                                                     ! density anomaly (masked) 
    457                   zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
    458                   prd(ji,jj,jk) = zn * r1_rau0 * ztm 
    459                   ! 
    460                END DO 
    461             END DO 
    462          END DO 
    463          ! 
    464       END SELECT 
    465       ! 
    466       IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
    467       ! 
    468       IF( ln_timing )   CALL timing_stop('eos-pot') 
    469       ! 
    470    END SUBROUTINE eos_insitu_pot 
    471  
    472  
    473    SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
    474       !!---------------------------------------------------------------------- 
    475       !!                  ***  ROUTINE eos_insitu_2d  *** 
    476       !! 
    477       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    478       !!      potential temperature and salinity using an equation of state 
    479       !!      selected in the nameos namelist. * 2D field case 
    480       !! 
    481       !! ** Action  : - prd , the in situ density (no units) (unmasked) 
    482       !! 
    483       !!---------------------------------------------------------------------- 
    484       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
    485       !                                                           ! 2 : salinity               [psu] 
    486       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
    487       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    488       ! 
    489       INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    490       REAL(wp) ::   zt , zh , zs              ! local scalars 
    491       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
    492       !!---------------------------------------------------------------------- 
    493       ! 
    494       IF( ln_timing )   CALL timing_start('eos2d') 
    495       ! 
    496       prd(:,:) = 0._wp 
    497       ! 
    498       SELECT CASE( neos ) 
    499       ! 
    500       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    501          ! 
    502          DO jj = 1, jpjm1 
    503             DO ji = 1, fs_jpim1   ! vector opt. 
    504                ! 
    505                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    506                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    507                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     389            DO_3D_11_11( 1, jpkm1 ) 
     390               ! 
     391               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     392               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     393               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     394               ztm = tmask(ji,jj,jk)                                         ! tmask 
    508395               ! 
    509396               zn3 = EOS013*zt   & 
    … …  
    530417               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    531418               ! 
    532                prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly 
    533                ! 
    534             END DO 
    535          END DO 
    536          ! 
    537          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
    538          ! 
     419               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     420               ! 
     421               prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     422            END_3D 
     423         ENDIF 
     424          
    539425      CASE( np_seos )                !==  simplified EOS  ==! 
    540426         ! 
    541          DO jj = 1, jpjm1 
    542             DO ji = 1, fs_jpim1   ! vector opt. 
    543                ! 
    544                zt    = pts  (ji,jj,jp_tem)  - 10._wp 
    545                zs    = pts  (ji,jj,jp_sal)  - 35._wp 
    546                zh    = pdep (ji,jj)                         ! depth at the partial step level 
    547                ! 
    548                zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
    549                   &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
    550                   &  - rn_nu * zt * zs 
    551                   ! 
    552                prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly 
    553                ! 
    554             END DO 
    555          END DO 
    556          ! 
    557          CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions 
     427         DO_3D_11_11( 1, jpkm1 ) 
     428            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     429            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     430            zh  = pdep (ji,jj,jk) 
     431            ztm = tmask(ji,jj,jk) 
     432            !                                                     ! potential density referenced at the surface 
     433            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
     434               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
     435               &  - rn_nu * zt * zs 
     436            prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 
     437            !                                                     ! density anomaly (masked) 
     438            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
     439            prd(ji,jj,jk) = zn * r1_rau0 * ztm 
     440            ! 
     441         END_3D 
    558442         ! 
    559443      END SELECT 
    560444      ! 
    561       IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
     445      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk ) 
     446      ! 
     447      IF( ln_timing )   CALL timing_stop('eos-pot') 
     448      ! 
     449   END SUBROUTINE eos_insitu_pot 
     450 
     451 
     452   SUBROUTINE eos_insitu_2d( pts, pdep, prd ) 
     453      !!---------------------------------------------------------------------- 
     454      !!                  ***  ROUTINE eos_insitu_2d  *** 
     455      !! 
     456      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
     457      !!      potential temperature and salinity using an equation of state 
     458      !!      selected in the nameos namelist. * 2D field case 
     459      !! 
     460      !! ** Action  : - prd , the in situ density (no units) (unmasked) 
     461      !! 
     462      !!---------------------------------------------------------------------- 
     463      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     464      !                                                           ! 2 : salinity               [psu] 
     465      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m] 
     466      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
     467      ! 
     468      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     469      REAL(wp) ::   zt , zh , zs              ! local scalars 
     470      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     471      !!---------------------------------------------------------------------- 
     472      ! 
     473      IF( ln_timing )   CALL timing_start('eos2d') 
     474      ! 
     475      prd(:,:) = 0._wp 
     476      ! 
     477      SELECT CASE( neos ) 
     478      ! 
     479      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     480         ! 
     481         DO_2D_11_11 
     482            ! 
     483            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     484            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     485            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     486            ! 
     487            zn3 = EOS013*zt   & 
     488               &   + EOS103*zs+EOS003 
     489               ! 
     490            zn2 = (EOS022*zt   & 
     491               &   + EOS112*zs+EOS012)*zt   & 
     492               &   + (EOS202*zs+EOS102)*zs+EOS002 
     493               ! 
     494            zn1 = (((EOS041*zt   & 
     495               &   + EOS131*zs+EOS031)*zt   & 
     496               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     497               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     498               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     499               ! 
     500            zn0 = (((((EOS060*zt   & 
     501               &   + EOS150*zs+EOS050)*zt   & 
     502               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     503               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     504               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     505               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     506               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     507               ! 
     508            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     509            ! 
     510            prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly 
     511            ! 
     512         END_2D 
     513         ! 
     514      CASE( np_seos )                !==  simplified EOS  ==! 
     515         ! 
     516         DO_2D_11_11 
     517            ! 
     518            zt    = pts  (ji,jj,jp_tem)  - 10._wp 
     519            zs    = pts  (ji,jj,jp_sal)  - 35._wp 
     520            zh    = pdep (ji,jj)                         ! depth at the partial step level 
     521            ! 
     522            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     523               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     524               &  - rn_nu * zt * zs 
     525               ! 
     526            prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly 
     527            ! 
     528         END_2D 
     529         ! 
     530      END SELECT 
     531      ! 
     532      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' ) 
    562533      ! 
    563534      IF( ln_timing )   CALL timing_stop('eos2d') 
    … …  
    566537 
    567538 
    568    SUBROUTINE rab_3d( pts, pab ) 
     539   SUBROUTINE rab_3d( pts, pab, Kmm ) 
    569540      !!---------------------------------------------------------------------- 
    570541      !!                 ***  ROUTINE rab_3d  *** 
    … …  
    576547      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    577548      !!---------------------------------------------------------------------- 
     549      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    578550      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
    579551      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
    … …  
    590562      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    591563         ! 
    592          DO jk = 1, jpkm1 
    593             DO jj = 1, jpj 
    594                DO ji = 1, jpi 
    595                   ! 
    596                   zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth 
    597                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    598                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    599                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    600                   ! 
    601                   ! alpha 
    602                   zn3 = ALP003 
    603                   ! 
    604                   zn2 = ALP012*zt + ALP102*zs+ALP002 
    605                   ! 
    606                   zn1 = ((ALP031*zt   & 
    607                      &   + ALP121*zs+ALP021)*zt   & 
    608                      &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    609                      &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    610                      ! 
    611                   zn0 = ((((ALP050*zt   & 
    612                      &   + ALP140*zs+ALP040)*zt   & 
    613                      &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    614                      &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    615                      &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    616                      &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    617                      ! 
    618                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    619                   ! 
    620                   pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 
    621                   ! 
    622                   ! beta 
    623                   zn3 = BET003 
    624                   ! 
    625                   zn2 = BET012*zt + BET102*zs+BET002 
    626                   ! 
    627                   zn1 = ((BET031*zt   & 
    628                      &   + BET121*zs+BET021)*zt   & 
    629                      &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    630                      &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    631                      ! 
    632                   zn0 = ((((BET050*zt   & 
    633                      &   + BET140*zs+BET040)*zt   & 
    634                      &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    635                      &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    636                      &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    637                      &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    638                      ! 
    639                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    640                   ! 
    641                   pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 
    642                   ! 
    643                END DO 
    644             END DO 
    645          END DO 
     564         DO_3D_11_11( 1, jpkm1 ) 
     565            ! 
     566            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     567            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     568            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     569            ztm = tmask(ji,jj,jk)                                         ! tmask 
     570            ! 
     571            ! alpha 
     572            zn3 = ALP003 
     573            ! 
     574            zn2 = ALP012*zt + ALP102*zs+ALP002 
     575            ! 
     576            zn1 = ((ALP031*zt   & 
     577               &   + ALP121*zs+ALP021)*zt   & 
     578               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     579               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     580               ! 
     581            zn0 = ((((ALP050*zt   & 
     582               &   + ALP140*zs+ALP040)*zt   & 
     583               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     584               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     585               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     586               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     587               ! 
     588            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     589            ! 
     590            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 
     591            ! 
     592            ! beta 
     593            zn3 = BET003 
     594            ! 
     595            zn2 = BET012*zt + BET102*zs+BET002 
     596            ! 
     597            zn1 = ((BET031*zt   & 
     598               &   + BET121*zs+BET021)*zt   & 
     599               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     600               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     601               ! 
     602            zn0 = ((((BET050*zt   & 
     603               &   + BET140*zs+BET040)*zt   & 
     604               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     605               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     606               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     607               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     608               ! 
     609            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     610            ! 
     611            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 
     612            ! 
     613         END_3D 
    646614         ! 
    647615      CASE( np_seos )                  !==  simplified EOS  ==! 
    648616         ! 
    649          DO jk = 1, jpkm1 
    650             DO jj = 1, jpj 
    651                DO ji = 1, jpi 
    652                   zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    653                   zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    654                   zh  = gdept_n(ji,jj,jk)                ! depth in meters at t-point 
    655                   ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
    656                   ! 
    657                   zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    658                   pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha 
    659                   ! 
    660                   zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    661                   pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta 
    662                   ! 
    663                END DO 
    664             END DO 
    665          END DO 
     617         DO_3D_11_11( 1, jpkm1 ) 
     618            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     619            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     620            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
     621            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
     622            ! 
     623            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     624            pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha 
     625            ! 
     626            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     627            pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta 
     628            ! 
     629         END_3D 
    666630         ! 
    667631      CASE DEFAULT 
    … …  
    671635      END SELECT 
    672636      ! 
    673       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
    674          &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
     637      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', & 
     638         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk ) 
    675639      ! 
    676640      IF( ln_timing )   CALL timing_stop('rab_3d') 
    … …  
    679643 
    680644 
    681    SUBROUTINE rab_2d( pts, pdep, pab ) 
     645   SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 
    682646      !!---------------------------------------------------------------------- 
    683647      !!                 ***  ROUTINE rab_2d  *** 
    … …  
    687651      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    688652      !!---------------------------------------------------------------------- 
     653      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    689654      REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    690655      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m] 
    … …  
    704669      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    705670         ! 
    706          DO jj = 1, jpjm1 
    707             DO ji = 1, fs_jpim1   ! vector opt. 
    708                ! 
    709                zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    710                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    711                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    712                ! 
    713                ! alpha 
    714                zn3 = ALP003 
    715                ! 
    716                zn2 = ALP012*zt + ALP102*zs+ALP002 
    717                ! 
    718                zn1 = ((ALP031*zt   & 
    719                   &   + ALP121*zs+ALP021)*zt   & 
    720                   &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
    721                   &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
    722                   ! 
    723                zn0 = ((((ALP050*zt   & 
    724                   &   + ALP140*zs+ALP040)*zt   & 
    725                   &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
    726                   &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
    727                   &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
    728                   &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
    729                   ! 
    730                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    731                ! 
    732                pab(ji,jj,jp_tem) = zn * r1_rau0 
    733                ! 
    734                ! beta 
    735                zn3 = BET003 
    736                ! 
    737                zn2 = BET012*zt + BET102*zs+BET002 
    738                ! 
    739                zn1 = ((BET031*zt   & 
    740                   &   + BET121*zs+BET021)*zt   & 
    741                   &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
    742                   &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
    743                   ! 
    744                zn0 = ((((BET050*zt   & 
    745                   &   + BET140*zs+BET040)*zt   & 
    746                   &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
    747                   &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
    748                   &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
    749                   &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
    750                   ! 
    751                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    752                ! 
    753                pab(ji,jj,jp_sal) = zn / zs * r1_rau0 
    754                ! 
    755                ! 
    756             END DO 
    757          END DO 
    758          !                            ! Lateral boundary conditions 
    759          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     671         DO_2D_11_11 
     672            ! 
     673            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     674            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     675            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     676            ! 
     677            ! alpha 
     678            zn3 = ALP003 
     679            ! 
     680            zn2 = ALP012*zt + ALP102*zs+ALP002 
     681            ! 
     682            zn1 = ((ALP031*zt   & 
     683               &   + ALP121*zs+ALP021)*zt   & 
     684               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     685               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     686               ! 
     687            zn0 = ((((ALP050*zt   & 
     688               &   + ALP140*zs+ALP040)*zt   & 
     689               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     690               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     691               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     692               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     693               ! 
     694            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     695            ! 
     696            pab(ji,jj,jp_tem) = zn * r1_rau0 
     697            ! 
     698            ! beta 
     699            zn3 = BET003 
     700            ! 
     701            zn2 = BET012*zt + BET102*zs+BET002 
     702            ! 
     703            zn1 = ((BET031*zt   & 
     704               &   + BET121*zs+BET021)*zt   & 
     705               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     706               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     707               ! 
     708            zn0 = ((((BET050*zt   & 
     709               &   + BET140*zs+BET040)*zt   & 
     710               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     711               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     712               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     713               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     714               ! 
     715            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     716            ! 
     717            pab(ji,jj,jp_sal) = zn / zs * r1_rau0 
     718            ! 
     719            ! 
     720         END_2D 
    760721         ! 
    761722      CASE( np_seos )                  !==  simplified EOS  ==! 
    762723         ! 
    763          DO jj = 1, jpjm1 
    764             DO ji = 1, fs_jpim1   ! vector opt. 
    765                ! 
    766                zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
    767                zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
    768                zh    = pdep (ji,jj)                   ! depth at the partial step level 
    769                ! 
    770                zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    771                pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha 
    772                ! 
    773                zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    774                pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta 
    775                ! 
    776             END DO 
    777          END DO 
    778          !                            ! Lateral boundary conditions 
    779          CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                     
     724         DO_2D_11_11 
     725            ! 
     726            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     727            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     728            zh    = pdep (ji,jj)                   ! depth at the partial step level 
     729            ! 
     730            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     731            pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha 
     732            ! 
     733            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     734            pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta 
     735            ! 
     736         END_2D 
    780737         ! 
    781738      CASE DEFAULT 
    … …  
    785742      END SELECT 
    786743      ! 
    787       IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
    788          &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
     744      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', & 
     745         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' ) 
    789746      ! 
    790747      IF( ln_timing )   CALL timing_stop('rab_2d') 
    … …  
    793750 
    794751 
    795    SUBROUTINE rab_0d( pts, pdep, pab ) 
     752   SUBROUTINE rab_0d( pts, pdep, pab, Kmm ) 
    796753      !!---------------------------------------------------------------------- 
    797754      !!                 ***  ROUTINE rab_0d  *** 
    … …  
    801758      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    802759      !!---------------------------------------------------------------------- 
     760      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    803761      REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    804762      REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m] 
    … …  
    889847 
    890848 
    891    SUBROUTINE bn2( pts, pab, pn2 ) 
     849   SUBROUTINE bn2( pts, pab, pn2, Kmm ) 
    892850      !!---------------------------------------------------------------------- 
    893851      !!                  ***  ROUTINE bn2  *** 
    … …  
    903861      !! 
    904862      !!---------------------------------------------------------------------- 
     863      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    905864      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
    906865      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1] 
    … …  
    913872      IF( ln_timing )   CALL timing_start('bn2') 
    914873      ! 
    915       DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 ) 
    916          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90 
    917             DO ji = 1, jpi 
    918                zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   & 
    919                   &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )  
    920                   ! 
    921                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
    922                zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
    923                ! 
    924                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
    925                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
    926                   &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    927             END DO 
    928          END DO 
    929       END DO 
    930       ! 
    931       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
     874      DO_3D_11_11( 2, jpkm1 ) 
     875         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   & 
     876            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )  
     877            ! 
     878         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw  
     879         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw 
     880         ! 
     881         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     & 
     882            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  & 
     883            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 
     884      END_3D 
     885      ! 
     886      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk ) 
    932887      ! 
    933888      IF( ln_timing )   CALL timing_stop('bn2') 
    … …  
    965920      z1_T0   = 1._wp/40._wp 
    966921      ! 
    967       DO jj = 1, jpj 
    968          DO ji = 1, jpi 
    969             ! 
    970             zt  = ctmp   (ji,jj) * z1_T0 
    971             zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
    972             ztm = tmask(ji,jj,1) 
    973             ! 
    974             zn = ((((-2.1385727895e-01_wp*zt   & 
    975                &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
    976                &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
    977                &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
    978                &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
    979                &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
    980                &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
    981                &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
    982                ! 
    983             zd = (2.0035003456_wp*zt   & 
    984                &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
    985                &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
    986                ! 
    987             ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
    988                ! 
    989          END DO 
    990       END DO 
     922      DO_2D_11_11 
     923         ! 
     924         zt  = ctmp   (ji,jj) * z1_T0 
     925         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 ) 
     926         ztm = tmask(ji,jj,1) 
     927         ! 
     928         zn = ((((-2.1385727895e-01_wp*zt   & 
     929            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   & 
     930            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   & 
     931            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   & 
     932            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   & 
     933            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   & 
     934            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   & 
     935            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp 
     936            ! 
     937         zd = (2.0035003456_wp*zt   & 
     938            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   & 
     939            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp 
     940            ! 
     941         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm 
     942            ! 
     943      END_2D 
    991944      ! 
    992945      IF( ln_timing )   CALL timing_stop('eos_pt_from_ct') 
    … …  
    1020973         ! 
    1021974         z1_S0 = 1._wp / 35.16504_wp 
    1022          DO jj = 1, jpj 
    1023             DO ji = 1, jpi 
    1024                zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
    1025                ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
    1026                   &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
    1027             END DO 
    1028          END DO 
     975         DO_2D_11_11 
     976            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity 
     977            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs & 
     978               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp 
     979         END_2D 
    1029980         ptf(:,:) = ptf(:,:) * psal(:,:) 
    1030981         ! 
    … …  
    10931044 
    10941045 
    1095    SUBROUTINE eos_pen( pts, pab_pe, ppen ) 
     1046   SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm ) 
    10961047      !!---------------------------------------------------------------------- 
    10971048      !!                 ***  ROUTINE eos_pen  *** 
    … …  
    11131064      !!                    pab_pe(:,:,:,jp_sal) is beta_pe 
    11141065      !!---------------------------------------------------------------------- 
     1066      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index 
    11151067      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity 
    11161068      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe 
    … …  
    11281080      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    11291081         ! 
    1130          DO jk = 1, jpkm1 
    1131             DO jj = 1, jpj 
    1132                DO ji = 1, jpi 
    1133                   ! 
    1134                   zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth 
    1135                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    1136                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    1137                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    1138                   ! 
    1139                   ! potential energy non-linear anomaly 
    1140                   zn2 = (PEN012)*zt   & 
    1141                      &   + PEN102*zs+PEN002 
    1142                      ! 
    1143                   zn1 = ((PEN021)*zt   & 
    1144                      &   + PEN111*zs+PEN011)*zt   & 
    1145                      &   + (PEN201*zs+PEN101)*zs+PEN001 
    1146                      ! 
    1147                   zn0 = ((((PEN040)*zt   & 
    1148                      &   + PEN130*zs+PEN030)*zt   & 
    1149                      &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
    1150                      &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
    1151                      &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
    1152                      ! 
    1153                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1154                   ! 
    1155                   ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm 
    1156                   ! 
    1157                   ! alphaPE non-linear anomaly 
    1158                   zn2 = APE002 
    1159                   ! 
    1160                   zn1 = (APE011)*zt   & 
    1161                      &   + APE101*zs+APE001 
    1162                      ! 
    1163                   zn0 = (((APE030)*zt   & 
    1164                      &   + APE120*zs+APE020)*zt   & 
    1165                      &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
    1166                      &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
    1167                      ! 
    1168                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1169                   !                               
    1170                   pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 
    1171                   ! 
    1172                   ! betaPE non-linear anomaly 
    1173                   zn2 = BPE002 
    1174                   ! 
    1175                   zn1 = (BPE011)*zt   & 
    1176                      &   + BPE101*zs+BPE001 
    1177                      ! 
    1178                   zn0 = (((BPE030)*zt   & 
    1179                      &   + BPE120*zs+BPE020)*zt   & 
    1180                      &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
    1181                      &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
    1182                      ! 
    1183                   zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    1184                   !                               
    1185                   pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 
    1186                   ! 
    1187                END DO 
    1188             END DO 
    1189          END DO 
     1082         DO_3D_11_11( 1, jpkm1 ) 
     1083            ! 
     1084            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     1085            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     1086            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1087            ztm = tmask(ji,jj,jk)                                         ! tmask 
     1088            ! 
     1089            ! potential energy non-linear anomaly 
     1090            zn2 = (PEN012)*zt   & 
     1091               &   + PEN102*zs+PEN002 
     1092               ! 
     1093            zn1 = ((PEN021)*zt   & 
     1094               &   + PEN111*zs+PEN011)*zt   & 
     1095               &   + (PEN201*zs+PEN101)*zs+PEN001 
     1096               ! 
     1097            zn0 = ((((PEN040)*zt   & 
     1098               &   + PEN130*zs+PEN030)*zt   & 
     1099               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   & 
     1100               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   & 
     1101               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000 
     1102               ! 
     1103            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1104            ! 
     1105            ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm 
     1106            ! 
     1107            ! alphaPE non-linear anomaly 
     1108            zn2 = APE002 
     1109            ! 
     1110            zn1 = (APE011)*zt   & 
     1111               &   + APE101*zs+APE001 
     1112               ! 
     1113            zn0 = (((APE030)*zt   & 
     1114               &   + APE120*zs+APE020)*zt   & 
     1115               &   + (APE210*zs+APE110)*zs+APE010)*zt   & 
     1116               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000 
     1117               ! 
     1118            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1119            !                               
     1120            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 
     1121            ! 
     1122            ! betaPE non-linear anomaly 
     1123            zn2 = BPE002 
     1124            ! 
     1125            zn1 = (BPE011)*zt   & 
     1126               &   + BPE101*zs+BPE001 
     1127               ! 
     1128            zn0 = (((BPE030)*zt   & 
     1129               &   + BPE120*zs+BPE020)*zt   & 
     1130               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   & 
     1131               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000 
     1132               ! 
     1133            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
     1134            !                               
     1135            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 
     1136            ! 
     1137         END_3D 
    11901138         ! 
    11911139      CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==! 
    11921140         ! 
    1193          DO jk = 1, jpkm1 
    1194             DO jj = 1, jpj 
    1195                DO ji = 1, jpi 
    1196                   zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
    1197                   zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
    1198                   zh  = gdept_n(ji,jj,jk)              ! depth in meters  at t-point 
    1199                   ztm = tmask(ji,jj,jk)                ! tmask 
    1200                   zn  = 0.5_wp * zh * r1_rau0 * ztm 
    1201                   !                                    ! Potential Energy 
    1202                   ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
    1203                   !                                    ! alphaPE 
    1204                   pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
    1205                   pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
    1206                   ! 
    1207                END DO 
    1208             END DO 
    1209          END DO 
     1141         DO_3D_11_11( 1, jpkm1 ) 
     1142            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0) 
     1143            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0) 
     1144            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
     1145            ztm = tmask(ji,jj,jk)                ! tmask 
     1146            zn  = 0.5_wp * zh * r1_rau0 * ztm 
     1147            !                                    ! Potential Energy 
     1148            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
     1149            !                                    ! alphaPE 
     1150            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn 
     1151            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn 
     1152            ! 
     1153         END_3D 
    12101154         ! 
    12111155      CASE DEFAULT 
    … …  
    12351179      !!---------------------------------------------------------------------- 
    12361180      ! 
    1237       REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state 
    12381181      READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 ) 
    12391182901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' ) 
    12401183      ! 
    1241       REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state 
    12421184      READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 ) 
    12431185902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' ) 
Note: See TracChangeset for help on using the changeset viewer.