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.

Location:
NEMO/trunk/src/OCE/ZDF
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/ZDF/zdfdrg.F90

    r12377 r12489  
    165165      !!--------------------------------------------------------------------- 
    166166      ! 
    167 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 
    168       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     167!!gm bug : time step is only rn_Dt (not 2 rn_Dt if euler start !) 
     168      zm1_2dt = - 1._wp / ( 2._wp * rn_Dt ) 
    169169 
    170170      IF( l_trddyn ) THEN      ! trends: store the input trends 
  • NEMO/trunk/src/OCE/ZDF/zdfgls.F90

    r12377 r12489  
    170170         ! 
    171171         ! surface friction 
    172          ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 
     172         ustar2_surf(ji,jj) = r1_rho0 * taum(ji,jj) * tmask(ji,jj,1) 
    173173         !    
    174174!!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that... 
     
    267267         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    268268         !                                        ! diagonal 
    269          zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk)  
     269         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rn_Dt * zdiss * wmask(ji,jj,jk)  
    270270         !                                        ! right hand side in en 
    271          en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     271         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 
    272272      END_3D 
    273273      ! 
     
    477477         zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk,Kmm) ) 
    478478         !                                               ! diagonal 
    479          zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) 
     479         zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 
    480480         !                                               ! right hand side in psi 
    481          psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk) 
     481         psi(ji,jj,jk) = psi(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) 
    482482      END_3D 
    483483      ! 
     
    10071007      rc04  = rc03 * rc0 
    10081008      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1009       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
     1009      rsbc_tke2 = rn_Dt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
    10101010      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    10111011      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
    10121012      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    10131013      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
    1014       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 
    1015       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
    1016       ! 
    1017       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke 
    1018       rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke 
     1014      rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi 
     1015      rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking  
     1016      ! 
     1017      rfact_tke = -0.5_wp / rsc_tke * rn_Dt                                ! Cst used for the Diffusion term of tke 
     1018      rfact_psi = -0.5_wp / rsc_psi * rn_Dt                                ! Cst used for the Diffusion term of tke 
    10191019      ! 
    10201020      !                                !* Wall proximity function 
  • NEMO/trunk/src/OCE/ZDF/zdfiwm.F90

    r12377 r12489  
    8787      !!              This is divided into three components: 
    8888      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
    89       !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 
     89      !!                     zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm ) 
    9090      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 
    9191      !!              where hcri_iwm is the characteristic length scale of the bottom  
    9292      !!              intensification, ecri_iwm a map of available power, and H the ocean depth. 
    9393      !!                 2. Pycnocline-intensified low-mode dissipation 
    94       !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     94      !!                     zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
    9595      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
    9696      !!              where epyc_iwm is a map of available power, and nn_zpyc 
     
    9898      !!              energy dissipation. 
    9999      !!                 3. WKB-height dependent high mode dissipation 
    100       !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
     100      !!                     zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
    101101      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 
    102102      !!              where hbot_iwm is the characteristic length scale of the WKB bottom  
     
    151151      DO_2D_11_11 
    152152         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    153          zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     153         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
    154154         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 
    155155      END_2D 
     
    181181         ! 
    182182         DO_2D_11_11 
    183             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     183            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    184184         END_2D 
    185185         ! 
     
    196196         ! 
    197197         DO_2D_11_11 
    198             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     198            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    199199         END_2D 
    200200         ! 
     
    243243      ! 
    244244      DO_2D_11_11 
    245          IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     245         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    246246      END_2D 
    247247      ! 
     
    255255      ! Calculate molecular kinematic viscosity 
    256256      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  & 
    257          &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rau0 
     257         &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rho0 
    258258      DO jk = 2, jpkm1 
    259259         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
     
    293293         END_3D 
    294294         CALL mpp_sum( 'zdfiwm', zztmp ) 
    295          zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     295         zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    296296         ! 
    297297         IF(lwp) THEN 
     
    337337                                    !* output useful diagnostics: Kz*N^2 ,  
    338338!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
    339                                     !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm) 
     339                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    340340      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    341341         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
     
    345345            z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 
    346346         END DO 
    347          z2d(:,:) = rau0 * z2d(:,:) 
     347         z2d(:,:) = rho0 * z2d(:,:) 
    348348         CALL iom_put( "bflx_iwm", z3d ) 
    349349         CALL iom_put( "pcmap_iwm", z2d ) 
  • NEMO/trunk/src/OCE/ZDF/zdfmxl.F90

    r12377 r12489  
    9797      nmln(:,:)  = nlb10               ! Initialization to the number of w ocean point 
    9898      hmlp(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2 
    99       zN2_c = grav * rho_c * r1_rau0   ! convert density criteria into N^2 criteria 
     99      zN2_c = grav * rho_c * r1_rho0   ! convert density criteria into N^2 criteria 
    100100      DO_3D_11_11( nlb10, jpkm1 ) 
    101101         ikt = mbkt(ji,jj) 
  • NEMO/trunk/src/OCE/ZDF/zdfosm.F90

    r12377 r12489  
    300300     DO_2D_00_00 
    301301        ! Surface downward irradiance (so always +ve) 
    302         zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp 
     302        zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp 
    303303        ! Downwards irradiance at base of boundary layer 
    304304        zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) ) 
     
    312312        zbeta    = rab_n(ji,jj,1,jp_sal) 
    313313        ! Upwards surface Temperature flux for non-local term 
    314         zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 
     314        zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1) 
    315315        ! Upwards surface salinity flux for non-local term 
    316         zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 
     316        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1) 
    317317        ! Non radiative upwards surface buoyancy flux 
    318318        zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj) 
     
    324324        zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj) 
    325325        ! Surface upward velocity fluxes 
    326         zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
    327         zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * tmask(ji,jj,1) 
     326        zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
     327        zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1) 
    328328        ! Friction velocity (zustar), at T-point : LMD94 eq. 2 
    329329        zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 ) 
     
    441441                        &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) ) 
    442442 
    443                    zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     443                   zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    444444                        &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    445445! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment. 
     
    447447!                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) ) 
    448448 
    449 !                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
     449!                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / & 
    450450!                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    451451                   zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) ) 
     
    458458                   IF ( zzdhdt < 0._wp ) THEN 
    459459                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL 
    460                       zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
     460                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) 
    461461                   ELSE 
    462                       zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
     462                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) & 
    463463                           &  + MAX( zdb_bl(ji,jj), 0.0 ) 
    464464                   ENDIF 
     
    472472      ibld(:,:) = 3 
    473473 
    474       zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 
     474      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it 
    475475      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 
    476       zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
     476      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 
    477477 
    478478      DO_3D_00_00( 4, jpkm1 ) 
     
    496496            IF ( lconv(ji,jj) ) THEN 
    497497!unstable 
    498                zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) & 
     498               zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) & 
    499499                    &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 
    500500 
     
    503503                       & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 
    504504 
    505                   zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 
     505                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 
    506506                  zhbl_s = MIN(zhbl_s, ht(ji,jj)) 
    507507 
     
    12501250            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift 
    12511251            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift 
    1252             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
     1252            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    12531253         ! Stokes drift read in from sbcwave  (=2). 
    12541254         CASE(2) 
    12551255            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift 
    12561256            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift 
    1257             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* & 
     1257            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* & 
    12581258                 & SQRT(ut0sd**2 + vt0sd**2 ) ) 
    12591259         END SELECT 
     
    12711271         IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale 
    12721272         IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale 
    1273          IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
    1274          IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*tmask(:,:,1)*zustar**2*zustke ) 
     1273         IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine 
     1274         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke ) 
    12751275         IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine 
    12761276         IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine 
     
    15001500     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point 
    15011501     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
    1502      zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria 
     1502     zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria 
    15031503     ! 
    15041504     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2 
  • NEMO/trunk/src/OCE/ZDF/zdfric.F90

    r12377 r12489  
    174174         ! 
    175175         DO_2D_00_00 
    176             zustar = SQRT( taum(ji,jj) * r1_rau0 ) 
     176            zustar = SQRT( taum(ji,jj) * r1_rho0 ) 
    177177            zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth 
    178178            zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range 
  • NEMO/trunk/src/OCE/ZDF/zdftke.F90

    r12377 r12489  
    206206      !!-------------------------------------------------------------------- 
    207207      ! 
    208       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    209       zfact1 = -.5_wp * rdt  
    210       zfact2 = 1.5_wp * rdt * rn_ediss 
     208      zbbrau = rn_ebb / rho0       ! Local constant initialisation 
     209      zfact1 = -.5_wp * rn_Dt  
     210      zfact2 = 1.5_wp * rn_Dt * rn_ediss 
    211211      zfact3 = 0.5_wp       * rn_ediss 
    212212      ! 
     
    228228      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    229229      ! 
    230       !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
     230      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
    231231      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 
    232232      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 
     
    237237            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 
    238238            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 
    239             !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
     239            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 
    240240            zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  & 
    241241               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) 
     
    246246               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 
    247247               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 
    248                !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
     248               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0) 
    249249               zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  & 
    250250                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  ) 
     
    288288                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
    289289                  !                                           ! TKE Langmuir circulation source term 
    290                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     290                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    291291               ENDIF 
    292292            ENDIF 
     
    325325         ! 
    326326         !                                   ! right hand side in en 
    327          en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear 
     327         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear 
    328328            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification 
    329329            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation 
     
    439439      zmxld(:,:,:)  = rmxl_min 
    440440      ! 
    441       IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    442          zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     441      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 
     442         zraug = vkarmn * 2.e5_wp / ( rho0 * grav ) 
    443443         DO_2D_00_00 
    444444            zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
Note: See TracChangeset for help on using the changeset viewer.