Adaptive time step for small grid distance

dyn_em/adapt_timestep_em.F

@@ -123,7 +123,12 @@ RECURSIVE SUBROUTINE adapt_timestep(grid, config_flags)
   !
   ! Else, calculate the time step based on cfl.
   !
-  if ( ( domain_get_advanceCount ( grid ) .EQ. 1 ) .AND. ( .NOT. config_flags%restart ) ) then
+  !BPR BEGIN
+  !At the initial time advanceCount == 0, but the following line instead looked
+  !for advanceCount == 1
+  !if ( ( domain_get_advanceCount ( grid ) .EQ. 1 ) .AND. ( .NOT. config_flags%restart ) ) then
+  if ( ( domain_get_advanceCount ( grid ) .EQ. 0 ) .AND. ( .NOT. config_flags%restart ) ) then
+  !BPR END
      if ( grid%starting_time_step_den .EQ. 0 ) then
         CALL WRFU_TimeIntervalSet(dtInterval, Sn=grid%starting_time_step, Sd=1)
      else
@@ -471,6 +476,15 @@ SUBROUTINE calc_dt(dtInterval, max_cfl, max_increase_factor, precision, &
         !
 
         factor = ( target_cfl - 0.5 * (max_cfl - target_cfl) ) / max_cfl
+
+        ! BPR BEGIN
+        ! Factor can be negative in some cases so prevent factor from being
+        ! lower than 0.1
+        ! Otherwise model crashes can occur in normalize_basetime noting that
+        ! denominator of seconds cannot be negative
+        factor = MAX(0.1,factor)
+        ! BPR END
+
         num = INT(factor * precision + 0.5)
         den = precision
 

dyn_em/start_em.F

@@ -900,6 +900,23 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
         ( ( grid%dfi_opt .EQ. DFI_NODFI ) .OR. ( grid%dfi_stage .EQ. DFI_FST ) ) ) THEN
 
 !  Calculate any variables that were not set
+!BPR BEGIN
+!     This subroutine is called more than once at the first time step for a
+!     given domain.  The following if statement is to prevent the code in
+!     the if statement from executing more than once per domain at the
+!     beginning of the model run (since last_step_update=-1 the first time
+!     this is reached and should be =0 after this).
+!     Without this if statement, when this code was executed for the second
+!     time it can result in grid%dt being set incorrectly.
+!     -This is because grid%dt will be set equal to grid%starting_time_step
+!      which ignores a possible denominator in the starting time step.
+!     -The first time this code is reached is also does this, but then the
+!      call to adapt_timestep correct this
+!     -Subsequent times this code is reached adapt_timestep will not correct
+!      this because it will recognize that it has already been executed for
+!      this timestep and exit out before doing the calculation.
+
+      if (grid%last_step_updated .NE. grid%itimestep) then
 
       if (grid%starting_time_step == -1) then
          grid%starting_time_step = NINT(4 * MIN(grid%dx,grid%dy) / 1000)
@@ -937,6 +954,8 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
       CALL wrf_dm_maxval(grid%max_msftx, idex, jdex)
       CALL wrf_dm_maxval(grid%max_msfty, idex, jdex)
 #endif
+      end if
+!BPR END
 
 !     This first call just initializes variables.
 !     If a restart, get initialized variables from restart file

external/esmf_time_f90/Meat.F90

@@ -9,6 +9,11 @@ SUBROUTINE normalize_basetime( basetime )
   USE esmf_basetimemod
   IMPLICIT NONE
   TYPE(ESMF_BaseTime), INTENT(INOUT) :: basetime
+  !BPR BEGIN
+  INTEGER(ESMF_KIND_I8) :: Sn_simplified, Sd_simplified
+  INTEGER :: primes_to_check
+  !BPR END
+
 !PRINT *,'DEBUG:  BEGIN normalize_basetime()'
   ! Consistency check...  
   IF ( basetime%Sd < 0 ) THEN
@@ -41,6 +46,30 @@ SUBROUTINE normalize_basetime( basetime )
 !PRINT *,'DEBUG:  normalize_basetime() C2:  S,Sn,Sd = ',basetime%S,basetime%Sn,basetime%Sd
     ENDIF
   ENDIF
+
+  !BPR BEGIN
+  !Simplify the fraction -- otherwise the fraction can get needlessly complicated and
+  !cause WRF to crash
+  IF ( ( basetime%Sd > 0 ) .AND. (basetime%Sn > 0 ) ) THEN
+    CALL simplify( basetime%Sn, basetime%Sd, Sn_simplified, Sd_simplified )
+    basetime%Sn = Sn_simplified
+    basetime%Sd = Sd_simplified
+    !If the numerator and denominator are both larger than 10000, after simplification
+    !using the first 9 primes, the chances increase that there is a common prime factor other
+    !than the 9 searched for in the standard simplify
+    !By only searching for more than 9 primes when the numerator and denominator are
+    !large, we avoid the additional computational expense of checking additional primes
+    !for a large number of cases
+    IF ( ( basetime%Sd > 10000 ) .AND. (basetime%Sn > 10000 ) ) THEN
+      primes_to_check = 62
+      CALL simplify_numprimes( basetime%Sn, basetime%Sd, Sn_simplified, Sd_simplified, &
+                               primes_to_check )
+      basetime%Sn = Sn_simplified
+      basetime%Sd = Sd_simplified
+    ENDIF
+  ENDIF
+  !BPR END
+
 !PRINT *,'DEBUG:  END normalize_basetime()'
 END SUBROUTINE normalize_basetime
 
