Ignore:
Timestamp:
2020-02-28T16:55:11+01:00 (11 months ago)
Author:
davestorkey
Message:

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/TRA/eosbn2.F90

    r12377 r12489  
    191191      !!                   ***  ROUTINE eos_insitu  *** 
    192192      !! 
    193       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
     193      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
    194194      !!       potential temperature and salinity using an equation of state 
    195195      !!       selected in the nameos namelist 
    196196      !! 
    197       !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0 
     197      !! ** Method  :   prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 
    198198      !!         with   prd    in situ density anomaly      no units 
    199199      !!                t      TEOS10: CT or EOS80: PT      Celsius 
     
    201201      !!                z      depth                        meters 
    202202      !!                rho    in situ density              kg/m^3 
    203       !!                rau0   reference density            kg/m^3 
     203      !!                rho0   reference density            kg/m^3 
    204204      !! 
    205205      !!     ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 
     
    210210      !! 
    211211      !!     ln_seos : simplified equation of state 
    212       !!              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 
    213213      !!              linear case function of T only: rn_alpha<>0, other coefficients = 0 
    214214      !!              linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 
     
    267267            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    268268            ! 
    269             prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm  ! density anomaly (masked) 
     269            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked) 
    270270            ! 
    271271         END_3D 
     
    283283               &  - rn_nu * zt * zs 
    284284               !                                  
    285             prd(ji,jj,jk) = zn * r1_rau0 * ztm                ! density anomaly (masked) 
     285            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked) 
    286286         END_3D 
    287287         ! 
     
    299299      !!                  ***  ROUTINE eos_insitu_pot  *** 
    300300      !! 
    301       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
     301      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the 
    302302      !!      potential volumic mass (Kg/m3) from potential temperature and 
    303303      !!      salinity fields using an equation of state selected in the 
     
    379379                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
    380380                  ! 
    381                   prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
     381                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked) 
    382382               END DO 
    383383               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     
    419419               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    420420               ! 
    421                prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     421               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked) 
    422422            END_3D 
    423423         ENDIF 
     
    434434               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   & 
    435435               &  - rn_nu * zt * zs 
    436             prhop(ji,jj,jk) = ( rau0 + zn ) * ztm 
     436            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm 
    437437            !                                                     ! density anomaly (masked) 
    438438            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh 
    439             prd(ji,jj,jk) = zn * r1_rau0 * ztm 
     439            prd(ji,jj,jk) = zn * r1_rho0 * ztm 
    440440            ! 
    441441         END_3D 
     
    454454      !!                  ***  ROUTINE eos_insitu_2d  *** 
    455455      !! 
    456       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
     456      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from 
    457457      !!      potential temperature and salinity using an equation of state 
    458458      !!      selected in the nameos namelist. * 2D field case 
     
    508508            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    509509            ! 
    510             prd(ji,jj) = zn * r1_rau0 - 1._wp               ! unmasked in situ density anomaly 
     510            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly 
    511511            ! 
    512512         END_2D 
     
    524524               &  - rn_nu * zt * zs 
    525525               ! 
    526             prd(ji,jj) = zn * r1_rau0               ! unmasked in situ density anomaly 
     526            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly 
    527527            ! 
    528528         END_2D 
     
    588588            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    589589            ! 
    590             pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm 
     590            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm 
    591591            ! 
    592592            ! beta 
     
    609609            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    610610            ! 
    611             pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm 
     611            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm 
    612612            ! 
    613613         END_3D 
     
    622622            ! 
    623623            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    624             pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm   ! alpha 
     624            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha 
    625625            ! 
    626626            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    627             pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm   ! beta 
     627            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta 
    628628            ! 
    629629         END_3D 
     
    694694            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    695695            ! 
    696             pab(ji,jj,jp_tem) = zn * r1_rau0 
     696            pab(ji,jj,jp_tem) = zn * r1_rho0 
    697697            ! 
    698698            ! beta 
     
    715715            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    716716            ! 
    717             pab(ji,jj,jp_sal) = zn / zs * r1_rau0 
     717            pab(ji,jj,jp_sal) = zn / zs * r1_rho0 
    718718            ! 
    719719            ! 
     
    729729            ! 
    730730            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    731             pab(ji,jj,jp_tem) = zn * r1_rau0   ! alpha 
     731            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha 
    732732            ! 
    733733            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    734             pab(ji,jj,jp_sal) = zn * r1_rau0   ! beta 
     734            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta 
    735735            ! 
    736736         END_2D 
     
    799799         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    800800         ! 
    801          pab(jp_tem) = zn * r1_rau0 
     801         pab(jp_tem) = zn * r1_rho0 
    802802         ! 
    803803         ! beta 
     
    820820         zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    821821         ! 
    822          pab(jp_sal) = zn / zs * r1_rau0 
     822         pab(jp_sal) = zn / zs * r1_rho0 
    823823         ! 
    824824         ! 
     
    831831         ! 
    832832         zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs 
    833          pab(jp_tem) = zn * r1_rau0   ! alpha 
     833         pab(jp_tem) = zn * r1_rho0   ! alpha 
    834834         ! 
    835835         zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt 
    836          pab(jp_sal) = zn * r1_rau0   ! beta 
     836         pab(jp_sal) = zn * r1_rho0   ! beta 
    837837         ! 
    838838      CASE DEFAULT 
     
    10521052      !! ** Method  :   PE is defined analytically as the vertical  
    10531053      !!                   primitive of EOS times -g integrated between 0 and z>0. 
    1054       !!                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 
    10551055      !!                                                      = 1/z * /int_0^z rd dz - rd  
    10561056      !!                                where rd is the density anomaly (see eos_rhd function) 
    10571057      !!                ab_pe are partial derivatives of PE anomaly with respect to T and S: 
    1058       !!                    ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT 
    1059       !!                    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 
    10601060      !! 
    10611061      !! ** Action  : - pen         : PE anomaly given at T-points 
     
    11031103            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    11041104            ! 
    1105             ppen(ji,jj,jk)  = zn * zh * r1_rau0 * ztm 
     1105            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm 
    11061106            ! 
    11071107            ! alphaPE non-linear anomaly 
     
    11181118            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    11191119            !                               
    1120             pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm 
     1120            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm 
    11211121            ! 
    11221122            ! betaPE non-linear anomaly 
     
    11331133            zn  = ( zn2 * zh + zn1 ) * zh + zn0 
    11341134            !                               
    1135             pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm 
     1135            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm 
    11361136            ! 
    11371137         END_3D 
     
    11441144            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point 
    11451145            ztm = tmask(ji,jj,jk)                ! tmask 
    1146             zn  = 0.5_wp * zh * r1_rau0 * ztm 
     1146            zn  = 0.5_wp * zh * r1_rho0 * ztm 
    11471147            !                                    ! Potential Energy 
    11481148            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn 
     
    11861186      IF(lwm) WRITE( numond, nameos ) 
    11871187      ! 
    1188       rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     1188      rho0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
    11891189      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
    11901190      ! 
     
    15981598            WRITE(numout,*) '   ==>>>   use of simplified eos:    ' 
    15991599            WRITE(numout,*) '              rhd(dT=T-10,dS=S-35,Z) = [-a0*(1+lambda1/2*dT+mu1*Z)*dT ' 
    1600             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' 
    16011601            WRITE(numout,*) '              with the following coefficients :' 
    16021602            WRITE(numout,*) '                 thermal exp. coef.    rn_a0      = ', rn_a0 
     
    16171617      END SELECT 
    16181618      ! 
    1619       rau0_rcp    = rau0 * rcp  
    1620       r1_rau0     = 1._wp / rau0 
     1619      rho0_rcp    = rho0 * rcp  
     1620      r1_rho0     = 1._wp / rho0 
    16211621      r1_rcp      = 1._wp / rcp 
    1622       r1_rau0_rcp = 1._wp / rau0_rcp  
     1622      r1_rho0_rcp = 1._wp / rho0_rcp  
    16231623      ! 
    16241624      IF(lwp) THEN 
     
    16351635      IF(lwp) WRITE(numout,*) 
    16361636      IF(lwp) WRITE(numout,*) '   Associated physical constant' 
    1637       IF(lwp) WRITE(numout,*) '      volumic mass of reference           rau0  = ', rau0   , ' kg/m^3' 
    1638       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' 
    16391639      IF(lwp) WRITE(numout,*) '      ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
    1640       IF(lwp) WRITE(numout,*) '      rau0 * rcp                       rau0_rcp = ', rau0_rcp 
    1641       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 
    16421642      ! 
    16431643   END SUBROUTINE eos_init 
Note: See TracChangeset for help on using the changeset viewer.