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 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/eosbn2.F90 – NEMO

Ignore:
Timestamp:
2020-05-14T21:46:00+02:00 (4 years ago)
Author:
smueller
Message:

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

        old new  
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@HEAD         sette 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/eosbn2.F90

    r12178 r12928  
    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) 
     
    192191      !!                   ***  ROUTINE eos_insitu  *** 
    193192      !! 
    194       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
     193      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
    195194      !!       potential temperature and salinity using an equation of state 
    196195      !!       selected in the nameos namelist 
    197196      !! 
    198       !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
     197      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 
    199198      !!         with   prd    in situ density anomaly      no units 
    200199      !!                t      TEOS10: CT or EOS80: PT      Celsius 
     
    202201      !!                z      depth                        meters 
    203202      !!                rho    in situ density              kg/m^3 
    204       !!                rau0   reference density            kg/m^3 
     203      !!                rho0   reference density            kg/m^3 
    205204      !! 
    206205      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     
    211210      !! 
    212211      !!     ln_seos : simplified equation of state 
    213       !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rau0 
     212      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 
    214213      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
    215214      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     
    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_rho0 - 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_rho0 * 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/rho0) 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_rho0 - 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_rho0 - 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) = ( rho0 + 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_rho0 * 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/rho0) 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_rho0 - 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_rho0               ! 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_rho0 * 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_rho0 * 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_rho0 * 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_rho0 * 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_rho0 
     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_rho0 
     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_rho0   ! 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_rho0   ! 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] 
     
    841799         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    842800         ! 
    843          pab(jp_tem) = zn * r1_rau0 
     801         pab(jp_tem) = zn * r1_rho0 
    844802         ! 
    845803         ! beta 
     
    862820         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    863821         ! 
    864          pab(jp_sal) = zn / zs * r1_rau0 
     822         pab(jp_sal) = zn / zs * r1_rho0 
    865823         ! 
    866824         ! 
     
    873831         ! 
    874832         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    875          pab(jp_tem) = zn * r1_rau0   ! alpha 
     833         pab(jp_tem) = zn * r1_rho0   ! alpha 
    876834         ! 
    877835         zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    878          pab(jp_sal) = zn * r1_rau0   ! beta 
     836         pab(jp_sal) = zn * r1_rho0   ! beta 
    879837         ! 
    880838      CASE DEFAULT 
     
    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  *** 
     
    11011052      !! ** Method  :   PE is defined analytically as the vertical  
    11021053      !!                   primitive of EOS times -g integrated between 0 and z>0. 
    1103       !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd 
     1054      !!                pen is the nonlinear bsq-PE anomaly: pen = ( PE - rho0 gz ) / rho0 gz - rd 
    11041055      !!                                                      = 1/z * /int_0^z rd dz - rd  
    11051056      !!                                where rd is the density anomaly (see eos_rhd function) 
    11061057      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S: 
    1107       !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT 
    1108       !!                    ab_pe(2) =   1/(rau0 gz) * dPE/dS + drd/dS =   d(pen)/dS 
     1058      !!                    ab_pe(1) = - 1/(rho0 gz) * dPE/dT + drd/dT = - d(pen)/dT 
     1059      !!                    ab_pe(2) =   1/(rho0 gz) * dPE/dS + drd/dS =   d(pen)/dS 
    11091060      !! 
    11101061      !! ** Action  : - pen         : PE anomaly given at T-points 
     
    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_rho0 * 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_rho0 * 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_rho0 * 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_rho0 * 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' ) 
    12441186      IF(lwm) WRITE( numond, nameos ) 
    12451187      ! 
    1246       rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     1188      rho0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
    12471189      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
    12481190      ! 
     
    16561598            WRITE(numout,*) '   ==>>>   use of simplified eos:    ' 
    16571599            WRITE(numout,*) '              rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 
    1658             WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rau0' 
     1600            WRITE(numout,*) '                                       + b0*(1+lambda2/2*dT+mu2*Z)*dS - nu*dT*dS] / rho0' 
    16591601            WRITE(numout,*) '              with the following coefficients :' 
    16601602            WRITE(numout,*) '                 thermal exp. coef.    rn_a0      = ', rn_a0 
     
    16751617      END SELECT 
    16761618      ! 
    1677       rau0_rcp    = rau0 * rcp  
    1678       r1_rau0     = 1._wp / rau0 
     1619      rho0_rcp    = rho0 * rcp  
     1620      r1_rho0     = 1._wp / rho0 
    16791621      r1_rcp      = 1._wp / rcp 
    1680       r1_rau0_rcp = 1._wp / rau0_rcp  
     1622      r1_rho0_rcp = 1._wp / rho0_rcp  
    16811623      ! 
    16821624      IF(lwp) THEN 
     
    16931635      IF(lwp) WRITE(numout,*) 
    16941636      IF(lwp) WRITE(numout,*) '   Associated physical constant' 
    1695       IF(lwp) WRITE(numout,*) '      volumic mass of reference           rau0  = ', rau0   , ' kg/m^3' 
    1696       IF(lwp) WRITE(numout,*) '      1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg' 
     1637      IF(lwp) WRITE(numout,*) '      volumic mass of reference           rho0  = ', rho0   , ' kg/m^3' 
     1638      IF(lwp) WRITE(numout,*) '      1. / rho0                        r1_rho0  = ', r1_rho0, ' m^3/kg' 
    16971639      IF(lwp) WRITE(numout,*) '      ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
    1698       IF(lwp) WRITE(numout,*) '      rau0 * rcp                       rau0_rcp = ', rau0_rcp 
    1699       IF(lwp) WRITE(numout,*) '      1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
     1640      IF(lwp) WRITE(numout,*) '      rho0 * rcp                       rho0_rcp = ', rho0_rcp 
     1641      IF(lwp) WRITE(numout,*) '      1. / ( rho0 * rcp )           r1_rho0_rcp = ', r1_rho0_rcp 
    17001642      ! 
    17011643   END SUBROUTINE eos_init 
Note: See TracChangeset for help on using the changeset viewer.