Corrections for tipping bucket and nudging in very long simulations (#2063)

TYPE: bug fix

KEYWORDS: precipitation, tipping bucket, nudging, spectral nudging, analysis nudging, grid nudging, regional climate, dynamical downscaling, downscaling

SOURCE: Tanya Spero (U.S. EPA)

DESCRIPTION OF CHANGES:
Problem:
Several processes in WRF are currently triggered at periodic intervals
(such as reading/writing files and managing certain bookkeeping
processes). Some of those processes had been coded to identify those
time triggers by referencing a variable that contains the time elapsed
since the model simulation was initialized. That variable, XTIME, is a
Fortran single-precision real variable that counts the number of elapsed
minutes since initialization. However, single-precision real numbers
become imprecise (i.e., cannot accurately resolve "whole" numbers) after
they exceed 2^24, which is 16,777,216. In long simulations, that occurs
just before 32 years of simulation period.

Solution:

Two new variables are introduced based on existing variable CURR_SECS2: CURR_SECS2_R8, and double precision version of CURR_SECS2, to address bucket tipping and CURR_MINS2, based on CURR_SECS2_R8, to address issues in nudging code.

LIST OF MODIFIED FILES: 
M     phys/module_diag_misc.F
M     dyn_em/module_first_rk_step_part1.F 
M     dyn_em/solve_em.F,
M     phys/module_fdda_psufddagd.F
M     phys/module_fdda_spnudging.F,
M     phys/module_fddagd_driver.F
M     wrftladj/wrftladj/solve_em_ad.F

TESTS CONDUCTED: 
1. Modified code has been tested extensively using a long simulation (using a restart at more than 31 years into the simulation) with WRFv4.5.1 and with shorter (3-day) simulations using WRFv4.6. 
2. The Jenkins tests are all passing.

RELEASE NOTE: Corrected algorithms in the tipping bucket for precipitation and in the nudging routines to adjust for imprecision in single-precision real numbers exceeding the resolvable values in long (>23-year) continuous simulations.

Author: tlspero

dyn_em/module_first_rk_step_part1.F

@@ -19,6 +19,7 @@ SUBROUTINE first_rk_step_part1 (   grid , config_flags              &
                              , ph_tendf, mu_tendf               &
                              , tke_tend                         &
                              , adapt_step_flag , curr_secs      &
+                             , curr_mins2                       &
                              , psim , psih , gz1oz0 , chklowq   &
                              , cu_act_flag , hol , th_phy       &
                              , pi_phy , p_phy , t_phy           &
@@ -77,7 +78,7 @@ SUBROUTINE first_rk_step_part1 (   grid , config_flags              &
 
 
     LOGICAL ,INTENT(IN)                        :: adapt_step_flag
-    REAL, INTENT(IN)                           :: curr_secs
+    REAL, INTENT(IN)                           :: curr_secs, curr_mins2
 
     REAL    ,DIMENSION(ims:ime,kms:kme,jms:jme,num_moist),INTENT(INOUT)   :: moist
     REAL    ,DIMENSION(ims:ime,kms:kme,jms:jme,num_moist),INTENT(INOUT)   :: moist_tend
@@ -1770,7 +1771,7 @@ SUBROUTINE first_rk_step_part1 (   grid , config_flags              &
 
 BENCH_START(fdda_driver_tim)
       CALL fddagd_driver(itimestep=grid%itimestep,dt=grid%dt,xtime=grid%XTIME,      &
-                  id=grid%id,                                                       &
+                  curr_mins2=curr_mins2, id=grid%id,                                &
                   RUNDGDTEN=grid%rundgdten,RVNDGDTEN=grid%rvndgdten,                &
                   RTHNDGDTEN=grid%rthndgdten,RPHNDGDTEN=grid%rphndgdten,            &
                   RQVNDGDTEN=grid%rqvndgdten,RMUNDGDTEN=grid%rmundgdten,            &

dyn_em/solve_em.F

@@ -183,7 +183,8 @@ SUBROUTINE solve_em ( grid , config_flags  &
    LOGICAL :: leapfrog
    INTEGER :: l,kte,kk
    LOGICAL :: f_flux  ! flag for computing averaged fluxes in cu_gd
-   REAL    :: curr_secs, curr_secs2
+   REAL    :: curr_secs, curr_secs2, curr_mins2
+   REAL(8) :: curr_secs_r8, curr_secs2_r8
    INTEGER :: num_sound_steps
    INTEGER :: idex, jdex
    REAL    :: max_msft
@@ -198,6 +199,7 @@ SUBROUTINE solve_em ( grid , config_flags  &
 
    TYPE(WRFU_TimeInterval)                    :: tmpTimeInterval, tmpTimeInterval2
    REAL                                       :: real_time
+   REAL(8)                                    :: real_time_r8
    LOGICAL                                    :: adapt_step_flag
    LOGICAL                                    :: fill_w_flag
 
@@ -331,6 +333,9 @@ END SUBROUTINE CMAQ_DRIVER
    tmpTimeInterval2 = domain_get_current_time ( grid ) - domain_get_start_time ( grid )
    curr_secs  = real_time(tmpTimeInterval)
    curr_secs2 = real_time(tmpTimeInterval2)
+   curr_secs_r8  = real_time_r8(tmpTimeInterval)
+   curr_secs2_r8 = real_time_r8(tmpTimeInterval2)
+   curr_mins2 = REAL( curr_secs2_r8 / 60.0d0 )
 
    old_dt = grid%dt   ! store old time step for flux averaging code at end of RK loop
 
@@ -811,7 +816,7 @@ END SUBROUTINE CMAQ_DRIVER
                              , ph_tendf, mu_tendf               &
                              , tke_tend                         &
                              , config_flags%use_adaptive_time_step &
-                             , curr_secs                        &
+                             , curr_secs, curr_mins2            &
                              , psim , psih , gz1oz0             &
                              , chklowq                          &
                              , cu_act_flag , hol , th_phy       &

phys/module_diag_misc.F

@@ -288,7 +288,11 @@ SUBROUTINE diagnostic_output_calc(                                 &
 !-----------------------------------------------------------------
 ! Handle accumulations with buckets to prevent round-off truncation in long runs
 ! This is done every 360 minutes assuming time step fits exactly into 360 minutes
-   IF(bucket_mm .gt. 0. .AND. MOD(NINT(XTIME),360) .EQ. 0)THEN
+
+!!!~~ CURR_SECS2 is elapsed seconds since restart. It is preferred to
+!!!~~ XTIME here because XTIME goes imprecise at 2^24, just under 32 years.
+
+   IF(bucket_mm .gt. 0. .AND. MOD(NINT(CURR_SECS2),3600) .EQ. 0)THEN
 ! SET START AND END POINTS FOR TILES
 !  !$OMP PARALLEL DO   &
 !  !$OMP PRIVATE ( ij )

phys/module_fdda_psufddagd.F

@@ -9,7 +9,17 @@
 ! surfance reanalsys 
 !Reference: Alapaty et al., 2008: Development of the flux-adjusting surface
 ! data assimilation system for mesoscale models. JAMC, 47, 2331-2350
-
+!
+! Changed logic for determining next nudging time to rely on minutes elapsed
+! since restart (CURR_MINS2) rather than on minutes since initialization
+! (XTIME). The REAL variable, XTIME, will become imprecise after 1 X 2^24
+! minutes (just under 32 years of continuous simulation). Cannot remove all
+! reliance on XTIME because actual end time is in absolute minutes. Using XTIME
+! results in spectral nudging analyses ingested at the wrong times, beginning
+! 23 years and 3.5 months into a continous simulation. Purposefully not
+! modifying the ramping function because pragmatically we will not get
+! very large XTIME values in any situation where the off-ramp for nudging would
+! be used. (Tanya Spero, U.S. Environmental Protection Agency -- June 2024)
 !
 !
 MODULE module_fdda_psufddagd
@@ -18,7 +28,7 @@ MODULE module_fdda_psufddagd
 !
 !-------------------------------------------------------------------
 !
-   SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
+   SUBROUTINE fddagd(itimestep,dx,dt,xtime, curr_mins2, &
                id,analysis_interval, end_fdda_hour, &
                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
                if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
@@ -101,7 +111,7 @@ SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
    INTEGER,  INTENT(IN)   ::      if_ramping
 
    INTEGER , INTENT(IN)   ::      id
-   REAL,     INTENT(IN)   ::      DT, dx, xtime, dtramp_min
+   REAL,     INTENT(IN)   ::      DT, dx, xtime, dtramp_min, curr_mins2
 
    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
                                   ims,ime, jms,jme, kms,kme, &
@@ -243,10 +253,10 @@ SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
    ENDIF
 
    IF( analysis_interval <= 0 )CALL wrf_error_fatal('In grid FDDA, gfdda_interval_m must be > 0')
-   xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
+   xtime_old = FLOOR(curr_mins2/analysis_interval) * analysis_interval * 1.0
    xtime_new = xtime_old + analysis_interval
    IF( int4 == 1 ) THEN
-     coef = (xtime-xtime_old)/(xtime_new-xtime_old)
+     coef = (curr_mins2-xtime_old)/(xtime_new-xtime_old)
    ELSE
      coef = 1.0          ! Nudging toward a target value (*_ndg_new values)
    ENDIF
@@ -255,7 +265,7 @@ SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
 
      CALL get_wrf_debug_level( dbg_level )
 
-     IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
+     IF( curr_mins2-xtime_old < 0.5*dt/60.0 ) THEN
 
        IF( xtime < end_fdda_hour*60.0 ) THEN
          WRITE(message,'(a,i1,a,f10.3,a)') &
@@ -578,7 +588,7 @@ SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
 ! Surface Analysis Nudging
 !
    IF( grid_sfdda >= 1 ) THEN
-     CALL SFDDAGD(itimestep,dx,dt,xtime, id, &
+     CALL SFDDAGD(itimestep,dx,dt,xtime, curr_mins2, id, &
      analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
      rinblw, &
                u3d,v3d,th3d,t3d,                 &
@@ -680,7 +690,7 @@ SUBROUTINE fddagd(itimestep,dx,dt,xtime,  &
    END SUBROUTINE fddagd
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   SUBROUTINE sfddagd(itimestep,dx,dt,xtime,  &
+   SUBROUTINE sfddagd(itimestep,dx,dt,xtime, curr_mins2,  &
                id, analysis_interval_sfc, end_fdda_hour_sfc, &
                guv_sfc, gt_sfc, gq_sfc, rinblw,  &
                u3d,v3d,th3d,t3d,                 &
@@ -758,7 +768,7 @@ SUBROUTINE sfddagd(itimestep,dx,dt,xtime,  &
    INTEGER,  INTENT(IN)   ::      itimestep, analysis_interval_sfc, end_fdda_hour_sfc
 
    INTEGER , INTENT(IN)   ::      id
-   REAL,     INTENT(IN)   ::      dx,DT, xtime
+   REAL,     INTENT(IN)   ::      dx,DT, xtime, curr_mins2
 
    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
                                   ims,ime, jms,jme, kms,kme, &
@@ -862,10 +872,10 @@ SUBROUTINE sfddagd(itimestep,dx,dt,xtime,  &
    int4 = 1      ! 1: temporal ionterpolation. else: target nudging toward *_ndg_new values
 
    IF( analysis_interval_sfc <= 0 )CALL wrf_error_fatal('In grid sfc FDDA, sgfdda_interval_m must be > 0')
-   xtime_old_sfc = FLOOR(xtime/analysis_interval_sfc) * analysis_interval_sfc * 1.0
+   xtime_old_sfc = FLOOR(curr_mins2/analysis_interval_sfc) * analysis_interval_sfc * 1.0
    xtime_new_sfc = xtime_old_sfc + analysis_interval_sfc
    IF( int4 == 1 ) THEN
-     coef = (xtime-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc)  ! Temporal interpolation
+     coef = (curr_mins2-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc)  ! Temporal interpolation
    ELSE
      coef = 1.0          ! Nudging toward a target value (*_ndg_new values)
    ENDIF
@@ -874,7 +884,7 @@ SUBROUTINE sfddagd(itimestep,dx,dt,xtime,  &
 
      CALL get_wrf_debug_level( dbg_level )
 
-     IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
+     IF( curr_mins2-xtime_old_sfc < 0.5*dt/60.0 ) THEN
 
        IF( xtime < end_fdda_hour_sfc*60.0 ) THEN
          WRITE(message,'(a,i1,a,f10.3,a)') &

phys/module_fdda_spnudging.F

@@ -5,6 +5,17 @@
 ! Added capability to spectrally nudge water vapor mixing ratio, and added
 ! user-definable lid for nudging potential temperature and water vapor mixing
 ! ratio. (Tanya Spero, U.S. Environmental Protection Agency -- October 2017)
+!
+! Changed logic for determining next nudging time to rely on minutes elapsed
+! since restart (CURR_MINS2) rather than on minutes since initialization
+! (XTIME). The REAL variable, XTIME, will become imprecise after 1 X 2^24
+! minutes (just under 32 years of continuous simulation). Cannot remove all
+! reliance on XTIME because actual end time is in absolute minutes. Using XTIME
+! results in spectral nudging analyses ingested at the wrong times, beginning
+! 23 years and 3.5 months into a continous simulation. Purposefully not
+! modifying the ramping function because pragmatically we will not get
+! very large XTIME values in any situation where the off-ramp for nudging would
+! be used. (Tanya Spero, U.S. Environmental Protection Agency -- June 2024)
 
 MODULE module_fdda_spnudging
 
@@ -17,7 +28,8 @@ MODULE module_fdda_spnudging
 !
 !-------------------------------------------------------------------
 !
-   SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fdda_hour, &
+   SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,curr_mins2, &
+               id,analysis_interval, end_fdda_hour, &
                if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_ph,if_no_pbl_nudging_q,&
                if_zfac_uv, k_zfac_uv, dk_zfac_uv, & 
                if_zfac_t, k_zfac_t, dk_zfac_t, &
@@ -95,7 +107,7 @@ SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fd
    INTEGER,  INTENT(IN)   ::      xwavenum,ywavenum
 
    INTEGER , INTENT(IN)   ::      id
-   REAL,     INTENT(IN)   ::      DT, xtime, dtramp_min
+   REAL,     INTENT(IN)   ::      DT, xtime, dtramp_min, curr_mins2
 
    INTEGER,  INTENT(IN)   ::      ids,ide, jds,jde, kds,kde, &
                                   ims,ime, jms,jme, kms,kme, &
@@ -202,15 +214,15 @@ SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fd
 !  IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
 !      actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
 
-   xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
+   xtime_old = FLOOR(curr_mins2/analysis_interval) * analysis_interval * 1.0
    xtime_new = xtime_old + analysis_interval
-   coef = (xtime-xtime_old)/(xtime_new-xtime_old)
+   coef = (curr_mins2-xtime_old)/(xtime_new-xtime_old)
 
    IF ( wrf_dm_on_monitor()) THEN
 
      CALL get_wrf_debug_level( dbg_level )
 
-     IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
+     IF( curr_mins2-xtime_old < 0.5*dt/60.0 ) THEN
 
        IF( xtime < end_fdda_hour*60.0 ) THEN
          WRITE(wrf_err_message,FMT='(a,i2.2,a,f15.3,a)') &
@@ -549,7 +561,7 @@ SUBROUTINE spectral_nudging(grid,itimestep,dt,xtime,id,analysis_interval, end_fd
        tfac = 1.0
      ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
        tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
-       IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/analysis_interval
+       IF( dtramp_min > 0.0 ) coef = (curr_mins2-xtime_old+analysis_interval)/analysis_interval
      ELSE                                                     
        tfac = 0.0
      ENDIF

phys/module_fddagd_driver.F

@@ -6,7 +6,7 @@ MODULE module_fddagd_driver
 
 !------------------------------------------------------------------
    SUBROUTINE fddagd_driver(itimestep,dt,xtime,                   &
-                  id,  &
+                  curr_mins2, id,  &
                   RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,                 &
                   RPHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,               &
                   SDA_HFX, SDA_QFX,                               & !fasdas
@@ -143,7 +143,7 @@ SUBROUTINE fddagd_driver(itimestep,dt,xtime,                   &
 
    INTEGER,    INTENT(IN   )    ::     itimestep,STEPFG
 !
-   REAL,       INTENT(IN   )    ::     DT,DX,XTIME
+   REAL,       INTENT(IN   )    ::     DT,DX,XTIME, curr_mins2
 
 
 !
@@ -521,7 +521,7 @@ SUBROUTINE fddagd_driver(itimestep,dt,xtime,                   &
            ENDIF
 
            CALL FDDAGD(itimestep,dx,dt,xtime, &
-               id, &
+               curr_mins2, id, &
                config_flags%auxinput10_interval_m, &
                config_flags%auxinput10_end_h, &
                config_flags%if_no_pbl_nudging_uv, &
@@ -570,7 +570,7 @@ SUBROUTINE fddagd_driver(itimestep,dt,xtime,                   &
       CASE (SPNUDGING)
         CALL wrf_debug(100,'in SPECTRAL NUDGING scheme')
            CALL SPECTRAL_NUDGING(grid,itimestep,dt,xtime, &
-               id, &
+               curr_mins2, id, &
                config_flags%auxinput10_interval_m, &
                config_flags%auxinput10_end_h, &
                config_flags%if_no_pbl_nudging_uv, &

wrftladj/solve_em_ad.F

@@ -178,7 +178,8 @@ SUBROUTINE solve_em_ad ( grid , config_flags  &
    LOGICAL :: leapfrog
    INTEGER :: l,kte,kk
    LOGICAL :: f_flux  ! flag for computing averaged fluxes in cu_gd
-   REAL    :: curr_secs
+   REAL    :: curr_secs, curr_secs2, curr_mins2
+   REAL(8) :: curr_secs_r8, curr_secs2_r8
    INTEGER :: num_sound_steps
    INTEGER :: idex, jdex
    REAL    :: max_msft
@@ -191,8 +192,9 @@ SUBROUTINE solve_em_ad ( grid , config_flags  &
 ! urban related variables
    INTEGER :: NUM_ROOF_LAYERS, NUM_WALL_LAYERS, NUM_ROAD_LAYERS   ! urban
 
-   TYPE(WRFU_TimeInterval)                    :: tmpTimeInterval
+   TYPE(WRFU_TimeInterval)                    :: tmpTimeInterval, tmpTimeInterval2
    REAL                                       :: real_time
+   REAL(8)                                    :: real_time_r8
    LOGICAL                                    :: adapt_step_flag
    LOGICAL                                    :: fill_w_flag
 
@@ -299,7 +301,12 @@ SUBROUTINE solve_em_ad ( grid , config_flags  &
   !   calculate it here--but, this is not clean!!
   !
    tmpTimeInterval = domain_get_current_time ( grid ) - domain_get_sim_start_time ( grid )
+   tmpTimeInterval2 = domain_get_current_time ( grid ) - domain_get_start_time ( grid )
    curr_secs = real_time(tmpTimeInterval)
+   curr_secs2 = real_time(tmpTimeInterval2)
+   curr_secs_r8  = real_time_r8(tmpTimeInterval)
+   curr_secs2_r8 = real_time_r8(tmpTimeInterval2)
+   curr_mins2 = REAL( curr_secs2_r8 / 60.0d0 )
 
    old_dt = grid%dt   ! store old time step for flux averaging code at end of RK loop
 !-----------------------------------------------------------------------------
@@ -672,7 +679,7 @@ SUBROUTINE solve_em_ad ( grid , config_flags  &
                              , ph_tendf, mu_tendf               &
                              , tke_tend                         &
                              , config_flags%use_adaptive_time_step &
-                             , curr_secs                        &
+                             , curr_secs, curr_mins2            &
                              , psim , psih , gz1oz0             &
                              , chklowq                          &
                              , cu_act_flag , hol , th_phy       &