Changeset 12489 for NEMO/trunk/src/OCE/DOM

Timestamp:

2020-02-28T16:55:11+01:00 (4 years 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/DOM

Files:

• daymod.F90 (modified) (7 diffs)
• dom_oce.F90 (modified) (2 diffs)
• domain.F90 (modified) (4 diffs)
• domvvl.F90 (modified) (12 diffs)
• istate.F90 (modified) (1 diff)
• phycst.F90 (modified) (1 diff)

NEMO/trunk/src/OCE/DOM/daymod.F90

@@ -20 +20 @@
    !!                    -------------------------------
    !!   sbcmod assume that the time step is dividing the number of second of 
-   !!   in a day, i.e. ===> MOD( rday, rdt ) == 0 
+   !!   in a day, i.e. ===> MOD( rday, rn_Dt ) == 0 
    !!   except when user defined forcing is used (see sbcmod.F90)
    !!----------------------------------------------------------------------
@@ -73 +73 @@
       !
       ! max number of seconds between each restart
-      IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN
+      IF( REAL( nitend - nit000 + 1 ) * rn_Dt > REAL( HUGE( nsec1jan000 ) ) ) THEN
          CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ',   &
             &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' )
@@ -79 +79 @@
       nsecd   = NINT(       rday )
       nsecd05 = NINT( 0.5 * rday )
-      ndt     = NINT(       rdt  )
-      ndt05   = NINT( 0.5 * rdt  )
+      ndt     = NINT(       rn_Dt  )
+      ndt05   = NINT( 0.5 * rn_Dt  )
 
       IF( .NOT. l_offline )   CALL day_rst( nit000, 'READ' )
@@ -239 +239 @@
       nsec_monday  = nsec_monday  + ndt
       nsec_day   = nsec_day   + ndt
-      adatrj  = adatrj  + rdt / rday
-      fjulday = fjulday + rdt / rday
+      adatrj  = adatrj  + rn_Dt / rday
+      fjulday = fjulday + rn_Dt / rday
       IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error
       IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error
@@ -309 +309 @@
       !!       In both those options, the  exact duration of the experiment
       !!       since the beginning (cumulated duration of all previous restart runs)
-      !!       is not stored in the restart and is assumed to be (nit000-1)*rdt.
+      !!       is not stored in the restart and is assumed to be (nit000-1)*rn_Dt.
       !!       This is valid is the time step has remained constant.
       !!
@@ -379 +379 @@
                isecond = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss)
                IF( isecond - ndt05 .lt. 0 )   ndastp = ndastp - 1      ! Start hour is specified in the namelist (default 0)
-               adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
+               adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday
                ! note this is wrong if time step has changed during run
             ENDIF
@@ -389 +389 @@
             isecond = ( nhour * NINT(rhhmm) + nminute ) * NINT(rmmss)
             IF( isecond - ndt05 .LT. 0 )   ndastp = ndastp - 1         ! Start hour is specified in the namelist (default 0)
-            adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday
+            adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday
          ENDIF
          IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error

NEMO/trunk/src/OCE/DOM/dom_oce.F90

@@ -33 +33 @@
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
    LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc)
-   REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer
+   REAL(wp), PUBLIC ::   rn_Dt          !: time step for the dynamics and tracer
    REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter
-   INTEGER , PUBLIC ::   nn_euler       !: =0 start with forward time step or not (=1)
+   LOGICAL , PUBLIC ::   ln_1st_euler   !: =T start with forward time step or not (=F)
    LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers
 
@@ -49 +49 @@
    LOGICAL,  PUBLIC :: ln_bt_auto       !: Set number of barotropic iterations automatically
    INTEGER,  PUBLIC :: nn_bt_flt        !: Filter choice
-   INTEGER,  PUBLIC :: nn_baro          !: Number of barotropic iterations during one baroclinic step (rdt)
+   INTEGER,  PUBLIC :: nn_e          !: Number of barotropic iterations during one baroclinic step (rn_Dt)
    REAL(wp), PUBLIC :: rn_bt_cmax       !: Maximum allowed courant number (used if ln_bt_auto=T)
    REAL(wp), PUBLIC :: rn_bt_alpha      !: Time stepping diffusion parameter
 
 
-   !                                   !! old non-DOCTOR names still used in the model
-   REAL(wp), PUBLIC ::   atfp           !: asselin time filter parameter
-   REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer
-
    !                                   !!! associated variables
-   INTEGER , PUBLIC ::   neuler         !: restart euler forward option (0=Euler)
-   REAL(wp), PUBLIC ::   r2dt           !: = 2*rdt except at nit000 (=rdt) if neuler=0
+   LOGICAL , PUBLIC ::   l_1st_euler    !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T)
+   REAL(wp), PUBLIC ::   rDt, r1_Dt     !: Current model timestep and reciprocal
+                                        !: rDt = 2 * rn_Dt if leapfrog and l_1st_euler = F
+                                        !:     =     rn_Dt if leapfrog and l_1st_euler = T
+                                        !:     =     rn_Dt if RK3
 
    !!----------------------------------------------------------------------

NEMO/trunk/src/OCE/DOM/domain.F90

@@ -287 +287 @@
          &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
