Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add capability to read IAU increments in cubed sphere format #274

Closed
wants to merge 5 commits into from
Closed

Add capability to read IAU increments in cubed sphere format #274

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d16baac
Meged in netcdf_read
mark-a-potts Apr 17, 2024
b4d6609
added missing use
mark-a-potts Apr 17, 2024
554b42b
fixed checksum changes
mark-a-potts Apr 17, 2024
54990d5
removed hard-coded TESTING
mark-a-potts Apr 17, 2024
4b65db2
removed is_tracer function
mark-a-potts Apr 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 39 additions & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,21 @@ project(FV3
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)

include(GNUInstallDirs)
include(CTest)
include(FetchContent)

function(download_file url hash)

message("downloading test data")
FetchContent_Declare(download_${hash}
URL ${url}
URL_HASH SHA256=${hash}
DOWNLOAD_NO_EXTRACT true
TLS_VERIFY true
)
FetchContent_Populate(download_${hash})

endfunction(download_file)

if(NOT CMAKE_BUILD_TYPE MATCHES "^(Debug|Release|Repro|MinSizeRel|RelWithDebInfo)$")
message(STATUS "No build type specified.")
Expand Down Expand Up @@ -104,7 +119,8 @@ list(APPEND tools_srcs
tools/test_cases.F90)

list(APPEND tools_srcs_extra
tools/fv_iau_mod.F90)
tools/fv_iau_mod.F90
tools/module_get_cubed_sphere_inc.F90)

