Ignore:
Timestamp:
2020-02-28T16:55:11+01:00 (12 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.

Location:
NEMO/trunk/src/OCE/TRA
Files:
12 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 
  • NEMO/trunk/src/OCE/TRA/traadv.F90

    r12377 r12489  
    9292      IF( ln_timing )   CALL timing_start('tra_adv') 
    9393      ! 
    94       !                                          ! set time step 
    95       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =         rdt   ! at nit000             (Euler) 
    96       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp * rdt   ! at nit000 or nit000+1 (Leapfrog) 
    97       ENDIF 
    98       ! 
    9994      !                                         !==  effective transport  ==! 
    10095      zuu(:,:,jpk) = 0._wp 
     
    149144         CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 
    150145      CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order 
    151          CALL tra_adv_fct    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 
     146         CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 
    152147      CASE ( np_MUS )                                 ! MUSCL 
    153          CALL tra_adv_mus    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )  
     148         CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )  
    154149      CASE ( np_UBS )                                 ! UBS 
    155          CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   ) 
     150         CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   ) 
    156151      CASE ( np_QCK )                                 ! QUICKEST 
    157          CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 
     152         CALL tra_adv_qck    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 
    158153      ! 
    159154      END SELECT 
  • NEMO/trunk/src/OCE/TRA/traadv_fct.F90

    r12377 r12489  
    2020   USE diaptr         ! poleward transport diagnostics 
    2121   USE diaar5         ! AR5 diagnostics 
    22    USE phycst  , ONLY : rau0_rcp 
     22   USE phycst  , ONLY : rho0_rcp 
    2323   USE zdf_oce , ONLY : ln_zad_Aimp 
    2424   ! 
  • NEMO/trunk/src/OCE/TRA/traatf.F90

    r12377 r12489  
    113113      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries 
    114114  
    115       ! set time step size (Euler/Leapfrog) 
    116       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler) 
    117       ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog) 
    118       ENDIF 
    119  
    120115      ! trends computation initialisation 
    121116      IF( l_trdtra )   THEN                     
     
    128123         ENDIF 
    129124         ! total trend for the non-time-filtered variables.  
    130          zfact = 1.0 / rdt 
     125         zfact = 1.0 / rn_Dt 
    131126         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms 
    132127         DO jk = 1, jpkm1 
     
    144139      ENDIF 
    145140 
    146       IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping  
     141      IF( l_1st_euler ) THEN       ! Euler time-stepping  
    147142         ! 
    148143         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl 
     
    156151      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    157152         ! 
    158          IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface  
    159          ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
     153         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,        'TRA', pts, jpts )  ! linear free surface  
     154         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rn_Dt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface 
    160155         ENDIF 
    161156         ! 
     
    167162      ! 
    168163      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt      
    169          zfact = 1._wp / r2dt              
     164         zfact = 1._wp / rDt              
    170165         DO jk = 1, jpkm1 
    171166            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kmm) - ztrdt(:,:,jk) ) * zfact 
     
    219214            ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers 
    220215            ! 
    221             pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! pt <-- filtered pt 
     216            pt(ji,jj,jk,jn,Kmm) = ztn + rn_atfp * ztd                      ! pt <-- filtered pt 
    222217         END_3D 
    223218         ! 
     
    234229      !!  
    235230      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields. 
    236       !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) 
    237       !!                       /( e3t(Kmm)         + atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] ) 
     231      !!             pt(Kmm)  = ( e3t(Kmm)*pt(Kmm) + rn_atfp*[ e3t(Kbb)*pt(Kbb) - 2 e3t(Kmm)*pt(Kmm) + e3t_a*pt(Kaa) ] ) 
     232      !!                       /( e3t(Kmm)         + rn_atfp*[ e3t(Kbb)         - 2 e3t(Kmm)         + e3t(Kaa)    ] ) 
    238233      !! 
    239234      !! ** Action  : - pt(Kmm) ready for the next time step 
     
    277272      ENDIF 
    278273      zfact = 1._wp / p2dt 
    279       zfact1 = atfp * p2dt 
    280       zfact2 = zfact1 * r1_rau0 
     274      zfact1 = rn_atfp * p2dt 
     275      zfact2 = zfact1 * r1_rho0 
    281276      DO jn = 1, kjpt       
    282277         DO_3D_00_00( 1, jpkm1 ) 
     
    292287            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b 
    293288            ! 
    294             ze3t_f = ze3t_n + atfp * ze3t_d 
    295             ztc_f  = ztc_n  + atfp * ztc_d 
     289            ze3t_f = ze3t_n + rn_atfp * ze3t_d 
     290            ztc_f  = ztc_n  + rn_atfp * ztc_d 
    296291            ! 
    297292            ! Add asselin correction on scale factors: 
  • NEMO/trunk/src/OCE/TRA/trabbc.F90

    r12377 r12489  
    6666      !!       ocean bottom can be computed once and is added to the temperature 
    6767      !!       trend juste above the bottom at each time step: 
    68       !!            ta = ta + Qsf / (rau0 rcp e3T) for k= mbkt 
     68      !!            ta = ta + Qsf / (rho0 rcp e3T) for k= mbkt 
    6969      !!       Where Qsf is the geothermal heat flux. 
    7070      !! 
     
    102102      ENDIF 
    103103      ! 
    104       CALL iom_put ( "hfgeou" , rau0_rcp * qgh_trd0(:,:) ) 
     104      CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 
    105105      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
    106106      ! 
     
    162162         CASE ( 1 )                          !* constant flux 
    163163            IF(lwp) WRITE(numout,*) '   ==>>>   constant heat flux  =   ', rn_geoflx_cst 
    164             qgh_trd0(:,:) = r1_rau0_rcp * rn_geoflx_cst 
     164            qgh_trd0(:,:) = r1_rho0_rcp * rn_geoflx_cst 
    165165            ! 
    166166         CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 
     
    179179 
    180180            CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data 
    181             qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 
     181            qgh_trd0(:,:) = r1_rho0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 
    182182            ! 
    183183         CASE DEFAULT 
  • NEMO/trunk/src/OCE/TRA/traldf_iso.F90

    r12377 r12489  
    109109      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
    110110      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
    111       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      - 
     111      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      - 
    112112      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d 
    113113      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw  
     
    129129         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    130130      ! 
    131       !                                            ! set time step size (Euler/Leapfrog) 
    132       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    133       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
    134       ENDIF 
    135       z1_2dt = 1._wp / z2dt 
    136131      ! 
    137132      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     
    178173               DO_3D_10_10( 2, jpkm1 ) 
    179174                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    180                   zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    181                   akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     175                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     176                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 
    182177               END_3D 
    183178           ENDIF 
  • NEMO/trunk/src/OCE/TRA/traldf_triad.F90

    r12377 r12489  
    8686      INTEGER  ::  ip,jp,kp         ! dummy loop indices 
    8787      INTEGER  ::  ierr            ! local integer 
    88       REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars 
    89       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      - 
    90       REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      - 
     88      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars 
     89      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      - 
     90      REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      - 
    9191      ! 
    9292      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv 
     
    111111      l_hst = .FALSE. 
    112112      l_ptr = .FALSE. 
    113       IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )      l_ptr = .TRUE.  
    114       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    115          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE. 
    116       ! 
    117       !                                                        ! set time step size (Euler/Leapfrog) 
    118       IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler) 
    119       ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog) 
     113      IF( cdtype == 'TRA' ) THEN 
     114         IF( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf') )      l_ptr = .TRUE.  
     115         IF( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.                   & 
     116         &   iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  )   l_hst = .TRUE. 
    120117      ENDIF 
    121       z1_2dt = 1._wp / z2dt 
    122118      ! 
    123119      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
     
    189185               DO_3D_10_10( 2, jpkm1 ) 
    190186                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    191                   zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    192                   akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 
     187                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     188                  akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_Dt 
    193189               END_3D 
    194190           ENDIF 
  • NEMO/trunk/src/OCE/TRA/tramle.F90

    r12377 r12489  
    4141 
    4242   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation 
    43    REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld 
     43   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld 
    4444   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_mle=1 case 
    4545 
     
    112112         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
    113113         zmld(ji,jj) = zmld(ji,jj) + zc 
    114          zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0 
     114         zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0 
    115115         zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp 
    116116      END_3D 
     
    273273      IF( ln_mle ) THEN                ! MLE initialisation 
    274274         ! 
    275          rb_c = grav * rn_rho_c_mle /rau0        ! Mixed Layer buoyancy criteria 
     275         rb_c = grav * rn_rho_c_mle /rho0        ! Mixed Layer buoyancy criteria 
    276276         IF(lwp) WRITE(numout,*) 
    277277         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 ' 
  • NEMO/trunk/src/OCE/TRA/tranpc.F90

    r12377 r12489  
    6767      LOGICAL  ::   l_bottom_reached, l_column_treated 
    6868      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
    69       REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     69      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_rDt 
    7070      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0) 
    7171      REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point... 
     
    301301         ! 
    302302         IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic 
    303             z1_r2dt = 1._wp / (2._wp * rdt) 
    304             ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_r2dt 
    305             ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_r2dt 
     303            z1_rDt = 1._wp / (2._wp * rn_Dt) 
     304            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt 
     305            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt 
    306306            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt ) 
    307307            CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds ) 
  • NEMO/trunk/src/OCE/TRA/traqsr.F90

    r12377 r12489  
    8787      !!         I(k) = Qsr*( rn_abs*EXP(z(k)/rn_si0) + (1.-rn_abs)*EXP(z(k)/rn_si1) ) 
    8888      !!         The temperature trend associated with the solar radiation penetration  
    89       !!         is given by : zta = 1/e3t dk[ I ] / (rau0*Cp) 
     89      !!         is given by : zta = 1/e3t dk[ I ] / (rho0*Cp) 
    9090      !!         At the bottom, boudary condition for the radiation is no flux : 
    9191      !!      all heat which has not been absorbed in the above levels is put 
     
    135135      !                         !-----------------------------------! 
    136136      IF( kt == nit000 ) THEN          !==  1st time step  ==! 
    137 !!gm case neuler  not taken into account.... 
    138          IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart 
     137         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart 
    139138            IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file' 
    140139            z1_2 = 0.5_wp 
     
    156155         ! 
    157156         DO jk = 1, nksr 
    158             qsr_hc(:,:,jk) = r1_rau0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 
     157            qsr_hc(:,:,jk) = r1_rho0_rcp * ( etot3(:,:,jk) - etot3(:,:,jk+1) ) 
    159158         END DO 
    160159         ! 
     
    228227         ! 
    229228         DO_3D_00_00( 1, nksr ) 
    230             qsr_hc(ji,jj,jk) = r1_rau0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
     229            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( zea(ji,jj,jk) - zea(ji,jj,jk+1) ) 
    231230         END_3D 
    232231         ! 
     
    235234      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
    236235         ! 
    237          zz0 =        rn_abs   * r1_rau0_rcp      ! surface equi-partition in 2-bands 
    238          zz1 = ( 1. - rn_abs ) * r1_rau0_rcp 
     236         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands 
     237         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 
    239238         DO_3D_00_00( 1, nksr ) 
    240239            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r ) 
     
    253252      ! sea-ice: store the 1st ocean level attenuation coefficient 
    254253      DO_2D_00_00 
    255          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) 
     254         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 
    256255         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
    257256         ENDIF 
     
    263262         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero 
    264263         DO jk = nksr, 1, -1 
    265             zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rau0_rcp 
     264            zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 
    266265         END DO          
    267266         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation 
  • NEMO/trunk/src/OCE/TRA/trasbc.F90

    r12377 r12489  
    124124      !                             !==  Now sbc tracer content fields  ==! 
    125125      DO_2D_01_00 
    126          sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)   ! non solar heat flux 
    127          sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
     126         sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux 
     127         sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
    128128      END_2D 
    129129      IF( ln_linssh ) THEN                !* linear free surface   
    130130         DO_2D_01_00 
    131             sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 
    132             sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 
     131            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 
     132            sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 
    133133         END_2D 
    134134         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 
  • NEMO/trunk/src/OCE/TRA/trazdf.F90

    r12377 r12489  
    6666      ENDIF 
    6767      ! 
    68       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping) 
    69       ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog) 
    70       ENDIF 
    71       ! 
    7268      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
    7369         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
     
    7773      ! 
    7874      !                                      !* compute lateral mixing trend and add it to the general trend 
    79       CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )  
     75      CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )  
    8076 
    8177!!gm WHY here !   and I don't like that ! 
     
    8985         DO jk = 1, jpkm1 
    9086            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 
    91                &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 
     87               &          / (e3t(:,:,jk,Kmm)*rDt) ) - ztrdt(:,:,jk) 
    9288            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 
    93               &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 
     89              &           / (e3t(:,:,jk,Kmm)*rDt) ) - ztrds(:,:,jk) 
    9490         END DO 
    9591!!gm this should be moved in trdtra.F90 and done on all trends 
Note: See TracChangeset for help on using the changeset viewer.