-         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
+         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler  , &
          &             ln_cfmeta, ln_xios_read, nn_wxios
-      NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask
+      NAMELIST/namdom/ ln_linssh, rn_Dt, rn_atfp, ln_crs, ln_meshmask
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -317 +317 @@
          WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir )
          WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart
-         WRITE(numout,*) '      start with forward time step    nn_euler        = ', nn_euler
+         WRITE(numout,*) '      start with forward time step    ln_1st_euler    = ', ln_1st_euler
          WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl
          WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000
@@ -353 +353 @@
       nleapy = nn_leapy
       ninist = nn_istate
-      neuler = nn_euler
-      IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN
+      l_1st_euler = ln_1st_euler
+      IF( .NOT. l_1st_euler .AND. .NOT. ln_rstart ) THEN
          IF(lwp) WRITE(numout,*)  
          IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
-         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : nn_euler is forced to 0 '   
-         neuler = 0
+         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : l_1st_euler is forced to .true. '   
+         l_1st_euler = .true.
       ENDIF
       !                             ! control of output frequency
@@ -408 +408 @@
          WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
          WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt
+         WRITE(numout,*) '      ocean time step                         rn_Dt       = ', rn_Dt
          WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
          WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
       ENDIF
       !
-      !          ! conversion DOCTOR names into model names (this should disappear soon)
-      atfp = rn_atfp
-      rdt  = rn_rdt
+      !! Initialise current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
+      rDt  = 2._wp * rn_Dt
+      r1_Dt = 1._wp / rDt
 
       IF( TRIM(Agrif_CFixed()) == '0' ) THEN

NEMO/trunk/src/OCE/DOM/domvvl.F90

@@ -235 +235 @@
          IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
             frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
-            frq_rst_hdv(:,:) = 1._wp / rdt
+            frq_rst_hdv(:,:) = 1._wp / rn_Dt
          ENDIF
          IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
@@ -247 +247 @@
                   ! values inside the equatorial band (ztilde as zstar)
                   frq_rst_e3t(ji,jj) =  0.0_wp
-                  frq_rst_hdv(ji,jj) =  1.0_wp / rdt
+                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt
                ELSE                                      ! transition band (2.5 to 6 degrees N/S)
                   !                                      ! (linearly transition from z-tilde to z-star)
@@ -253 +253 @@
                      &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
                      &                                          * 180._wp / 3.5_wp ) )
-                  frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
-                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
+                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                &
+                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   &
                      &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
                      &                                          * 180._wp / 3.5_wp ) )
@@ -264 +264 @@
                   ij0 = 128   ;   ij1 = 135   ;   
                   frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
-                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
+                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rn_Dt
                ENDIF
             ENDIF
@@ -319 +319 @@
       INTEGER                ::   ji, jj, jk            ! dummy loop indices
       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
-      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
+      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
       LOGICAL                ::   ll_do_bclinic         ! local logical
       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
@@ -373 +373 @@
             IF( kt > nit000 ) THEN
                DO jk = 1, jpkm1
-                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
+                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
                      &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
                END DO
@@ -430 +430 @@
          ! 4 - Time stepping of baroclinic scale factors
          ! ---------------------------------------------
-         ! Leapfrog time stepping
-         ! ~~~~~~~~~~~~~~~~~~~~~~
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            z2dt =  rdt
-         ELSE
-            z2dt = 2.0_wp * rdt
-         ENDIF
          CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
+         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
 
          ! Maximum deformation control
@@ -624 +617 @@
       ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
+         IF( l_1st_euler ) THEN
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
          ELSE
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
@@ -821 +814 @@
                   e3t(:,:,:,Kbb) = e3t_0(:,:,:)
                END WHERE
-               IF( neuler == 0 ) THEN
+               IF( l_1st_euler ) THEN
                   e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                ENDIF
@@ -827 +820 @@
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
                e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
-               neuler = 0
+               l_1st_euler = .true.
             ELSE IF( id2 > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-               neuler = 0
+               l_1st_euler = .true.
             ELSE
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                DO jk = 1, jpk
                   e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
@@ -848 +841 @@
                END DO
                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
-               neuler = 0
+               l_1st_euler = .true.
             ENDIF
             !                             ! ----------- !
@@ -1008 +1001 @@
             WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
             WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
-            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
+            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
          ELSE
             WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t

NEMO/trunk/src/OCE/DOM/istate.F90

@@ -92 +92 @@
          !                                    ! ---------------
          numror = 0                           ! define numror = 0 -> no restart file to read
-         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
+         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !                                    ! Initialization of ocean to zero

NEMO/trunk/src/OCE/DOM/phycst.F90

@@ -39 +39 @@
    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin]
 
-   REAL(wp), PUBLIC ::   rau0                        !: volumic mass of reference     [kg/m3]
-   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg]
+   REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference     [kg/m3]
+   REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                   [m3/kg]
    REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin]
    REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J]
-   REAL(wp), PUBLIC ::   rau0_rcp                    !: = rau0 * rcp 
-   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp )
+   REAL(wp), PUBLIC ::   rho0_rcp                    !: = rho0 * rcp 
+   REAL(wp), PUBLIC ::   r1_rho0_rcp                 !: = 1. / ( rho0 * rcp )
 
    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice (not used?)