list(APPEND driver_srcs
driver/fvGFS/DYCORE_typedefs.F90
Expand Down Expand Up @@ -234,3 +250,25 @@ install(EXPORT FV3Exports
FILE fv3-targets.cmake
DESTINATION ${CONFIG_INSTALL_DESTINATION})

if(${CTESTS})
set(hash a71a77233ee56f8533df1fc5c6d551472d42165d86313004edb3b70ddb5e82c7)
download_file(
https://s3.amazonaws.com/epic.sandbox.srw/inc-test-files.tar.gz
${hash}
)
file(ARCHIVE_EXTRACT INPUT ${CMAKE_BINARY_DIR}/_deps/download_${hash}-src/inc-test-files.tar.gz DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
list(APPEND READ_CUBED_SPHERE_SRC tools/test/test_read.F90 )
list(APPEND TEST_IAU_SRC tools/test/test_iau.F90 )

add_executable(test_read ${READ_CUBED_SPHERE_SRC})
add_executable(test_iau ${TEST_IAU_SRC})

target_link_libraries(test_read esmf NetCDF::NetCDF_Fortran FMS::fms_r8 fv3 fv3ccpp ccpp_physics w3emc::w3emc_d bacio::bacio_4 sp::sp_d fv3atm fv3 )
target_link_libraries(test_iau esmf NetCDF::NetCDF_Fortran FMS::fms_r8 fv3 fv3atm fv3ccpp ccpp_physics sp::sp_d w3emc::w3emc_d bacio::bacio_4 )

target_include_directories(test_read PRIVATE ${ESMF_F90COMPILEPATHS} $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/mod )
target_include_directories(test_iau PRIVATE ${ESMF_F90COMPILEPATHS} $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/mod )

add_test(NAME incRead COMMAND test_read )
add_test(NAME iauTest COMMAND test_iau )
endif()
175 changes: 105 additions & 70 deletions tools/fv_iau_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
module fv_iau_mod

use fms2_io_mod, only: file_exists
use fms_mod, only: file_exist
use mpp_mod, only: mpp_error, FATAL, NOTE, mpp_pe
use mpp_domains_mod, only: domain2d

Expand All @@ -53,6 +54,7 @@ module fv_iau_mod
get_var1_double, &
get_var3_r4, &
get_var1_real, check_var_exists
use fv_arrays_mod, only: fv_atmos_type
#ifdef GFS_TYPES
use GFS_typedefs, only: IPD_init_type => GFS_init_type, &
IPD_control_type => GFS_control_type, &
Expand All @@ -66,6 +68,7 @@ module fv_iau_mod
use fv_treat_da_inc_mod, only: remap_coef
use tracer_manager_mod, only: get_tracer_names,get_tracer_index, get_number_tracers
use field_manager_mod, only: MODEL_ATMOS
use module_get_cubed_sphere_inc, only: read_netcdf_inc, iau_internal_data_type
implicit none

private
Expand All @@ -84,14 +87,6 @@ module fv_iau_mod
integer, allocatable :: tracer_indicies(:)

real(kind=4), allocatable:: wk3(:,:,:)
type iau_internal_data_type
real,allocatable :: ua_inc(:,:,:)
real,allocatable :: va_inc(:,:,:)
real,allocatable :: temp_inc(:,:,:)
real,allocatable :: delp_inc(:,:,:)
real,allocatable :: delz_inc(:,:,:)
real,allocatable :: tracer_inc(:,:,:,:)
end type iau_internal_data_type
type iau_external_data_type
real,allocatable :: ua_inc(:,:,:)
real,allocatable :: va_inc(:,:,:)
Expand All @@ -111,13 +106,17 @@ module fv_iau_mod
real(kind=kind_phys) :: wt_normfact
end type iau_state_type
type(iau_state_type) :: IAU_state
public iau_external_data_type,IAU_initialize,getiauforcing
public IAU_initialize,getiauforcing,iau_external_data_type

contains
subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
subroutine IAU_initialize (IPD_Control, IAU_Data, Atm, mygrid, Init_parm, testing )
!subroutine IAU_initialize (IPD_Control)
type (IPD_control_type), intent(in) :: IPD_Control
type (IAU_external_data_type), intent(inout) :: IAU_Data
type (fv_atmos_type),allocatable, intent(inout) :: Atm(:)
integer, intent(in) :: mygrid
type (IPD_init_type), intent(in) :: Init_parm
logical, optional, intent(in) :: testing
! local

character(len=128) :: fname
Expand All @@ -129,21 +128,40 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
integer:: i1, i2, j1
integer:: jbeg, jend

logical:: found
logical:: found, runTests
integer nfilesall
integer, allocatable :: idt(:)

if(present(testing)) then
runTests = testing
else
runTests = .false.
endif
is = IPD_Control%isc
ie = is + IPD_Control%nx-1
js = IPD_Control%jsc
je = js + IPD_Control%ny-1
call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)
if(runTests) then
ntracers = 4
else
call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers)
endif
allocate (tracer_names(ntracers))
allocate (tracer_indicies(ntracers))
do i = 1, ntracers
call get_tracer_names(MODEL_ATMOS, i, tracer_names(i))
tracer_indicies(i) = get_tracer_index(MODEL_ATMOS,tracer_names(i))
enddo
if(runTests) then
tracer_names(1) = "sphum"
tracer_names(2) = "liq_wat"
tracer_names(3) = "ice_wat"
tracer_names(4) = "o3mr"
tracer_indicies(1) = 1
tracer_indicies(2) = 2
tracer_indicies(3) = 3
tracer_indicies(4) = 4
else
do i = 1, ntracers
call get_tracer_names(MODEL_ATMOS, i, tracer_names(i))
tracer_indicies(i) = get_tracer_index(MODEL_ATMOS,tracer_names(i))
enddo
endif
allocate(s2c(is:ie,js:je,4))
allocate(id1(is:ie,js:je))
allocate(id2(is:ie,js:je))
Expand Down Expand Up @@ -184,60 +202,62 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
deg2rad = pi/180.

npz = IPD_Control%levs
km = npz
fname = 'INPUT/'//trim(IPD_Control%iau_inc_files(1))

