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 15066 for NEMO – NEMO

Changeset 15066 for NEMO


Ignore:
Timestamp:
2021-07-01T13:17:50+02:00 (3 years ago)
Author:
sparonuz
Message:

Added interface in eosbn2.F90 to ease the weight of casts necessary in mixed precision + Removed unnecessary casts

Location:
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/C1D/step_c1d.F90

    r14644 r15066  
    7272                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points 
    7373                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points 
    74                          CALL bn2( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
    75                          CALL bn2( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency 
     74                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     75                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency 
    7676 
    7777      !  VERTICAL PHYSICS 
     
    108108      IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes 
    109109                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing 
    110                         CALL eos( CASTEWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl 
     110                        CALL eos( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl 
    111111      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s) 
    112112                        CALL tra_atf( kstp, Nbb, Nnn, Naa, ts )                 ! time filtering of "now" tracer arrays 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DOM/domqco.F90

    r14986 r15066  
    118118      !                                ! Horizontal interpolation of e3t 
    119119#if defined key_RK3 
    120       CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) ) 
    121       CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           ) 
     120      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) ) 
     121      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           ) 
    122122#else 
    123       CALL dom_qco_r3c( CASTWP(ssh(:,:,Kbb)), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           ) 
    124       CALL dom_qco_r3c( CASTWP(ssh(:,:,Kmm)), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 
     123      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           ) 
     124      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 
    125125#endif 
    126126      ! dom_qco_r3c defines over [nn_hls, nn_hls-1, nn_hls, nn_hls-1] 
     
    142142      !!              - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional) 
    143143      !!---------------------------------------------------------------------- 
    144       REAL(wp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m] 
     144      REAL(dp), DIMENSION(:,:)          , INTENT(in   )  ::   pssh               ! sea surface height   [m] 
    145145      REAL(wp), DIMENSION(:,:)          , INTENT(  out)  ::   pr3t, pr3u, pr3v   ! ssh/h0 ratio at t-, u-, v-,points  [-] 
    146146      REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT(  out)  ::   pr3f               ! ssh/h0 ratio at f-point   [-] 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynspg_ts.F90

    r14986 r15066  
    8787 
    8888 
    89    INTERFACE dyn_drg_init 
    90       MODULE PROCEDURE dyn_drg_init_wp, dyn_drg_init_mixed 
    91    END INTERFACE 
    92  
    9389   !! * Substitutions 
    9490#  include "do_loop_substitute.h90" 
     
    288284         END_2D 
    289285      ELSE                !* remove baroclinic drag AND provide the barotropic drag coefficients 
    290          CALL dyn_drg_init( Kbb, Kmm, CASTWP(puu), CASTWP(pvv), puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 
     286         CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 
    291287      ENDIF 
    292288      ! 
     
    13471343      
    13481344 
    1349    SUBROUTINE dyn_drg_init_wp( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
     1345   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
    13501346      !!---------------------------------------------------------------------- 
    13511347      !!                  ***  ROUTINE dyn_drg_init  *** 
     
    13581354      !!---------------------------------------------------------------------- 
    13591355      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices 
    1360       REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation 
     1356      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation 
    13611357      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels 
    13621358      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS 
     
    14511447      ENDIF 
    14521448      ! 
    1453    END SUBROUTINE dyn_drg_init_wp 
    1454  
    1455    SUBROUTINE dyn_drg_init_mixed( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 
    1456       !!---------------------------------------------------------------------- 
    1457       !!                  ***  ROUTINE dyn_drg_init  *** 
    1458       !!                     
    1459       !! ** Purpose : - add the baroclinic top/bottom drag contribution to  
    1460       !!              the baroclinic part of the barotropic RHS 
    1461       !!              - compute the barotropic drag coefficients 
    1462       !! 
    1463       !! ** Method  :   computation done over the INNER domain only  
    1464       !!---------------------------------------------------------------------- 
    1465       INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices 
    1466       REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation 
    1467       REAL(sp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels 
    1468       REAL(sp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS 
    1469       REAL(sp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients 
    1470       ! 
    1471       INTEGER  ::   ji, jj   ! dummy loop indices 
    1472       INTEGER  ::   ikbu, ikbv, iktu, iktv 
    1473       REAL(sp) ::   zztmp 
    1474       REAL(sp), DIMENSION(jpi,jpj) ::   zu_i, zv_i 
    1475       !!---------------------------------------------------------------------- 
    1476       ! 
    1477       !                    !==  Set the barotropic drag coef.  ==! 
    1478       ! 
    1479       IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities) 
    1480           
    1481          DO_2D( 0, 0, 0, 0 ) 
    1482             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 
    1483             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 
    1484          END_2D 
    1485       ELSE                          ! bottom friction only 
    1486          DO_2D( 0, 0, 0, 0 ) 
    1487             pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 
    1488             pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 
    1489          END_2D 
    1490       ENDIF 
    1491       ! 
    1492       !                    !==  BOTTOM stress contribution from baroclinic velocities  ==! 
    1493       ! 
    1494       IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities 
    1495           
    1496          DO_2D( 0, 0, 0, 0 ) 
    1497             ikbu = mbku(ji,jj)        
    1498             ikbv = mbkv(ji,jj)     
    1499             zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 
    1500             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 
    1501          END_2D 
    1502       ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities 
    1503           
    1504          DO_2D( 0, 0, 0, 0 ) 
    1505             ikbu = mbku(ji,jj)        
    1506             ikbv = mbkv(ji,jj)     
    1507             zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 
    1508             zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 
    1509          END_2D 
    1510       ENDIF 
    1511       ! 
    1512       IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please ! 
    1513          zztmp = -1._wp / rDt_e 
    1514          DO_2D( 0, 0, 0, 0 ) 
    1515             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 &  
    1516                  &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  ) 
    1517             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 &  
    1518                  &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  ) 
    1519          END_2D 
    1520       ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 
    1521           
    1522          DO_2D( 0, 0, 0, 0 ) 
    1523             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 
    1524             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 
    1525          END_2D 
    1526       END IF 
    1527       ! 
    1528       !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case) 
    1529       ! 
    1530       IF( ln_isfcav.OR.ln_drgice_imp ) THEN 
    1531          ! 
    1532          IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity 
    1533              
    1534             DO_2D( 0, 0, 0, 0 ) 
    1535                iktu = miku(ji,jj) 
    1536                iktv = mikv(ji,jj) 
    1537                zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) 
    1538                zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 
    1539             END_2D 
    1540          ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity 
    1541              
    1542             DO_2D( 0, 0, 0, 0 ) 
    1543                iktu = miku(ji,jj) 
    1544                iktv = mikv(ji,jj) 
    1545                zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) 
    1546                zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 
    1547             END_2D 
    1548          ENDIF 
    1549          ! 
    1550          !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 
    1551           
    1552          DO_2D( 0, 0, 0, 0 ) 
    1553             pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 
    1554             pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 
    1555          END_2D 
    1556          ! 
    1557       ENDIF 
    1558       ! 
    1559    END SUBROUTINE dyn_drg_init_mixed 
     1449   END SUBROUTINE dyn_drg_init 
    15601450 
    15611451   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90

    r14986 r15066  
    5656   !                  !! * Interface 
    5757   INTERFACE eos 
    58       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d 
     58      MODULE PROCEDURE eos_insitu_wp, eos_insitu_pot_wp, eos_insitu_2d, eos_insitu_pot_2d 
     59      MODULE PROCEDURE eos_insitu_mixed, eos_insitu_pot_mixed 
    5960   END INTERFACE 
    6061   ! 
    6162   INTERFACE eos_rab 
    62       MODULE PROCEDURE rab_3d, rab_2d, rab_0d 
     63      MODULE PROCEDURE rab_3d_wp, rab_2d, rab_0d 
     64      MODULE PROCEDURE rab_3d_mixed 
    6365   END INTERFACE 
    6466   ! 
     
    190192CONTAINS 
    191193 
    192    SUBROUTINE eos_insitu( pts, prd, pdep ) 
     194   SUBROUTINE eos_insitu_wp( pts, prd, pdep ) 
    193195      !! 
    194196      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     
    197199      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m] 
    198200      !! 
    199       CALL eos_insitu_t( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 
    200    END SUBROUTINE eos_insitu 
    201  
    202    SUBROUTINE eos_insitu_t( pts, ktts, prd, ktrd, pdep, ktdep ) 
    203       !!---------------------------------------------------------------------- 
    204       !!                   ***  ROUTINE eos_insitu  *** 
     201      CALL eos_insitu_t_wp( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 
     202   END SUBROUTINE eos_insitu_wp 
     203 
     204   SUBROUTINE eos_insitu_mixed( pts, prd, pdep ) 
     205      !! 
     206      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     207      !                                                      ! 2 : salinity               [psu] 
     208      REAL(sp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd   ! in situ density            [-] 
     209      REAL(sp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep  ! depth                      [m] 
     210      !! 
     211      CALL eos_insitu_t_mixed( pts, is_tile(pts), prd, is_tile(prd), pdep, is_tile(pdep) ) 
     212   END SUBROUTINE  
     213 
     214   SUBROUTINE eos_insitu_t_wp( pts, ktts, prd, ktrd, pdep, ktdep ) 
     215      !!---------------------------------------------------------------------- 
     216      !!                   ***  ROUTINE  
    205217      !! 
    206218      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
     
    306318      IF( ln_timing )   CALL timing_stop('eos-insitu') 
    307319      ! 
    308    END SUBROUTINE eos_insitu_t 
    309  
    310  
    311    SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep ) 
     320   END SUBROUTINE eos_insitu_t_wp 
     321   SUBROUTINE eos_insitu_t_mixed( pts, ktts, prd, ktrd, pdep, ktdep ) 
     322      !!---------------------------------------------------------------------- 
     323      !!                   ***  ROUTINE eos_insitu  *** 
     324      !! 
     325      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
     326      !!       potential temperature and salinity using an equation of state 
     327      !!       selected in the nameos namelist 
     328      !! 
     329      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 
     330      !!         with   prd    in situ density anomaly      no units 
     331      !!                t      TEOS10: CT or EOS80: PT      Celsius 
     332      !!                s      TEOS10: SA or EOS80: SP      TEOS10: g/kg or EOS80: psu 
     333      !!                z      depth                        meters 
     334      !!                rho    in situ density              kg/m^3 
     335      !!                rho0   reference density            kg/m^3 
     336      !! 
     337      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     338      !!         Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg 
     339      !! 
     340      !!     ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z). 
     341      !!         Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu 
     342      !! 
     343      !!     ln_seos : simplified equation of state 
     344      !!              prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 
     345      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
     346      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     347      !!              Vallis like equation: use default values of coefficients 
     348      !! 
     349      !! ** Action  :   compute prd , the in situ density (no units) 
     350      !! 
     351      !! References :   Roquet et al, Ocean Modelling, in preparation (2014) 
     352      !!                Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 
     353      !!                TEOS-10 Manual, 2010 
     354      !!---------------------------------------------------------------------- 
     355      INTEGER                                 , INTENT(in   ) ::   ktts, ktrd, ktdep 
     356      REAL(dp), DIMENSION(A2D_T(ktts) ,JPK,JPTS), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius] 
     357      !                                                                  ! 2 : salinity               [psu] 
     358      REAL(sp), DIMENSION(A2D_T(ktrd) ,JPK     ), INTENT(  out) ::   prd   ! in situ density            [-] 
     359      REAL(sp), DIMENSION(A2D_T(ktdep),JPK     ), INTENT(in   ) ::   pdep  ! depth                      [m] 
     360      ! 
     361      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     362      REAL(sp) ::   zt , zh , zs , ztm        ! local scalars 
     363      REAL(sp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     364      !!---------------------------------------------------------------------- 
     365      ! 
     366      IF( ln_timing )   CALL timing_start('eos-insitu') 
     367      ! 
     368      SELECT CASE( neos ) 
     369      ! 
     370      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     371         ! 
     372         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     373            ! 
     374            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     375            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     376            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     377            ztm = tmask(ji,jj,jk)                                         ! tmask 
     378            ! 
     379            zn3 = EOS013*zt   & 
     380               &   + EOS103*zs+EOS003 
     381               ! 
     382            zn2 = (EOS022*zt   & 
     383               &   + EOS112*zs+EOS012)*zt   & 
     384               &   + (EOS202*zs+EOS102)*zs+EOS002 
     385               ! 
     386            zn1 = (((EOS041*zt   & 
     387               &   + EOS131*zs+EOS031)*zt   & 
     388               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     389               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     390               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     391               ! 
     392            zn0 = (((((EOS060*zt   & 
     393               &   + EOS150*zs+EOS050)*zt   & 
     394               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     395               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     396               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     397               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     398               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     399               ! 
     400            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     401            ! 
     402            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     403            ! 
     404         END_3D 
     405         ! 
     406      CASE( np_seos )                !==  simplified EOS  ==! 
     407         ! 
     408         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     409            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     410            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     411            zh  = pdep (ji,jj,jk) 
     412            ztm = tmask(ji,jj,jk) 
     413            ! 
     414            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   & 
     415               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   & 
     416               &  - rn_nu * zt * zs 
     417               ! 
     418            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
     419         END_3D 
     420         ! 
     421      END SELECT 
     422      ! 
     423      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-insitu  : ', kdim=jpk ) 
     424      ! 
     425      IF( ln_timing )   CALL timing_stop('eos-insitu') 
     426      ! 
     427   END SUBROUTINE eos_insitu_t_mixed 
     428 
     429 
     430   SUBROUTINE eos_insitu_pot_wp( pts, prd, prhop, pdep ) 
    312431      !! 
    313432      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     
    317436      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep   ! depth                      [m] 
    318437      !! 
    319       CALL eos_insitu_pot_t( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 
    320    END SUBROUTINE eos_insitu_pot 
    321  
    322  
    323    SUBROUTINE eos_insitu_pot_t( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 
     438      CALL eos_insitu_pot_t_wp( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 
     439   END SUBROUTINE eos_insitu_pot_wp 
     440 
     441   SUBROUTINE eos_insitu_pot_mixed( pts, prd, prhop, pdep ) 
     442      !! 
     443      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     444      !                                                       ! 2 : salinity               [psu] 
     445      REAL(sp), DIMENSION(:,:,:)  , INTENT(  out) ::   prd    ! in situ density            [-] 
     446      REAL(dp), DIMENSION(:,:,:)  , INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     447      REAL(sp), DIMENSION(:,:,:)  , INTENT(in   ) ::   pdep   ! depth                      [m] 
     448      !! 
     449      CALL eos_insitu_pot_t_mixed( pts, is_tile(pts), prd, is_tile(prd), prhop, is_tile(prhop), pdep, is_tile(pdep) ) 
     450   END SUBROUTINE eos_insitu_pot_mixed 
     451 
     452 
     453   SUBROUTINE eos_insitu_pot_t_wp( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 
    324454      !!---------------------------------------------------------------------- 
    325455      !!                  ***  ROUTINE eos_insitu_pot  *** 
     
    475605      IF( ln_timing )   CALL timing_stop('eos-pot') 
    476606      ! 
    477    END SUBROUTINE eos_insitu_pot_t 
     607   END SUBROUTINE eos_insitu_pot_t_wp 
     608 
     609   SUBROUTINE eos_insitu_pot_t_mixed( pts, ktts, prd, ktrd, prhop, ktrhop, pdep, ktdep ) 
     610      !!---------------------------------------------------------------------- 
     611      !!                  ***  ROUTINE eos_insitu_pot  *** 
     612      !! 
     613      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
     614      !!      potential volumic mass (Kg/m3) from potential temperature and 
     615      !!      salinity fields using an equation of state selected in the 
     616      !!     namelist. 
     617      !! 
     618      !! ** Action  : - prd  , the in situ density (no units) 
     619      !!              - prhop, the potential volumic mass (Kg/m3) 
     620      !! 
     621      !!---------------------------------------------------------------------- 
     622      INTEGER                                  , INTENT(in   ) ::   ktts, ktrd, ktrhop, ktdep 
     623      REAL(dp), DIMENSION(A2D_T(ktts)  ,JPK,JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     624      !                                                                    ! 2 : salinity               [psu] 
     625      REAL(sp), DIMENSION(A2D_T(ktrd)  ,JPK     ), INTENT(  out) ::   prd    ! in situ density            [-] 
     626      REAL(dp), DIMENSION(A2D_T(ktrhop),JPK     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     627      REAL(sp), DIMENSION(A2D_T(ktdep) ,JPK     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
     628      ! 
     629      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     630      INTEGER  ::   jdof 
     631      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     632      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     633      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
     634      !!---------------------------------------------------------------------- 
     635      ! 
     636      IF( ln_timing )   CALL timing_start('eos-pot') 
     637      ! 
     638      SELECT CASE ( neos ) 
     639      ! 
     640      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     641         ! 
     642         ! Stochastic equation of state 
     643         IF ( ln_sto_eos ) THEN 
     644            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     645            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     646            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     647            DO jsmp = 1, 2*nn_sto_eos, 2 
     648              zsign(jsmp)   = 1._wp 
     649              zsign(jsmp+1) = -1._wp 
     650            END DO 
     651            ! 
     652            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     653               ! 
     654               ! compute density (2*nn_sto_eos) times: 
     655               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     656               ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     657               DO jsmp = 1, nn_sto_eos*2 
     658                  jdof   = (jsmp + 1) / 2 
     659                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     660                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     661                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     662                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     663                  ztm    = tmask(ji,jj,jk)                                         ! tmask 
     664                  ! 
     665                  zn3 = EOS013*zt   & 
     666                     &   + EOS103*zs+EOS003 
     667                     ! 
     668                  zn2 = (EOS022*zt   & 
     669                     &   + EOS112*zs+EOS012)*zt   & 
     670                     &   + (EOS202*zs+EOS102)*zs+EOS002 
     671                     ! 
     672                  zn1 = (((EOS041*zt   & 
     673                     &   + EOS131*zs+EOS031)*zt   & 
     674                     &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     675                     &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     676                     &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     677                     ! 
     678                  zn0_sto(jsmp) = (((((EOS060*zt   & 
     679                     &   + EOS150*zs+EOS050)*zt   & 
     680                     &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     681                     &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     682                     &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     683                     &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     684                     &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     685                     ! 
     686                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     687               END DO 
     688               ! 
     689               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     690               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     691               DO jsmp = 1, nn_sto_eos*2 
     692                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
     693                  ! 
     694                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
     695               END DO 
     696               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     697               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     698            END_3D 
     699            DEALLOCATE(zn0_sto,zn_sto,zsign) 
     700         ! Non-stochastic equation of state 
     701         ELSE 
     702            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     703               ! 
     704               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     705               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     706               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     707               ztm = tmask(ji,jj,jk)                                         ! tmask 
     708               ! 
     709               zn3 = EOS013*zt   & 
     710                  &   + EOS103*zs+EOS003 
     711                  ! 
     712               zn2 = (EOS022*zt   & 
     713                  &   + EOS112*zs+EOS012)*zt   & 
     714                  &   + (EOS202*zs+EOS102)*zs+EOS002 
     715                  ! 
     716               zn1 = (((EOS041*zt   & 
     717                  &   + EOS131*zs+EOS031)*zt   & 
     718                  &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     719                  &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     720                  &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     721                  ! 
     722               zn0 = (((((EOS060*zt   & 
     723                  &   + EOS150*zs+EOS050)*zt   & 
     724                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     725                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     726                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     727                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     728                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     729                  ! 
     730               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     731               ! 
     732               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     733               ! 
     734               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     735            END_3D 
     736         ENDIF 
     737 
     738      CASE( np_seos )                !==  simplified EOS  ==! 
     739         ! 
     740         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     741            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp 
     742            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp 
     743            zh  = pdep (ji,jj,jk) 
     744            ztm = tmask(ji,jj,jk) 
     745            !                                                     ! potential density referenced at the surface 
     746            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   & 
     747               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
     748               &  - rn_nu * zt * zs 
     749            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
     750            !                                                     ! density anomaly (masked) 
     751            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
     752            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
     753            ! 
     754         END_3D 
     755         ! 
     756      END SELECT 
     757      ! 
     758      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(prd, wp), clinfo1=' eos-pot: ', & 
     759         &                                  tab3d_2=REAL(prhop, wp), clinfo2=' pot : ', kdim=jpk ) 
     760      ! 
     761      IF( ln_timing )   CALL timing_stop('eos-pot') 
     762      ! 
     763   END SUBROUTINE eos_insitu_pot_t_mixed 
    478764 
    479765 
     
    661947 
    662948 
    663    SUBROUTINE rab_3d( pts, pab, Kmm ) 
     949   SUBROUTINE rab_3d_wp( pts, pab, Kmm ) 
    664950      !! 
    665951      INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index 
     
    667953      REAL(wp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
    668954      !! 
    669       CALL rab_3d_t( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 
    670    END SUBROUTINE rab_3d 
    671  
    672  
    673    SUBROUTINE rab_3d_t( pts, ktts, pab, ktab, Kmm ) 
     955      CALL rab_3d_t_wp( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 
     956   END SUBROUTINE rab_3d_wp 
     957 
     958   SUBROUTINE rab_3d_mixed( pts, pab, Kmm ) 
     959      !! 
     960      INTEGER                     , INTENT(in   ) ::   Kmm   ! time level index 
     961      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
     962      REAL(sp), DIMENSION(:,:,:,:), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
     963      !! 
     964      CALL rab_3d_t_mixed( pts, is_tile(pts), pab, is_tile(pab), Kmm ) 
     965   END SUBROUTINE rab_3d_mixed 
     966 
     967 
     968   SUBROUTINE rab_3d_t_wp( pts, ktts, pab, ktab, Kmm ) 
    674969      !!---------------------------------------------------------------------- 
    675970      !!                 ***  ROUTINE rab_3d  *** 
     
    7751070      IF( ln_timing )   CALL timing_stop('rab_3d') 
    7761071      ! 
    777    END SUBROUTINE rab_3d_t 
    778  
    779  
    780    SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 
    781       !! 
    782       INTEGER                   , INTENT(in   ) ::   Kmm   ! time level index 
    783       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    784       REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pdep   ! depth                  [m] 
    785       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pab    ! thermal/haline expansion ratio 
    786       !! 
    787       CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm) 
    788    END SUBROUTINE rab_2d 
    789  
    790  
    791    SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm ) 
    792       !!---------------------------------------------------------------------- 
    793       !!                 ***  ROUTINE rab_2d  *** 
    794       !! 
    795       !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked) 
     1072   END SUBROUTINE rab_3d_t_wp 
     1073 
     1074   SUBROUTINE rab_3d_t_mixed( pts, ktts, pab, ktab, Kmm ) 
     1075      !!---------------------------------------------------------------------- 
     1076      !!                 ***  ROUTINE rab_3d  *** 
     1077      !! 
     1078      !! ** Purpose :   Calculates thermal/haline expansion ratio at T-points 
     1079      !! 
     1080      !! ** Method  :   calculates alpha / beta at T-points 
    7961081      !! 
    7971082      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
    7981083      !!---------------------------------------------------------------------- 
    799       INTEGER                            , INTENT(in   ) ::   Kmm   ! time level index 
    800       INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktab 
    801       REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! pot. temperature & salinity 
    802       REAL(wp), DIMENSION(A2D_T(ktdep)    ), INTENT(in   ) ::   pdep   ! depth                  [m] 
    803       REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT(  out) ::   pab    ! thermal/haline expansion ratio 
     1084      INTEGER                                , INTENT(in   ) ::   Kmm   ! time level index 
     1085      INTEGER                                , INTENT(in   ) ::   ktts, ktab 
     1086      REAL(dp), DIMENSION(A2D_T(ktts),JPK,JPTS), INTENT(in   ) ::   pts   ! pot. temperature & salinity 
     1087      REAL(sp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(  out) ::   pab   ! thermal/haline expansion ratio 
    8041088      ! 
    8051089      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    806       REAL(wp) ::   zt , zh , zs              ! local scalars 
     1090      REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
    8071091      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
    8081092      !!---------------------------------------------------------------------- 
    8091093      ! 
    810       IF( ln_timing )   CALL timing_start('rab_2d') 
    811       ! 
    812       pab(:,:,:) = 0._wp 
     1094      IF( ln_timing )   CALL timing_start('rab_3d') 
    8131095      ! 
    8141096      SELECT CASE ( neos ) 
     
    8161098      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    8171099         ! 
    818          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    819             ! 
    820             zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
    821             zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    822             zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1100         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     1101            ! 
     1102            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth 
     1103            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     1104            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1105            ztm = tmask(ji,jj,jk)                                         ! tmask 
    8231106            ! 
    8241107            ! alpha 
     
    8411124            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    8421125            ! 
    843             pab(ji,jj,jp_tem) = zn * r1_rho0 
     1126            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
    8441127            ! 
    8451128            ! beta 
     
    8621145            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    8631146            ! 
     1147            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
     1148            ! 
     1149         END_3D 
     1150         ! 
     1151      CASE( np_seos )                  !==  simplified EOS  ==! 
     1152         ! 
     1153         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     1154            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0) 
     1155            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0) 
     1156            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point 
     1157            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask 
     1158            ! 
     1159            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
     1160            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
     1161            ! 
     1162            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
     1163            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
     1164            ! 
     1165         END_3D 
     1166         ! 
     1167      CASE DEFAULT 
     1168         WRITE(ctmp1,*) '          bad flag value for neos = ', neos 
     1169         CALL ctl_stop( 'rab_3d:', ctmp1 ) 
     1170         ! 
     1171      END SELECT 
     1172      ! 
     1173      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=REAL(pab(:,:,:,jp_tem), wp), clinfo1=' rab_3d_t: ', & 
     1174         &                                  tab3d_2=REAL(pab(:,:,:,jp_sal),  wp), clinfo2=' rab_3d_s : ', kdim=jpk ) 
     1175      ! 
     1176      IF( ln_timing )   CALL timing_stop('rab_3d') 
     1177      ! 
     1178   END SUBROUTINE rab_3d_t_mixed 
     1179 
     1180 
     1181   SUBROUTINE rab_2d( pts, pdep, pab, Kmm ) 
     1182      !! 
     1183      INTEGER                   , INTENT(in   ) ::   Kmm   ! time level index 
     1184      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! pot. temperature & salinity 
     1185      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pdep   ! depth                  [m] 
     1186      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pab    ! thermal/haline expansion ratio 
     1187      !! 
     1188      CALL rab_2d_t(pts, is_tile(pts), pdep, is_tile(pdep), pab, is_tile(pab), Kmm) 
     1189   END SUBROUTINE rab_2d 
     1190 
     1191 
     1192   SUBROUTINE rab_2d_t( pts, ktts, pdep, ktdep, pab, ktab, Kmm ) 
     1193      !!---------------------------------------------------------------------- 
     1194      !!                 ***  ROUTINE rab_2d  *** 
     1195      !! 
     1196      !! ** Purpose :   Calculates thermal/haline expansion ratio for a 2d field (unmasked) 
     1197      !! 
     1198      !! ** Action  : - pab     : thermal/haline expansion ratio at T-points 
     1199      !!---------------------------------------------------------------------- 
     1200      INTEGER                            , INTENT(in   ) ::   Kmm   ! time level index 
     1201      INTEGER                            , INTENT(in   ) ::   ktts, ktdep, ktab 
     1202      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! pot. temperature & salinity 
     1203      REAL(wp), DIMENSION(A2D_T(ktdep)    ), INTENT(in   ) ::   pdep   ! depth                  [m] 
     1204      REAL(wp), DIMENSION(A2D_T(ktab),JPTS), INTENT(  out) ::   pab    ! thermal/haline expansion ratio 
     1205      ! 
     1206      INTEGER  ::   ji, jj, jk                ! dummy loop indices 
     1207      REAL(wp) ::   zt , zh , zs              ! local scalars 
     1208      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     1209      !!---------------------------------------------------------------------- 
     1210      ! 
     1211      IF( ln_timing )   CALL timing_start('rab_2d') 
     1212      ! 
     1213      pab(:,:,:) = 0._wp 
     1214      ! 
     1215      SELECT CASE ( neos ) 
     1216      ! 
     1217      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
     1218         ! 
     1219         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     1220            ! 
     1221            zh  = pdep(ji,jj) * r1_Z0                                  ! depth 
     1222            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     1223            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     1224            ! 
     1225            ! alpha 
     1226            zn3 = ALP003 
     1227            ! 
     1228            zn2 = ALP012*zt + ALP102*zs+ALP002 
     1229            ! 
     1230            zn1 = ((ALP031*zt   & 
     1231               &   + ALP121*zs+ALP021)*zt   & 
     1232               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   & 
     1233               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001 
     1234               ! 
     1235            zn0 = ((((ALP050*zt   & 
     1236               &   + ALP140*zs+ALP040)*zt   & 
     1237               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   & 
     1238               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   & 
     1239               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   & 
     1240               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000 
     1241               ! 
     1242            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     1243            ! 
     1244            pab(ji,jj,jp_tem) = zn * r1_rho0 
     1245            ! 
     1246            ! beta 
     1247            zn3 = BET003 
     1248            ! 
     1249            zn2 = BET012*zt + BET102*zs+BET002 
     1250            ! 
     1251            zn1 = ((BET031*zt   & 
     1252               &   + BET121*zs+BET021)*zt   & 
     1253               &   + (BET211*zs+BET111)*zs+BET011)*zt   & 
     1254               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001 
     1255               ! 
     1256            zn0 = ((((BET050*zt   & 
     1257               &   + BET140*zs+BET040)*zt   & 
     1258               &   + (BET230*zs+BET130)*zs+BET030)*zt   & 
     1259               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   & 
     1260               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   & 
     1261               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000 
     1262               ! 
     1263            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     1264            ! 
    8641265            pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
    8651266            ! 
     
    9971398      !! 
    9981399      INTEGER                              , INTENT(in   ) ::  Kmm   ! time level index 
    999       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
     1400      REAL(dp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
    10001401      REAL(wp), DIMENSION(:,:,:,:)         , INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1] 
    10011402      REAL(wp), DIMENSION(:,:,:)           , INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
     
    10211422      INTEGER                                , INTENT(in   ) ::  Kmm   ! time level index 
    10221423      INTEGER                                , INTENT(in   ) ::  ktab, ktn2 
    1023       REAL(wp), DIMENSION(jpi,jpj,  jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
     1424      REAL(dp), DIMENSION(jpi,jpj,  jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu] 
    10241425      REAL(wp), DIMENSION(A2D_T(ktab),JPK,JPTS), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1] 
    10251426      REAL(wp), DIMENSION(A2D_T(ktn2),JPK     ), INTENT(  out) ::  pn2   ! Brunt-Vaisala frequency squared [1/s^2] 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90

    r14986 r15066  
    101101         ! 
    102102         CALL eos_rab( CASTWP(pts(:,:,:,:,Kaa)), zab, Kmm )         ! after alpha and beta (given on T-points) 
    103          CALL bn2    ( CASTWP(pts(:,:,:,:,Kaa)), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points) 
     103         CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points) 
    104104         ! 
    105105         IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdfosm.F90

    r14986 r15066  
    32953295      ! w-level of the mixing and mixed layers 
    32963296      CALL eos_rab( CASTWP(ts(:,:,:,:,Kmm)), rab_n, Kmm ) 
    3297       CALL bn2( CASTWP(ts(:,:,:,:,Kmm)), rab_n, rn2, Kmm ) 
     3297      CALL bn2( ts(:,:,:,:,Kmm), rab_n, rn2, Kmm ) 
    32983298      imld_rst(:,:) = nlb10            ! Initialization to the number of w ocean point 
    32993299      hbl(:,:) = 0.0_wp                ! Here hbl used as a dummy variable, integrating vertically N^2 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/step.F90

    r14986 r15066  
    171171                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
    172172                         CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
    173                          CALL bn2    ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
    174                          CALL bn2    ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
     173                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     174                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
    175175 
    176176      !  VERTICAL PHYSICS 
  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/stpmlf.F90

    r14986 r15066  
    176176      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    177177      !  THERMODYNAMICS 
    178                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nbb)), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
    179                          CALL eos_rab( CASTWP(ts(:,:,:,:,Nnn)), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
    180                          CALL bn2    ( CASTWP(ts(:,:,:,:,Nbb)), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
    181                          CALL bn2    ( CASTWP(ts(:,:,:,:,Nnn)), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
     178                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points 
     179                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points 
     180                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency 
     181                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency 
    182182 
    183183      !  VERTICAL PHYSICS 
     
    217217                         CALL ssh_nxt    ( kstp, Nbb, Nnn, ssh,  Naa )   ! after ssh (includes call to div_hor) 
    218218      IF( .NOT.lk_linssh ) THEN 
    219                          CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts 
     219                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa)           )   ! "after" ssh/h_0 ratio at t,u,v pts 
    220220         IF( ln_dynspg_exp )   & 
    221             &            CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point 
     221            &            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) )   ! spg_exp : needed only for "now" ssh/h_0 ratio at f point 
    222222      ENDIF 
    223223                         CALL wzv        ( kstp, Nbb, Nnn, Naa, ww  )    ! Nnn cross-level velocity 
     
    227227                            zgdept(:,:,jk) = gdept(:,:,jk,Nnn) 
    228228                         END DO 
    229                          CALL eos        ( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, zgdept ) ! now in situ density for hpg computation 
     229                         CALL eos        ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 
    230230                         DEALLOCATE( zgdept ) 
    231231 
     
    268268                                      ! as well as vertical scale factors and vertical velocity need to be updated 
    269269                            CALL div_hor    ( kstp, Nbb, Nnn )                  ! Horizontal divergence  (2nd call in time-split case) 
    270             IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Naa)), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts 
     270            IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! update ssh/h_0 ratio at t,u,v,f pts 
    271271         ENDIF 
    272272                            CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion 
     
    303303      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    304304                         CALL ssh_atf    ( kstp, Nbb, Nnn, Naa, ssh )            ! time filtering of "now" sea surface height 
    305       IF(.NOT.lk_linssh) CALL dom_qco_r3c( CASTWP(ssh(:,:,Nnn)), r3t_f, r3u_f, r3v_f )   ! "now" ssh/h_0 ratio from filtrered ssh 
     305      IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f )   ! "now" ssh/h_0 ratio from filtrered ssh 
    306306#if defined key_top 
    307307      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
Note: See TracChangeset for help on using the changeset viewer.