Modify FV3 code for unit testing

.github/workflows/GCC.yml

@@ -120,7 +120,7 @@ jobs:
         export CC=mpicc
         export CXX=mpicxx
         export FC=mpif90
-        cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON ${{ matrix.cmake_opts }} ${{ env.gcov_cmake }}
+        cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON -DFV3=ON ${{ matrix.cmake_opts }} ${{ env.gcov_cmake }}
         make -j2
 
     - name: run-tests

.github/workflows/docs.yml

@@ -119,7 +119,7 @@ jobs:
         export CC=mpicc
         export CXX=mpicxx
         export FC=mpif90
-        cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON ${{ matrix.cmake_opts }} -DENABLE_FV3ATM_DOCS=ON ${{ env.gcov_cmake }}
+        cmake ${GITHUB_WORKSPACE}/fv3atm -DBUILD_CI_TESTING=ON -DFV3=ON ${{ matrix.cmake_opts }} -DENABLE_FV3ATM_DOCS=ON ${{ env.gcov_cmake }}
         make doxygen_doc
 
     - name: upload-docs

ci/CMakeLists.txt

@@ -7,6 +7,15 @@
 
 list(APPEND CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMakeModules/Modules)
 
+# Switch to ESMF's cmake files. They have been removed from EMC's modules repo
+if(DEFINED ENV{ESMFMKFILE})
+  get_filename_component(TEMP_DIR $ENV{ESMFMKFILE} DIRECTORY)
+  get_filename_component(ESMF_ROOT ${TEMP_DIR} DIRECTORY)
+  list(APPEND CMAKE_MODULE_PATH ${ESMF_ROOT}/cmake)
+  include_directories(${ESMF_ROOT}/include)
+  link_directories(${TEMP_DIR})
+endif()
+
 if(${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU")
   if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 9.0.0)
     message(FATAL_ERROR "GNU Compiler >= 9 is required")

fv3/atmos_model.F90

@@ -134,6 +134,8 @@ module atmos_model_mod
 public atmos_model_get_nth_domain_info
 public addLsmask2grid
 public setup_exportdata
+public set_fhzero_loop, InitTimeFromIAUOffset
+public get_atmos_tracer_types
 !-----------------------------------------------------------------------
 
 !<PUBLICTYPE >
@@ -784,33 +786,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
    !--- set the initial diagnostic timestamp
    diag_time = Time
    call get_time (Atmos%Time - Atmos%Time_init, sec)
-   !--- Model should restart at the forecast hours that are multiples of fhzero.
-   !--- WARNING: For special cases that model needs to restart at non-multiple of fhzero
-   !--- the fields in first output files are not accumulated from the beginning of
-   !--- the bucket, but the restart time.
-   if( GFS_Control%fhzero_array(1) > 0. ) then
-     fhzero_loop: do i=1,size(GFS_Control%fhzero_array)
-       tmpflag_fhzero= .false.
-       if( GFS_Control%fhzero_array(i) > 0.) then
-         if( i == 1 ) then
-           if( sec <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true.
-         else if( i > 1 ) then
-           if( sec > GFS_Control%fhzero_fhour(i-1)*3600. .and. sec <=GFS_Control%fhzero_fhour(i)*3600. ) &
-             tmpflag_fhzero = .true.
-         endif
-         if( tmpflag_fhzero ) then
-           GFS_Control%fhzero = GFS_Control%fhzero_array(i)
-           if( GFS_Control%fhzero > 0) then
-             sec_lastfhzerofh = (int(sec/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600
-           else
-             sec_lastfhzerofh = 0
-           endif
-         endif
-       endif
-     enddo fhzero_loop
-   else
-     sec_lastfhzerofh = 0
-   endif
+   call set_fhzero_loop(sec, sec_lastfhzerofh)
    if (mpp_pe() == mpp_root_pe()) print *,'in atmos_model, fhzero=',GFS_Control%fhzero, 'fhour=',sec/3600.,sec_lastfhzerofh/3600
 
    if (mod((sec-sec_lastfhzerofh),int(GFS_Control%fhzero*3600.)) /= 0) then
@@ -858,6 +834,50 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
 end subroutine atmos_model_init
 ! </SUBROUTINE>
 
+  !> This will set the forecast hour based on the fhzero_array 
+  !> input. This should handle both an input of a full fhzero 
+  !> array and one that uses increments.
+  !>
+  !> @param tmpflag_fhzero logical if current timestep is in between output hours
+  !> @param[inout] sec time since model initialization, in sec
+  !> @param[inout] sec_lastfhzerofh time since last fhzero time, in sec
+  !>
+  !> @author Daniel Sarmiento @date May 16, 2025
+subroutine set_fhzero_loop(sec, sec_lastfhzerofh)
+
+   logical                      :: tmpflag_fhzero
+   integer                      :: i
+   integer, intent(inout)       :: sec, sec_lastfhzerofh
+
+   !--- Model should restart at the forecast hours that are multiples of fhzero.
+   !--- WARNING: For special cases that model needs to restart at non-multiple of fhzero
+   !--- the fields in first output files are not accumulated from the beginning of
+   !--- the bucket, but the restart time.
+   if( GFS_Control%fhzero_array(1) > 0. ) then
+     fhzero_loop: do i=1,size(GFS_Control%fhzero_array)
+       tmpflag_fhzero= .false.
+       if( GFS_Control%fhzero_array(i) > 0.) then
+         if( i == 1 ) then
+           if( sec <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true.
+         else if( i > 1 ) then
+           if( sec > GFS_Control%fhzero_fhour(i-1)*3600. .and. sec <=GFS_Control%fhzero_fhour(i)*3600. ) &
+             tmpflag_fhzero = .true.
+         endif
+         if( tmpflag_fhzero ) then
+           GFS_Control%fhzero = GFS_Control%fhzero_array(i)
+           if( GFS_Control%fhzero > 0) then
+             sec_lastfhzerofh = (int(sec/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600
+           else
+             sec_lastfhzerofh = 0
+           endif
+         endif
+       endif
+     enddo fhzero_loop
+   else
+     sec_lastfhzerofh = 0
+   endif
+
+end subroutine set_fhzero_loop
 
 !#######################################################################
 ! <SUBROUTINE NAME="update_atmos_model_dynamics"
@@ -1005,21 +1025,7 @@ subroutine update_atmos_model_state (Atmos, rc)
     if (ANY(nint(output_fh(:)*3600.0) == seconds) .or. (GFS_control%kdt == first_kdt)) then
       if (mpp_pe() == mpp_root_pe()) write(6,*) "---isec,seconds",isec,seconds
       time_int = real(isec)
-      if(Atmos%iau_offset > zero) then
-        if( time_int - Atmos%iau_offset*3600. > zero ) then
-          time_int = time_int - Atmos%iau_offset*3600.
-        else if(seconds == Atmos%iau_offset*3600) then
-          call get_time (Atmos%Time - diag_time_fhzero, isec_fhzero)
-          time_int = real(isec_fhzero)
-          if (mpp_pe() == mpp_root_pe()) write(6,*) "---iseczero",isec_fhzero
-        endif
-      endif
-      time_intfull = real(seconds)
-      if(Atmos%iau_offset > zero) then
-        if( time_intfull - Atmos%iau_offset*3600. > zero) then
-          time_intfull = time_intfull - Atmos%iau_offset*3600.
-        endif
-      endif
+      call InitTimeFromIAUOffset(Atmos, time_int, time_intfull, seconds)
       if (mpp_pe() == mpp_root_pe()) write(6,*) ' gfs diags time since last bucket empty: ',time_int/3600.,'hrs'
       call atmosphere_nggps_diag(Atmos%Time)
       call fv3atm_diag_output(Atmos%Time, GFS_Diag, Atm_block, GFS_control%nx, GFS_control%ny, &
@@ -1028,29 +1034,7 @@ subroutine update_atmos_model_state (Atmos, rc)
     endif
 
     !---  find current fhzero
-    if( GFS_Control%fhzero_array(1) > 0. ) then
-      fhzero_loop: do i=1,size(GFS_Control%fhzero_array)
-        tmpflag_fhzero = .false.
-        if( GFS_Control%fhzero_array(i) > 0.) then
-          if( i == 1 ) then
-            if( seconds <= GFS_Control%fhzero_fhour(i)*3600. ) tmpflag_fhzero = .true.
-          else if( i > 1 ) then
-            if( seconds > GFS_Control%fhzero_fhour(i-1)*3600. .and. seconds <= GFS_Control%fhzero_fhour(i)*3600. ) &
-              tmpflag_fhzero = .true.
-          endif
-          if( tmpflag_fhzero) then
-            GFS_Control%fhzero = GFS_Control%fhzero_array(i)
-            if( GFS_Control%fhzero > 0) then
-              sec_lastfhzerofh = (int(seconds/3600.)/int(GFS_Control%fhzero))*int(GFS_Control%fhzero)*3600
-            else
-              sec_lastfhzerofh = 0
-            endif
-          endif
-        endif
-      enddo fhzero_loop
-    else
-      sec_lastfhzerofh = 0
-    endif
+    call set_fhzero_loop(seconds,sec_lastfhzerofh)
     if (mpp_pe() == mpp_root_pe()) print *,'in atmos_model update, fhzero=',GFS_Control%fhzero, 'fhour=',seconds/3600.,sec_lastfhzerofh/3600.
 
     if (nint(GFS_Control%fhzero) > 0) then
@@ -1078,7 +1062,39 @@ subroutine update_atmos_model_state (Atmos, rc)
  end subroutine update_atmos_model_state
 ! </SUBROUTINE>
 
+  !> This will calculate time if an IAU offest has been defined
+  !> in the model configuration.
+  !>
+  !> @param[inout] atmos the main atmos model configurations 
+  !> @param[inout] time_init model initialization time
+  !> @param[inout] time_intfull model time remaining
+  !> @param seconds time since model initialization
+  !>
+  !> @author Daniel Sarmiento @date May 16, 2025
+ subroutine InitTimeFromIAUOffset(Atmos, time_int, time_intfull, seconds)
+
+   type (atmos_data_type),   intent(inout)  :: Atmos
+   real(kind=GFS_kind_phys), intent(inout)  :: time_int, time_intfull
+   integer,                  intent(inout)  :: seconds
+   integer                                  :: isec_fhzero
+
+   if(Atmos%iau_offset > zero) then
+     if( time_int - Atmos%iau_offset*3600. > zero ) then
+       time_int = time_int - Atmos%iau_offset*3600.
+     else if(seconds == Atmos%iau_offset*3600) then
+       call get_time (Atmos%Time - diag_time_fhzero, isec_fhzero)
+       time_int = real(isec_fhzero)
+       if (mpp_pe() == mpp_root_pe()) write(6,*) "---iseczero",isec_fhzero
+     endif
+   endif
+   time_intfull = real(seconds)
+   if(Atmos%iau_offset > zero) then
+     if( time_intfull - Atmos%iau_offset*3600. > zero) then
+       time_intfull = time_intfull - Atmos%iau_offset*3600.
+     endif
+   endif
 
+ end subroutine InitTimeFromIAUOffset
 
 !#######################################################################
 ! <SUBROUTINE NAME="atmos_model_end">

fv3/fv3_cap.F90

@@ -9,7 +9,8 @@
 ! 18 Apr 2017: J. Wang          set up fcst grid component and write grid components
 ! 24 Jul 2017: J. Wang          initialization and time stepping changes for coupling
 ! 02 Nov 2017: J. Wang          Use Gerhard's transferable RouteHandle
-!
+! 20 May 2025: D. Sarmiento     Handle output hour array in seperate subroutines
+! 
 
 module fv3atm_cap_mod
 
@@ -53,6 +54,7 @@ module fv3atm_cap_mod
   implicit none
   private
   public SetServices
+  public OutputHours_FrequencyInput, OutputHours_ArrayInput
 !
 !-----------------------------------------------------------------------
 !
@@ -943,27 +945,9 @@ subroutine InitializeAdvertise(gcomp, rc)
           call ESMF_ConfigGetAttribute(CF,valueList=outputfh2,label='output_fh:', &
              count=noutput_fh, rc=rc)
           if(outputfh2(2) == -1) then
-            !--- output_hf is output frequency, the second item is -1
+            !--- output_fh is output frequency, the second item is -1
             lfreq = .true.
-            nfh = 0
-            if( nfhmax>output_startfh) nfh = nint((nfhmax-output_startfh)/outputfh2(1)) + 1
-            if( nfh > 0) then
-              allocate(output_fh(nfh))
-              if( output_startfh == 0) then
-                output_fh(1) = dt_atmos/3600.
-              else
-                output_fh(1) = output_startfh
-              endif
-              do i=2,nfh
-                output_fh(i) = (i-1)*outputfh2(1) + output_startfh
-                ! Except fh000, which is the first time output, if any other of the
-                ! output time is not integer hour, set lflname_fulltime to be true, so the
-                ! history file names will contain the full time stamp (HHH-MM-SS).
-                if(.not.lflname_fulltime) then
-                  if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
-                endif
-              enddo
-            endif
+            call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2)
           endif
         endif
         if( noutput_fh /= 2 .or. .not. lfreq ) then
@@ -972,32 +956,7 @@ subroutine InitializeAdvertise(gcomp, rc)
           call ESMF_ConfigGetAttribute(CF,valueList=output_fh,label='output_fh:', &
              count=noutput_fh, rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-          if( output_startfh == 0) then
-            ! If the output time in output_fh array contains first time stamp output,
-            ! check the rest of output time, otherwise, check all the output time.
-            ! If any of them is not integer hour, the history file names will
-            ! contain the full time stamp (HHH-MM-SS)
-            ist = 1
-            if(output_fh(1)==0) then
-              output_fh(1) = dt_atmos/3600.
-              ist= 2
-            endif
-            do i=ist,noutput_fh
-              if(.not.lflname_fulltime) then
-                if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
-              endif
-            enddo
-          else
-            do i=1,noutput_fh
-              output_fh(i) = output_startfh + output_fh(i)
-              ! When output_startfh >0, check all the output time, if any of
-              ! them is not integer hour, set lflname_fulltime to be true. The
-              ! history file names will contain the full time stamp (HHH-MM-SS).
-              if(.not.lflname_fulltime) then
-                if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
-              endif
-            enddo
-          endif
+          call OutputHours_ArrayInput(noutput_fh,output_startfh)
         endif
       endif ! end loutput_fh
     endif
@@ -1030,6 +989,80 @@ subroutine InitializeAdvertise(gcomp, rc)
   end subroutine InitializeAdvertise
 
 !-----------------------------------------------------------------------------
+  !> This will calculate output hours if the user has stated a 
+  !> an fhzero frequency.
+  !>
+  !> @param[inout] nfhmax maximum number of forecast hours
+  !> @param[inout] output_startfh ouptut start time
+  !> @param[inout] outputfh2 user defined forecast hour configuration
+  !>
+  !> @author Daniel Sarmiento @date May 16, 2025
+  subroutine OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2)
+    integer                   :: nfh, i
+    real, intent(inout)       :: nfhmax, output_startfh, outputfh2(2)
+
+    nfh = 0
+    if( nfhmax>output_startfh) nfh = nint((nfhmax-output_startfh)/outputfh2(1)) + 1
+    if( nfh > 0) then
+      allocate(output_fh(nfh))
+      if( output_startfh == 0) then
+        output_fh(1) = dt_atmos/3600.
+      else
+        output_fh(1) = output_startfh
+      endif
+      do i=2,nfh
+        output_fh(i) = (i-1)*outputfh2(1) + output_startfh
+        ! Except fh000, which is the first time output, if any other of the
+        ! output time is not integer hour, set lflname_fulltime to be true, so the
+        ! history file names will contain the full time stamp (HHH-MM-SS).
+        if(.not.lflname_fulltime) then
+          if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
+        endif
+      enddo
+    endif
+  end subroutine OutputHours_FrequencyInput
+
+  !> This will calculate output hours if the user has stated a
+  !> an array of desired output hours.
+  !>
+  !> @param[inout] noutput_fh index of output hours array
+  !> @param[inout] output_startfh ouptut start time
+  !>
+  !> @author Daniel Sarmiento @date May 16, 2025
+  subroutine OutputHours_ArrayInput(noutput_fh,output_startfh)
+
+    integer                   :: ist, i
+    integer, intent(inout)    :: noutput_fh
+    real, intent(inout)       :: output_startfh
+
+    if( output_startfh == 0) then
+      ! If the output time in output_fh array contains first time stamp output,
+      ! check the rest of output time, otherwise, check all the output time.
+      ! If any of them is not integer hour, the history file names will
+      ! contain the full time stamp (HHH-MM-SS)
+      ist = 1
+      if(output_fh(1)==0) then
+        output_fh(1) = dt_atmos/3600.
+        ist= 2
+      endif
+      do i=ist,noutput_fh
+        if(.not.lflname_fulltime) then
+          if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
+        endif
+      enddo
+    else
+      do i=1,noutput_fh
+        output_fh(i) = output_startfh + output_fh(i)
+        ! When output_startfh >0, check all the output time, if any of
+        ! them is not integer hour, set lflname_fulltime to be true. The
+        ! history file names will contain the full time stamp (HHH-MM-SS).
+        if(.not.lflname_fulltime) then
+          if(mod(nint(output_fh(i)*3600.),3600) /= 0) lflname_fulltime = .true.
+        endif
+      enddo
+    endif
+
+  end subroutine OutputHours_ArrayInput
 
   subroutine InitializeRealize(gcomp, rc)
     type(ESMF_GridComp)  :: gcomp

fv3/io/module_wrt_grid_comp.F90

@@ -59,6 +59,7 @@ module module_wrt_grid_comp
 !-----------------------------------------------------------------------
 
       private
+      public get_outfile, lambert, rtll
 !
 !-----------------------------------------------------------------------
 !