if( file_exists(fname) ) then
call open_ncfile( fname, ncid ) ! open the file
call get_ncdim1( ncid, 'lon', im)
call get_ncdim1( ncid, 'lat', jm)
call get_ncdim1( ncid, 'lev', km)

if (km.ne.npz) then
if (is_master()) print *, 'km = ', km
call mpp_error(FATAL, &
'==> Error in IAU_initialize: km is not equal to npz')
if(IPD_Control%iau_gaussian) then
if( file_exist(fname) ) then
call open_ncfile( fname, ncid ) ! open the file
call get_ncdim1( ncid, 'lon', im)
call get_ncdim1( ncid, 'lat', jm)
call get_ncdim1( ncid, 'lev', km)

if (km.ne.npz) then
if (is_master()) print *, 'km = ', km
call mpp_error(FATAL, &
'==> Error in IAU_initialize: km is not equal to npz')
endif

if(is_master()) write(*,*) fname, ' DA increment dimensions:', im,jm,km

allocate ( lon(im) )
allocate ( lat(jm) )

call _GET_VAR1 (ncid, 'lon', im, lon )
call _GET_VAR1 (ncid, 'lat', jm, lat )
call close_ncfile(ncid)

! Convert to radians
do i=1,im
lon(i) = lon(i) * deg2rad
enddo
do j=1,jm
lat(j) = lat(j) * deg2rad
enddo

else
call mpp_error(FATAL,'==> Error in IAU_initialize: Expected file '&
//trim(fname)//' for DA increment does not exist')
endif

if(is_master()) write(*,*) fname, ' DA increment dimensions:', im,jm,km

allocate ( lon(im) )
allocate ( lat(jm) )

call _GET_VAR1 (ncid, 'lon', im, lon )
call _GET_VAR1 (ncid, 'lat', jm, lat )
call close_ncfile(ncid)

! Convert to radians
do i=1,im
lon(i) = lon(i) * deg2rad
enddo
do j=1,jm
lat(j) = lat(j) * deg2rad

! Initialize lat-lon to Cubed bi-linear interpolation coeff:
! populate agrid
! print*,'is,ie,js,je=',is,ie,js,ie
! print*,'size xlon=',size(Init_parm%xlon(:,1)),size(Init_parm%xlon(1,:))
! print*,'size agrid=',size(agrid(:,1,1)),size(agrid(1,:,1)),size(agrid(1,1,:))
do j = 1,size(Init_parm%xlon,2)
do i = 1,size(Init_parm%xlon,1)
! print*,i,j,is-1+j,js-1+j
agrid(is-1+i,js-1+j,1)=Init_parm%xlon(i,j)
agrid(is-1+i,js-1+j,2)=Init_parm%xlat(i,j)
enddo
enddo

else
call mpp_error(FATAL,'==> Error in IAU_initialize: Expected file '&
//trim(fname)//' for DA increment does not exist')
call remap_coef( is, ie, js, je, is, ie, js, je, &
mark-a-potts marked this conversation as resolved.
Show resolved Hide resolved
im, jm, lon, lat, id1, id2, jdc, s2c, &
agrid)
deallocate ( lon, lat,agrid )
endif

! Initialize lat-lon to Cubed bi-linear interpolation coeff:
! populate agrid
! print*,'is,ie,js,je=',is,ie,js,ie
! print*,'size xlon=',size(Init_parm%xlon(:,1)),size(Init_parm%xlon(1,:))
! print*,'size agrid=',size(agrid(:,1,1)),size(agrid(1,:,1)),size(agrid(1,1,:))
do j = 1,size(Init_parm%xlon,2)
do i = 1,size(Init_parm%xlon,1)
! print*,i,j,is-1+j,js-1+j
agrid(is-1+i,js-1+j,1)=Init_parm%xlon(i,j)
agrid(is-1+i,js-1+j,2)=Init_parm%xlat(i,j)
enddo
enddo
call remap_coef( is, ie, js, je, is, ie, js, je, &
im, jm, lon, lat, id1, id2, jdc, s2c, &
agrid)
deallocate ( lon, lat,agrid )


!mp probably need to do this
allocate(IAU_Data%ua_inc(is:ie, js:je, km))
allocate(IAU_Data%va_inc(is:ie, js:je, km))
allocate(IAU_Data%temp_inc(is:ie, js:je, km))
Expand Down Expand Up @@ -274,7 +294,12 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
enddo
iau_state%wt_normfact = (2*nstep+1)/normfact
endif
call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1)))
if(IPD_Control%iau_gaussian) then
call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1)))
else
write(6,*) 'calling read_netcdf_inc with ',trim(IPD_Control%iau_inc_files(1))
call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)),iau_state%inc1,Atm,mygrid,.false.)
endif
if (nfiles.EQ.1) then ! only need to get incrments once since constant forcing over window
call setiauforcing(IPD_Control,IAU_Data,iau_state%wt)
endif
Expand All @@ -286,18 +311,24 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
allocate (iau_state%inc2%delz_inc (is:ie, js:je, km))
allocate (iau_state%inc2%tracer_inc(is:ie, js:je, km,ntracers))
iau_state%hr2=IPD_Control%iaufhrs(2)
call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
if(IPD_Control%iau_gaussian) then
call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
else
call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)),iau_state%inc2,Atm,mygrid,.false.)
endif
endif
! print*,'in IAU init',dt,rdt
IAU_data%drymassfixer = IPD_control%iau_drymassfixer

