diff --git a/.github/workflows/GCC.yml b/.github/workflows/GCC.yml index 44d827db0..7e63c3c60 100644 --- a/.github/workflows/GCC.yml +++ b/.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 diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 1137d4b88..8d4a185f5 100644 --- a/.github/workflows/docs.yml +++ b/.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 diff --git a/ci/CMakeLists.txt b/ci/CMakeLists.txt index 71412b3af..43659bb48 100644 --- a/ci/CMakeLists.txt +++ b/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") diff --git a/fv3/atmos_model.F90 b/fv3/atmos_model.F90 index 159dc6f50..d5e10ef20 100644 --- a/fv3/atmos_model.F90 +++ b/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 !----------------------------------------------------------------------- ! @@ -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 ! + !> 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 !####################################################################### ! 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 ! + !> 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 !####################################################################### ! diff --git a/fv3/fv3_cap.F90 b/fv3/fv3_cap.F90 index 49f13d0d9..56a0c2ba4 100644 --- a/fv3/fv3_cap.F90 +++ b/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 diff --git a/fv3/io/module_wrt_grid_comp.F90 b/fv3/io/module_wrt_grid_comp.F90 index 7b2c1dd49..9baea5809 100644 --- a/fv3/io/module_wrt_grid_comp.F90 +++ b/fv3/io/module_wrt_grid_comp.F90 @@ -59,6 +59,7 @@ module module_wrt_grid_comp !----------------------------------------------------------------------- private + public get_outfile, lambert, rtll ! !----------------------------------------------------------------------- ! diff --git a/fv3/module_fcst_grid_comp.F90 b/fv3/module_fcst_grid_comp.F90 index aa8241ffe..cea2dfdf2 100644 --- a/fv3/module_fcst_grid_comp.F90 +++ b/fv3/module_fcst_grid_comp.F90 @@ -759,57 +759,10 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc) ! Time_step = set_time (dt_atmos,0) if (mype == 0) write(*,*)'time_init=', date_init,'time=',date,'time_end=',date_end,'dt_atmos=',dt_atmos - -! set up forecast time array that controls when to write out restart files - frestart = 0 - call get_time(Time_end - Time_init, total_inttime) -! set iau offset time - Atmos%iau_offset = iau_offset - if(iau_offset > 0 ) then - iautime = set_time(iau_offset * 3600, 0) - endif -! if the second item is -1, the first number is frequency - freq_restart = .false. - if(num_restart_interval == 2) then - if(restart_interval(2)== -1) freq_restart = .true. - endif - if(freq_restart) then - if(restart_interval(1) >= 0) then - tmpvar = restart_interval(1) * 3600 - Time_step_restart = set_time (tmpvar, 0) - if(iau_offset > 0 ) then - Time_restart = Time_init + iautime + Time_step_restart - frestart(1) = tmpvar + iau_offset *3600 - else - Time_restart = Time_init + Time_step_restart - frestart(1) = tmpvar - endif - if(restart_interval(1) > 0) then - i = 2 - do while ( Time_restart < Time_end ) - frestart(i) = frestart(i-1) + tmpvar - Time_restart = Time_restart + Time_step_restart - i = i + 1 - enddo - endif - endif -! otherwise it is an array with forecast time at which the restart files will be written out - else if(num_restart_interval >= 1) then - if(num_restart_interval == 1 .and. restart_interval(1) == 0 ) then - frestart(1) = total_inttime - else - if(iau_offset > 0 ) then - restart_starttime = iau_offset *3600 - else - restart_starttime = 0 - endif - do i=1,num_restart_interval - frestart(i) = restart_interval(i) * 3600. + restart_starttime - enddo - endif - endif -! if to write out restart at the end of forecast - if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'total_inttime=',total_inttime + + call fcst_time_array_setup(Time_init, Time_end, Time_step_restart, & + Time_restart, num_restart_interval, & + restart_interval) !------ initialize component models ------ @@ -1197,6 +1150,83 @@ subroutine create_bundle_and_add_it_to_state(name_fb, state, rc) end subroutine create_bundle_and_add_it_to_state end subroutine fcst_initialize + + !> Create forecast hour time array. This will be used + !> to dictate when restart files are going to be written. + !> + !> @param[inout] Time_init model initialization time + !> @param[inout] Time_end model end time + !> @param[inout] Time_step_restart restart time based on restart_interval + !> @param[inout] Time_restart calculated restart time + !> @param[inout] num_restart_interval user defined restart interval + !> @param[inout] restart_interval restart interval, allocatable + !> + !> @author Daniel Sarmiento @date May 16, 2025 + subroutine fcst_time_array_setup(Time_init, Time_end, Time_step_restart, & + Time_restart, num_restart_interval, & + restart_interval) + + type(time_type), intent(inout) :: Time_init, Time_end, & + Time_step_restart, & + Time_restart + type(time_type) :: iautime + integer, intent(inout) :: num_restart_interval + integer :: total_inttime, tmpvar, & + i, restart_starttime + logical :: freq_restart + real, dimension(:), allocatable, intent(inout) :: restart_interval + + ! set up forecast time array that controls when to write out restart files + frestart = 0 + call get_time(Time_end - Time_init, total_inttime) + ! set iau offset time + Atmos%iau_offset = iau_offset + if(iau_offset > 0 ) then + iautime = set_time(iau_offset * 3600, 0) + endif + ! if the second item is -1, the first number is frequency + freq_restart = .false. + if(num_restart_interval == 2) then + if(restart_interval(2)== -1) freq_restart = .true. + endif + if(freq_restart) then + if(restart_interval(1) >= 0) then + tmpvar = restart_interval(1) * 3600 + Time_step_restart = set_time (tmpvar, 0) + if(iau_offset > 0 ) then + Time_restart = Time_init + iautime + Time_step_restart + frestart(1) = tmpvar + iau_offset *3600 + else + Time_restart = Time_init + Time_step_restart + frestart(1) = tmpvar + endif + if(restart_interval(1) > 0) then + i = 2 + do while ( Time_restart < Time_end ) + frestart(i) = frestart(i-1) + tmpvar + Time_restart = Time_restart + Time_step_restart + i = i + 1 + enddo + endif + endif + ! otherwise it is an array with forecast time at which the restart files will be written out + else if(num_restart_interval >= 1) then + if(num_restart_interval == 1 .and. restart_interval(1) == 0 ) then + frestart(1) = total_inttime + else + if(iau_offset > 0 ) then + restart_starttime = iau_offset *3600 + else + restart_starttime = 0 + endif + do i=1,num_restart_interval + frestart(i) = restart_interval(i) * 3600. + restart_starttime + enddo + endif + endif + ! if to write out restart at the end of forecast + if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'total_inttime=',total_inttime + end subroutine fcst_time_array_setup ! !----------------------------------------------------------------------- !####################################################################### diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5667d3d61..1d00b4b76 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,6 +1,11 @@ # This file sets up CTest unit tests for fv3atm. # # Alex Richert, Jan 2024 +# +# Add More unit tests +# Daniel Sarmiento, Jun 2025 + +find_package(PIO REQUIRED COMPONENTS C Fortran) # Stage test data execute_process(COMMAND cmake -E create_symlink @@ -10,19 +15,19 @@ execute_process(COMMAND cmake -E create_symlink function(add_fv3atm_mpi_test TESTNAME) add_executable(${TESTNAME} ${TESTNAME}.F90) - target_link_libraries(${TESTNAME} PRIVATE fv3atm MPI::MPI_Fortran PIO::PIO_Fortran) + target_link_libraries(${TESTNAME} PRIVATE ufsatm_fv3 MPI::MPI_Fortran PIO::PIO_C) add_test(${TESTNAME} ${MPIEXEC_EXECUTABLE} -n 2 ${CMAKE_CURRENT_BINARY_DIR}/${TESTNAME}) endfunction() function(add_fv3atm_serial_test TESTNAME) add_executable(${TESTNAME} ${TESTNAME}.F90) - target_link_libraries(${TESTNAME} PRIVATE fv3atm) + target_link_libraries(${TESTNAME} PRIVATE ufsatm_fv3) add_test(${TESTNAME} ${CMAKE_CURRENT_BINARY_DIR}/${TESTNAME}) endfunction() -#foreach(testname test_write_netcdf) -# add_fv3atm_mpi_test(${testname}) -#endforeach() +foreach(testname test_atmos_model test_fv3_cap test_module_wrt_grid_comp) + add_fv3atm_mpi_test(${testname}) +endforeach() foreach(testname test_post_nems_routines) add_fv3atm_serial_test(${testname}) diff --git a/tests/test_atmos_model.F90 b/tests/test_atmos_model.F90 new file mode 100644 index 000000000..1b79838ba --- /dev/null +++ b/tests/test_atmos_model.F90 @@ -0,0 +1,192 @@ +program test_atmos_model + use atmos_model_mod, only: set_fhzero_loop, InitTimeFromIAUOffset, & + get_atmos_tracer_types, atmos_data_type + use GFS_typedefs, only: GFS_control_type, GFS_kind_phys => kind_phys + use time_manager_mod, only: time_type, set_time, get_time, operator(-) + use tracer_manager_mod, only: get_number_tracers + use field_manager_mod, only: MODEL_ATMOS + use mpp_mod, only: mpp_init, mpp_exit, FATAL, mpp_error + use CCPP_data, only: GFS_control + + implicit none + + integer :: test_passed, total_tests + integer :: suite_passed, suite_total + + ! Test Suite 1: set_fhzero_loop + call test_set_fhzero_loop_suite() + + ! Test Suite 2: get_atmos_tracer_types + call test_get_atmos_tracer_types_suite() + +contains + + !============================================================================ + ! TEST SUITE 1: set_fhzero_loop + !============================================================================ + subroutine test_set_fhzero_loop_suite() + integer :: sec, sec_lastfhzerofh + + ! Test 1: Basic functionality with single fhzero value + call test_single_fhzero() + + ! Test 2: Multiple fhzero array values + call test_multiple_fhzero() + + ! Test 3: Edge case with zero and negative values + call test_fhzero_edge_cases() + + end subroutine test_set_fhzero_loop_suite + + subroutine test_single_fhzero() + integer :: sec, sec_lastfhzerofh + + ! Setup + GFS_control%fhzero_array(1) = 6.0_GFS_kind_phys + GFS_control%fhzero_fhour(1) = 24.0_GFS_kind_phys + + ! Test case: Time within first interval + sec = 10800 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 6.0_GFS_kind_phys .or. sec_lastfhzerofh /= 0) then + print *, "Incorrect handling of single fhzero value" + stop 1 + end if + + end subroutine test_single_fhzero + + subroutine test_multiple_fhzero() + integer :: sec, sec_lastfhzerofh + + ! Setup + GFS_control%fhzero_array = [3.0_GFS_kind_phys, 6.0_GFS_kind_phys] + GFS_control%fhzero_fhour = [12.0_GFS_kind_phys, 24.0_GFS_kind_phys] + + ! Test first interval + sec = 7200 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 3.0_GFS_kind_phys) then + print *, "Incorrect handling of fhzero array" + stop 2 + end if + + ! Test second interval + sec = 50400 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (GFS_control%fhzero /= 6.0_GFS_kind_phys) then + print *, "Incorrect handling of fhzero array" + stop 3 + end if + + end subroutine test_multiple_fhzero + + subroutine test_fhzero_edge_cases() + integer :: sec, sec_lastfhzerofh + + ! Test zero fhzero value + GFS_control%fhzero_array = [0.0_GFS_kind_phys, 6.0_GFS_kind_phys] + GFS_control%fhzero_fhour = [6.0_GFS_kind_phys, 12.0_GFS_kind_phys] + + sec = 3600 + call set_fhzero_loop(sec, sec_lastfhzerofh) + + if (sec_lastfhzerofh /= 0) then + print *, "Incorrect handling of fh = 0 case" + stop 4 + end if + + end subroutine test_fhzero_edge_cases + + !============================================================================ + ! TEST SUITE 2: get_atmos_tracer_types + !============================================================================ + subroutine test_get_atmos_tracer_types_suite() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Test 1: Basic functionality with mock tracers + call test_tracer_basic_functionality() + + ! Test 2: Test with chemistry tracers + call test_chemistry_tracers() + + ! Test 3: Edge cases + call test_tracer_edge_cases() + + end subroutine test_get_atmos_tracer_types_suite + + subroutine test_tracer_basic_functionality() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! For this test, we'll simulate having 5 tracers + num_tracers = 5 + allocate(tracer_types(num_tracers)) + + ! Initialize all to zero (default) + tracer_types = 0 + + if (any(tracer_types /= 0)) then + print *, "Tracer type array being rewritten" + stop 5 + end if + + deallocate(tracer_types) + end subroutine test_tracer_basic_functionality + + subroutine test_chemistry_tracers() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Simulate having tracers with chemistry types + num_tracers = 8 + allocate(tracer_types(num_tracers)) + + ! Manually set tracer types to simulate: + ! Tracers 1-3: generic (0) + ! Tracers 4-6: chemistry prognostic (1) + ! Tracers 7-8: chemistry diagnostic (2) + tracer_types = [0, 0, 0, 1, 1, 1, 2, 2] + + ! Test generic tracers are contiguous + if (any(tracer_types(1:3) /= 0)) then + print *, "Tracer type array being rewritten or rearranged" + stop 6 + end if + + ! Test prognostic tracers are contiguous + if (any(tracer_types(4:6) /= 1)) then + print *, "Tracer type array being rewritten or rearranged" + stop 7 + end if + + ! Test diagnostic tracers are contiguous + if (any(tracer_types(7:8) /= 2)) then + print *, "Tracer type array being rewritten or rearranged" + stop 8 + end if + + deallocate(tracer_types) + end subroutine test_chemistry_tracers + + subroutine test_tracer_edge_cases() + integer, allocatable :: tracer_types(:) + integer :: num_tracers + + ! Test with large number of tracers + num_tracers = 100 + allocate(tracer_types(num_tracers)) + tracer_types = 0 + + if (size(tracer_types) /= 100) then + print *, "Tracer type array missing values when array is large" + stop 9 + end if + + deallocate(tracer_types) + end subroutine test_tracer_edge_cases + +end program test_atmos_model diff --git a/tests/test_fv3_cap.F90 b/tests/test_fv3_cap.F90 new file mode 100644 index 000000000..7d98fb571 --- /dev/null +++ b/tests/test_fv3_cap.F90 @@ -0,0 +1,229 @@ +program test_output_hours + use fv3atm_cap_mod, only: OutputHours_FrequencyInput, OutputHours_ArrayInput + use module_fv3_config, only: dt_atmos, output_fh + use module_fv3_io_def, only: lflname_fulltime + + implicit none + + ! Test variables + integer :: test_passed, test_failed + integer :: i, expected_size + logical :: test_result + + ! Variables for testing + real :: nfhmax, output_startfh, outputfh2(2) + integer :: noutput_fh + + dt_atmos = 1800 + + call test_frequency_input() + call test_array_input() + +contains + + ! Test OutputHours_FrequencyInput subroutine + subroutine test_frequency_input() + + !============================================ + ! Test 1: Basic frequency input with start at 0 + nfhmax = 24.0 + output_startfh = 0.0 + outputfh2(1) = 3.0 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + expected_size = 9 + + if (size(output_fh) /= expected_size) then + !Expected size should be 9 + print *, "Size of generated output_fh is incorrect" + stop 1 + end if + if (abs(output_fh(1) - dt_atmos/3600.0) > 1e-6) then + !First output should be dt_atmos/3600 + print *, "First output time is incorrect" + stop 2 + end if + ! lflname_fulltime should be false + if (lflname_fulltime) then + !Excluding first element, output_fh are all integers + print *, "lflname_fulltime bool set incorrectly" + stop 3 + end if + + !============================================ + ! Test 2: Frequency input with non-zero start + nfhmax = 48.0 + output_startfh = 6.0 + outputfh2(1) = 6.0 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + expected_size = 8 + + !Expected size should be 8 + if (size(output_fh) /= expected_size) then + print *, "Size of generated output_fh is incorrect" + stop 4 + end if + + ! First value (should be 6.0) + if (abs(output_fh(1) - 6.0) > 1e-6) then + print *, "First output time is incorrect" + stop 5 + end if + + ! lflname_fulltime should be false + if (lflname_fulltime) then + !Excluding first element, output_fh are all integers + print *, "lflname_fulltime bool set incorrectly" + stop 6 + end if + + !============================================ + ! Test 3: Frequency that creates non-integer hours + ! Only checking lflname_fulltime since other aspects of + ! array generation were already checked. + nfhmax = 10.0 + output_startfh = 0.0 + outputfh2(1) = 2.5 + outputfh2(2) = -1.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + ! Check lflname_fulltime, should be True since non-integer values exist + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 7 + end if + + !============================================ + ! Test 4: nfhmax equals output_startfh + nfhmax = 6.0 + output_startfh = 6.0 + outputfh2(1) = 3.0 + outputfh2(2) = -1.0 + + if (allocated(output_fh)) deallocate(output_fh) + call OutputHours_FrequencyInput(nfhmax, output_startfh, outputfh2) + + ! output_fh should not allocate when nfhmax == output_startfh + if (allocated(output_fh)) then + print *, "output_fh was allocated when output start time was equal to nfmax" + stop 8 + end if + + end subroutine test_frequency_input + + subroutine test_array_input() + + !============================================ + ! Test 1: Basic array input with start at 0 + noutput_fh = 5 + output_startfh = 0.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 3.0, 6.0, 9.0, 12.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check first value. Should be dt_atmos/3600 + if (abs(output_fh(1) - dt_atmos/3600.0) > 1e-6) then + print *, "First output time is incorrect" + stop 9 + end if + + ! Check lflname_fulltime, should be false + !Excluding first element, output_fh are all integers + if (lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 10 + end if + + !============================================ + ! Test 2: Array input with non-zero start + noutput_fh = 4 + output_startfh = 6.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 6.0, 12.0, 18.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check values (should be shifted by 6: 6, 12, 18, 24) + if (abs(output_fh(1) - 6.0) > 1e-6 .or. & + abs(output_fh(2) - 12.0) > 1e-6 .or. & + abs(output_fh(3) - 18.0) > 1e-6 .or. & + abs(output_fh(4) - 24.0) > 1e-6) then + print *, "output_fh array values were not shifted correctly or not allocated correctly" + stop 11 + end if + + ! Check lflname_fulltime (should be false) + if (lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 12 + end if + + !============================================ + ! Test 3: Array with non-integer hours + noutput_fh = 4 + output_startfh = 0.0 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 1.5, 3.0, 4.5, 6.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + test_result = .true. + print *, ' Output hours:', output_fh + + ! Check lflname_fulltime (should be true) + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 13 + end if + + !============================================ + ! Test 4: Array with fractional hours from non-zero start + noutput_fh = 3 + output_startfh = 0.5 + lflname_fulltime = .false. + + if (allocated(output_fh)) deallocate(output_fh) + allocate(output_fh(noutput_fh)) + output_fh = (/ 0.0, 3.0, 6.0 /) + + call OutputHours_ArrayInput(noutput_fh, output_startfh) + + ! Check values (should be 0.5, 3.5, 6.5) + if (abs(output_fh(1) - 0.5) > 1e-6 .or. & + abs(output_fh(2) - 3.5) > 1e-6 .or. & + abs(output_fh(3) - 6.5) > 1e-6) then + print *, "output_fh array values were not shifted correctly or not allocated correctly" + stop 14 + end if + + ! Check lflname_fulltime (should be true) + if (.not. lflname_fulltime) then + print *, "lflname_fulltime bool set incorrectly" + stop 15 + end if + + end subroutine test_array_input + +end program test_output_hours diff --git a/tests/test_module_wrt_grid_comp.F90 b/tests/test_module_wrt_grid_comp.F90 new file mode 100644 index 000000000..a3c67abad --- /dev/null +++ b/tests/test_module_wrt_grid_comp.F90 @@ -0,0 +1,181 @@ +program test_module_wrt_grid_comp + use module_wrt_grid_comp, only: get_outfile, lambert, rtll + implicit none + + call test_get_outfile() + call test_lambert() + call test_rtll() + +contains + + !--------------------------------------------------------------------------- + ! Test get_outfile subroutine + !--------------------------------------------------------------------------- + subroutine test_get_outfile() + character(len=128) :: filename(2000,3) + character(len=128) :: outfile_name(100) + integer :: noutfile + character(len=100) :: test_name + + ! Test 1: Basic test with unique filenames + filename = '' + filename(1,1) = 'file1.nc' + filename(2,1) = 'file2.nc' + filename(1,2) = 'file3.nc' + filename(1,3) = 'file4.nc' + + call get_outfile(3, filename, outfile_name, noutfile) + if ( trim(outfile_name(1)) /= "file1.nc" .or. & + trim(outfile_name(2)) /= "file2.nc" .or. & + trim(outfile_name(3)) /= "file3.nc" .or. & + trim(outfile_name(4)) /= "file4.nc") then + print *, "Filename trim function not working as intended" + stop 1 + + end if + + end subroutine test_get_outfile + + !--------------------------------------------------------------------------- + ! Test lambert subroutine + !--------------------------------------------------------------------------- + subroutine test_lambert() + real(8) :: stlat1, stlat2, c_lat, c_lon + real(8) :: glon, glat, x, y + real(8) :: glon_inv, glat_inv, x_out, y_out + real(8) :: true_x, true_y + real(8), parameter :: tol = 1.0e2 ! Difference tolerance + + ! Test 1: Forward transformation (glon,glat) -> (x,y) + ! National mall, Washington DC + stlat1 = 38.5_8 + stlat2 = 39.5_8 + c_lat = 39.0_8 + c_lon = -77.0_8 + glon = -77.0353_8 + glat = 38.8895_8 + + true_x = -3055.18 + true_y = -12286.37 + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( abs(true_x - x) > tol .or. abs(true_y - y) > tol ) then + print *, "Forward transformation incorrect" + stop 2 + end if + + ! Test inverse transformation + glon_inv = 0.0_8 + glat_inv = 0.0_8 + call lambert(stlat1, stlat2, c_lat, c_lon, glon_inv, glat_inv, x, y, -1) + + if ( (glon - glon_inv) > tol .or. (glat - glat_inv) > tol ) then + print *, "Results from glon,glat -> lon,lat -> glon,glat not identical" + stop 3 + end if + + ! Test 2: Special case where stlat1 == stlat2 + ! Italy + stlat1 = 45.5_8 + stlat2 = 45.5_8 + c_lat = 45.5_8 + c_lon = 11.0_8 + glat = 45.4642_8 + glon = 9.1900_8 + + true_x = -141149.15 + true_y = -2390.66 + + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( (true_x - x) > tol .or. (true_y - y) > tol ) then + print *, "Issue when stlat1 and stlat2 are equal" + stop 4 + end if + + ! Test 3: Point at projection center + stlat1 = 38.5_8 + stlat2 = 39.5_8 + c_lat = 39.0_8 + c_lon = -77.0_8 + glon = -77.0_8 + glat = 39.0_8 + + call lambert(stlat1, stlat2, c_lat, c_lon, glon, glat, x, y, 1) + + if ( abs(x) > 1e-3 .or. abs(y) > 1e-3 ) then + print *, "Issue when point is at the projection center" + stop 5 + end if + + end subroutine test_lambert + + !--------------------------------------------------------------------------- + ! Test rtll subroutine + !--------------------------------------------------------------------------- + subroutine test_rtll() + real(8) :: tlmd, tphd, almd, aphd, tlm0d, tph0d + real(8) :: true_almd, true_aphd + real(8), parameter :: tol = 1.0e-0 + + ! Test 1: No rotation + tlm0d = 0.0_8 + tph0d = 0.0_8 + tlmd = 0.0_8 + tphd = 0.0_8 + true_almd = 0.0_8 + true_aphd = 0.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Coordinates were rotated when no rotation was expected" + stop 6 + end if + + ! Test 2: North pole rotation + tlm0d = 0.0_8 + tph0d = 90.0_8 + tlmd = 0.0_8 + tphd = 0.0_8 + true_almd = 0.0_8 + true_aphd = 90.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Incorrect rotation when 90 deg rotation was expected" + stop 7 + end if + + ! Test 3: Central Europe + tlm0d = -170.0_8 + tph0d = 40.0_8 + tlmd = 10.0_8 + tphd = 50.0_8 + true_almd = -76.16_8 + true_aphd = 83.57_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( abs(true_almd - almd) > tol .or. abs(true_aphd - aphd) > tol) then + print *, "Incorrect rotation of coordiantes" + stop 8 + end if + + ! Test 4: Longitude wrapping + tlm0d = 0.0_8 + tph0d = 45.0_8 + tlmd = 179.0_8 + tphd = 0.0_8 + + call rtll(tlmd, tphd, almd, aphd, tlm0d, tph0d) + + if( almd > 180 .or. almd < -180) then + print *, "Incorrect rotation when longitude is at the edge of the domain" + stop 9 + end if + + end subroutine test_rtll + +end program test_module_wrt_grid_comp