@@ -754,6 +783,59 @@ SUBROUTINE simplify( ni, di, no, do )
     RETURN
 END SUBROUTINE simplify
 
+!BPR BEGIN
+! Same as simplify above, but allows user to choose the number of primes to check
+SUBROUTINE simplify_numprimes( ni, di, no, do, num_primes_to_check )
+  USE esmf_basemod
+    IMPLICIT NONE
+    INTEGER(ESMF_KIND_I8), INTENT(IN)  :: ni, di
+    INTEGER(ESMF_KIND_I8), INTENT(OUT) :: no, do
+    INTEGER, INTENT(IN) :: num_primes_to_check !Number of primes to check
+    INTEGER, PARAMETER ::  nprimes = 62
+    INTEGER(ESMF_KIND_I8), DIMENSION(nprimes), PARAMETER :: primes = (/2,3,5,7,11,13,17,&
+     19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,&
+     137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,&
+     251,257,263,269,271,277,281,283,293/)
+    INTEGER(ESMF_KIND_I8) :: pr, d, n
+    INTEGER :: np
+    LOGICAL keepgoing
+    INTEGER :: num_primes_to_check_final !Number of primes to check after being limited to max
+                                         !available number of primes
+
+    ! If the user chooses to check more primes than are currently specified in the subroutine
+    ! then use the maximum number of primes currently specified
+    num_primes_to_check_final = min(num_primes_to_check, nprimes)
+
+    IF ( ni .EQ. 0 ) THEN
+      do = 1
+      no = 0
+      RETURN
+    ENDIF
+    IF ( mod( di , ni ) .EQ. 0 ) THEN
+      do = di / ni
+      no = 1
+      RETURN
+    ENDIF
+    d = di
+    n = ni
+    DO np = 1, num_primes_to_check_final
+      pr = primes(np)
+      keepgoing = .TRUE.
+      DO WHILE ( keepgoing )
+        keepgoing = .FALSE.
+        IF ( d/pr .NE. 0 .AND. n/pr .NE. 0 .AND. MOD(d,pr) .EQ. 0 .AND. MOD(n,pr) .EQ. 0 ) THEN
+          d = d / pr
+          n = n / pr
+          keepgoing = .TRUE.
+        ENDIF
+      ENDDO
+    ENDDO
+    do = d
+    no = n
+    RETURN
+END SUBROUTINE simplify_numprimes
+!BPR END
+
 
 !$$$ this should be named "c_esmc_timesum" or something less misleading
 SUBROUTINE c_esmc_basetimesum( time1, timeinterval, timeOut )

frame/module_domain.F

@@ -768,11 +768,13 @@ SUBROUTINE alloc_and_configure_domain ( domain_id , active_this_task, grid , par
       new_grid%max_tiles   = 0
       new_grid%num_tiles_spec   = 0
       new_grid%nframes   = 0         ! initialize the number of frames per file (array assignment)
-#if (EM_CORE == 1)
-      new_grid%stepping_to_time = .FALSE.
-      new_grid%adaptation_domain = 1
-      new_grid%last_step_updated = -1
-#endif
+!BPR BEGIN
+!#if (EM_CORE == 1)
+!      new_grid%stepping_to_time = .FALSE.
+!      new_grid%adaptation_domain = 1
+!      new_grid%last_step_updated = -1
+!#endif
+!BPR BEGIN
 
 !      IF (active) THEN
         ! only allocate state if this set of tasks actually computes that domain, jm 20140822
@@ -790,6 +792,16 @@ SUBROUTINE alloc_and_configure_domain ( domain_id , active_this_task, grid , par
 !        WRITE (wrf_err_message,*)"Not allocating storage for domain ",domain_id," on this set of tasks"
 !        CALL wrf_message(TRIM(wrf_err_message))
 !      ENDIF
+
+!BPR BEGIN
+#if (EM_CORE == 1)
+!Set these here, after alloc_space_field, which initializes at least last_step_updated to zero
+      new_grid%stepping_to_time = .FALSE.
+      new_grid%adaptation_domain = 1
+      new_grid%last_step_updated = -1
+#endif
+!BPR END
+
 #if MOVE_NESTS
 !set these here, after alloc_space_field, which initializes vc_i, vc_j to zero
       new_grid%xi = -1.0