end subroutine IAU_initialize

subroutine getiauforcing(IPD_Control,IAU_Data)
subroutine getiauforcing(IPD_Control,IAU_Data,Atm,mygrid )

implicit none
type (IPD_control_type), intent(in) :: IPD_Control
type(IAU_external_data_type), intent(inout) :: IAU_Data
type (fv_atmos_type),allocatable, intent(inout) :: Atm(:)
integer, intent(in) :: mygrid
real(kind=kind_phys) t1,t2,sx,wx,wt,dtp
integer n,i,j,k,sphum,kstep,nstep,itnext

Expand Down Expand Up @@ -371,7 +402,11 @@ subroutine getiauforcing(IPD_Control,IAU_Data)
iau_state%hr2=IPD_Control%iaufhrs(itnext)
iau_state%inc1=iau_state%inc2
if (is_master()) print *,'reading next increment file',trim(IPD_Control%iau_inc_files(itnext))
call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext)))
if(IPD_Control%iau_gaussian) then
call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext)))
else
call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)),iau_state%inc2,Atm,mygrid,.false.)
endif
endif
call updateiauforcing(IPD_Control,IAU_Data,iau_state%wt)
endif
Expand Down
2 changes: 1 addition & 1 deletion tools/fv_io.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -441,7 +441,7 @@ subroutine fv_io_read_restart(fv_domain,Atm)

allocate(pes(mpp_npes()))
call mpp_get_current_pelist(pes)

Atm(1)%flagstruct%ignore_rst_cksum = .true.
suffix = ''
fname = 'INPUT/fv_core.res.nc'
Atm(1)%Fv_restart_is_open = open_file(Atm(1)%Fv_restart,fname,"read", is_restart=.true., pelist=pes)
Expand Down
24 changes: 24 additions & 0 deletions tools/fv_netcdf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
module module_get_cubed_sphere_inc

use module_read_netcdf, only : read_netcdf_inc
use fv_iau_mod, only: iau_external_data_type

implicit none


contains

subroutine get_cubed_sphere_inc(filename, increment_data, &
testing, tests_passed,rc)
!
character(*), intent(in) :: filename
type( IAU_external_data_type ) :: increment_data
logical, intent(in) :: testing
logical, optional,intent(out) :: tests_passed
integer, optional,intent(out) :: rc
!

end subroutine get_cubed_sphere_inc

!----------------------------------------------------------------------------------------
end module module_get_cubed_sphere_inc
Loading