From af48b75934a03d546ecac15c565b9a9fb876d0e0 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Mon, 29 Apr 2024 14:01:35 +0000 Subject: [PATCH 01/16] Initial commit --- CMakeLists.txt | 3 +- model/fv_arrays.F90 | 1 + model/fv_control.F90 | 4 +- tools/fv_iau_mod.F90 | 33 +++--- tools/fv_restart.F90 | 14 ++- tools/fv_treat_da_inc.F90 | 156 +++++++++++++++++++++++++- tools/module_get_cubed_sphere_inc.F90 | 84 ++++++++++++++ 7 files changed, 273 insertions(+), 22 deletions(-) create mode 100644 tools/module_get_cubed_sphere_inc.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 32c8075cf..c5b9e058d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,7 +101,8 @@ list(APPEND tools_srcs tools/module_diag_hailcast.F90 tools/sim_nc_mod.F90 tools/sorted_index.F90 - tools/test_cases.F90) + tools/test_cases.F90 + tools/module_get_cubed_sphere_inc.F90) list(APPEND tools_srcs_extra tools/fv_iau_mod.F90) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index b5dbf9e4f..27c933430 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -831,6 +831,7 @@ module fv_arrays_mod !< condition file (if nggps_ic or ecwmf_ic are .true.). This overrides the !< hard-coded levels in fv_eta. The default is .false. logical :: read_increment = .false. !< read in analysis increment and add to restart + logical :: gaussian_increment = .true. !< increment is on Gaussian grid ! following are namelist parameters for Stochastic Energy Baskscatter dissipation estimate logical :: do_skeb = .false. !< save dissipation estimate integer :: skeb_npass = 11 !< Filter dissipation estimate "skeb_npass" times diff --git a/model/fv_control.F90 b/model/fv_control.F90 index bf111afd2..31fd1574d 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -373,6 +373,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, logical , pointer :: external_ic logical , pointer :: external_eta logical , pointer :: read_increment + logical , pointer :: gaussian_increment logical , pointer :: hydrostatic logical , pointer :: phys_hydrostatic logical , pointer :: use_hydro_pressure @@ -949,6 +950,7 @@ subroutine set_namelist_pointers(Atm) external_ic => Atm%flagstruct%external_ic external_eta => Atm%flagstruct%external_eta read_increment => Atm%flagstruct%read_increment + gaussian_increment => Atm%flagstruct%gaussian_increment hydrostatic => Atm%flagstruct%hydrostatic phys_hydrostatic => Atm%flagstruct%phys_hydrostatic @@ -1070,7 +1072,7 @@ subroutine read_namelist_fv_core_nml(Atm) do_schmidt, do_cube_transform, & hord_mt, hord_vt, hord_tm, hord_dp, hord_tr, shift_fac, stretch_fac, target_lat, target_lon, & kord_mt, kord_wz, kord_tm, kord_tr, fv_debug, fv_timers, fv_land, nudge, do_sat_adj, do_inline_mp, do_f3d, & - external_ic, read_increment, ncep_ic, nggps_ic, hrrrv3_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, & + external_ic, read_increment, gaussian_increment, ncep_ic, nggps_ic, hrrrv3_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, & external_eta, res_latlon_dynamics, res_latlon_tracers, scale_z, w_max, z_min, lim_fac, & dddmp, d2_bg, d4_bg, vtdm4, trdm2, d_ext, delt_max, beta, non_ortho, n_sponge, & warm_start, adjust_dry_mass, mountain, d_con, ke_bg, nord, nord_tr, convert_ke, use_old_omega, & diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index e4ca56500..8e175ad11 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -66,6 +66,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 @@ -84,14 +85,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(:,:,:) @@ -114,10 +107,11 @@ module fv_iau_mod public iau_external_data_type,IAU_initialize,getiauforcing contains -subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) +subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) type (IPD_control_type), intent(in) :: IPD_Control type (IAU_external_data_type), intent(inout) :: IAU_Data type (IPD_init_type), intent(in) :: Init_parm + type (fv_atmos_type), intent(inout) :: Atm ! local character(len=128) :: fname @@ -274,7 +268,11 @@ 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 ( Atm%flagstruct%gaussian_increment ) then + call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1))) + else + call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm) + 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 @@ -286,18 +284,23 @@ 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 ( Atm%flagstruct%gaussian_increment ) 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) + 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) implicit none type (IPD_control_type), intent(in) :: IPD_Control type(IAU_external_data_type), intent(inout) :: IAU_Data + type (fv_atmos_type), intent(inout) :: Atm real(kind=kind_phys) t1,t2,sx,wx,wt,dtp integer n,i,j,k,sphum,kstep,nstep,itnext @@ -371,7 +374,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 ( Atm%flagstruct%gaussian_increment ) 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) + endif endif call updateiauforcing(IPD_Control,IAU_Data,iau_state%wt) endif diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index a3291ba91..c1ccbfa08 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -166,7 +166,7 @@ module fv_restart_mod use mpp_mod, only: mpp_send, mpp_recv, mpp_sync_self, mpp_set_current_pelist, mpp_get_current_pelist, mpp_npes, mpp_pe, mpp_sync use mpp_domains_mod, only: CENTER, CORNER, NORTH, EAST, mpp_get_C2F_index, WEST, SOUTH use mpp_domains_mod, only: mpp_global_field - use fv_treat_da_inc_mod, only: read_da_inc + use fv_treat_da_inc_mod, only: read_da_inc, read_da_inc_cubed_sphere use fms2_io_mod, only: file_exists, set_filename_appendix, FmsNetcdfFile_t, open_file, close_file use coarse_grained_restart_files_mod, only: fv_io_write_restart_coarse use fv_regional_mod, only: write_full_fields @@ -372,9 +372,15 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this j = (Atm(n)%bd%jsc + Atm(n)%bd%jec)/2 k = Atm(n)%npz/2 if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k) - call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & - Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & - isc, jsc, iec, jec ) + if ( Atm(n)%flagstruct%gaussian_increment ) then + call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & + Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & + isc, jsc, iec, jec ) + else + call read_da_inc_cubed_sphere(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & + Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & + isc, jsc, iec, jec ) + endif if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k) endif !====== end PJP added DA functionailty====== diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 4fbee842d..b7cd65d7d 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -126,12 +126,17 @@ module fv_treat_da_inc_mod use fv_grid_utils_mod, only: ptop_min, g_sum, & mid_pt_sphere, get_unit_vect2, & get_latlon_vector, inner_prod, & - cubed_to_latlon + cubed_to_latlon, & + update2d_dwinds_phys, & + update_dwinds_phys use fv_mp_mod, only: is_master, & fill_corners, & YDir, & mp_reduce_min, & - mp_reduce_max + mp_reduce_max, & + group_halo_update_type, & + start_group_halo_update, & + complete_group_halo_update use sim_nc_mod, only: open_ncfile, & close_ncfile, & get_ncdim1, & @@ -140,10 +145,12 @@ module fv_treat_da_inc_mod get_var3_r4, & get_var1_real, & check_var_exists + use module_get_cubed_sphere_inc, only: read_netcdf_inc, & + iau_internal_data_type implicit none private - public :: read_da_inc,remap_coef + public :: read_da_inc, read_da_inc_cubed_sphere, remap_coef contains !============================================================================= @@ -461,6 +468,149 @@ end subroutine apply_inc_on_3d_scalar !--------------------------------------------------------------------------- end subroutine read_da_inc + subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, u, v, q, delp, pt, delz, & + is_in, js_in, ie_in, je_in, isc_in, jsc_in, iec_in, jec_in) + type(fv_atmos_type), intent(in) :: Atm + type(domain2d), intent(inout) :: fv_domain + type(fv_grid_bounds_type), intent(in) :: bd + integer, intent(in) :: npz_in, nq, is_in, js_in, ie_in, je_in, isc_in, jsc_in, iec_in, jec_in + real, intent(inout) :: u(is_in:ie_in, js_in:je_in+1, npz_in) ! D grid zonal wind (m/s) + real, intent(inout) :: v(is_in:ie_in+1, js_in:je_in, npz_in) ! D grid meridional wind (m/s) + real, intent(inout) :: delp(is_in:ie_in, js_in:je_in, npz_in) ! pressure thickness (pascal) + real, intent(inout) :: pt( is_in:ie_in, js_in:je_in, npz_in) ! temperature (K) + real, intent(inout) :: q( is_in:ie_in, js_in:je_in, npz_in, nq) ! + real, intent(inout) :: delz(isc_in:iec_in, jsc_in:jec_in, npz_in) ! + + character(len=128) :: fname + integer :: sphum, liq_wat, spo, spo2, spo3, o3mr, i,j, k + type(iau_internal_data_type) :: increment_data + type(group_halo_update_type) :: i_pack(2) + + ! Get increment filename + fname = 'INPUT/'//Atm%flagstruct%res_latlon_dynamics + + ! Ensure file exists + if ( .not. file_exists(fname) ) then + call mpp_error(FATAL,'==> Error in read_da_inc: Expected file '& + //trim(fname)//' for DA increment does not exist') + endif + + ! Read increments + call read_netcdf_inc(fname, increment_data, Atm) + + ! Wind increments + ! --------------- + + ! Start halo update + if ( Atm%gridstruct%square_domain ) then + call start_group_halo_update(i_pack(1), increment_data%ua_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.false.) + call start_group_halo_update(i_pack(1), increment_data%va_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.true.) + else + call start_group_halo_update(i_pack(1), increment_data%ua_inc, fv_domain, complete=.false.) + call start_group_halo_update(i_pack(1), increment_data%va_inc, fv_domain, complete=.true.) + endif + + if ( Atm%flagstruct%dwind_2d ) then + ! Apply A-grid wind increments to D-grid winds for case of dwind_2d + call update2d_dwinds_phys(is_in, ie_in, js_in, je_in, bd%isd, bd%ied, bd%jsd, bd%jed, & + 1., increment_data%ua_inc, increment_data%va_inc, u, v, & + Atm%gridstruct, Atm%npx, Atm%npy, npz_in, fv_domain) + else + ! Complete halo update + call complete_group_halo_update(i_pack(1), fv_domain) + + ! Treat increment boundary conditions for case of regional grid + if (Atm%gridstruct%regional) then + ! Edges + if (is_in == 1) then + do k=1,npz_in + do j = js_in,je_in + increment_data%ua_inc(is_in-1,j,k) = increment_data%ua_inc(is_in,j,k) + increment_data%va_inc(is_in-1,j,k) = increment_data%va_inc(is_in,j,k) + enddo + enddo + endif + if (ie_in == Atm%npx) then + do k=1,npz_in + do j = js_in,je_in + increment_data%ua_inc(ie_in+1,j,k) = increment_data%ua_inc(ie_in,j,k) + increment_data%va_inc(ie_in+1,j,k) = increment_data%va_inc(ie_in,j,k) + enddo + enddo + endif + if (js_in == 1) then + do k=1,npz_in + do i = is_in,ie_in + increment_data%ua_inc(i,js_in-1,k) = increment_data%ua_inc(i,js_in,k) + increment_data%va_inc(i,js_in-1,k) = increment_data%va_inc(i,js_in,k) + enddo + enddo + endif + if (je_in == Atm%npy) then + do k=1,npz_in + do i = is_in,ie_in + increment_data%ua_inc(i,je_in+1,k) = increment_data%ua_inc(i,je_in,k) + increment_data%va_inc(i,je_in+1,k) = increment_data%va_inc(i,je_in,k) + enddo + enddo + endif + + ! corners + do k=1,npz_in + if (is_in == 1 .and. js_in == 1) then + increment_data%ua_inc(is_in-1,js_in-1,k) = increment_data%ua_inc(is_in,js_in,k) + increment_data%va_inc(is_in-1,js_in-1,k) = increment_data%va_inc(is_in,js_in,k) + elseif (is_in == 1 .and. je_in == Atm%npy) then + increment_data%ua_inc(is_in-1,je_in+1,k) = increment_data%ua_inc(is_in,je_in,k) + increment_data%va_inc(is_in-1,je_in+1,k) = increment_data%va_inc(is_in,je_in,k) + elseif (ie_in == Atm%npx .and. js_in == 1) then + increment_data%ua_inc(ie_in+1,js_in-1,k) = increment_data%ua_inc(ie_in,je_in,k) + increment_data%va_inc(ie_in+1,js_in-1,k) = increment_data%va_inc(ie_in,je_in,k) + elseif (ie_in == Atm%npx .and. je_in == Atm%npy) then + increment_data%ua_inc(ie_in+1,je_in+1,k) = increment_data%ua_inc(ie_in,je_in,k) + increment_data%va_inc(ie_in+1,je_in+1,k) = increment_data%va_inc(ie_in,je_in,k) + endif + enddo + endif + + ! Apply A-grid wind increments to D-grid winds + call update_dwinds_phys(is_in, ie_in, js_in, je_in, bd%isd, bd%ied, bd%jsd, bd%jed, & + 1., increment_data%ua_inc, increment_data%va_inc, u, v, & + Atm%gridstruct, Atm%npx, Atm%npy, npz_in, fv_domain) + endif + + ! Remaining increments + ! -------------------- + + ! Get tracer indices + sphum = get_tracer_index(MODEL_ATMOS, 'sphum') + liq_wat = get_tracer_index(MODEL_ATMOS, 'liq_wat') +#ifdef MULTI_GASES + spo = get_tracer_index(MODEL_ATMOS, 'spo') + spo2 = get_tracer_index(MODEL_ATMOS, 'spo2') + spo3 = get_tracer_index(MODEL_ATMOS, 'spo3') +#else + o3mr = get_tracer_index(MODEL_ATMOS, 'o3mr') +#endif + + ! Apply increments + pt = pt + increment_data%temp_inc + delp = delp + increment_data%delp_inc + if ( .not. Atm%flagstruct%hydrostatic ) then + delz = delz + increment_data%delz_inc + endif + q(:,:,:,sphum) = q(:,:,:,sphum) + increment_data%tracer_inc(:,:,:,sphum) + q(:,:,:,liq_wat) = q(:,:,:,liq_wat) + increment_data%tracer_inc(:,:,:,liq_wat) +#ifdef MULTI_GASES + q(:,:,:,spo) = q(:,:,:,spo) + increment_data%tracer_inc(:,:,:,spo) + q(:,:,:,spo2) = q(:,:,:,spo2) + increment_data%tracer_inc(:,:,:,spo2) + q(:,:,:,spo3) = q(:,:,:,spo3) + increment_data%tracer_inc(:,:,:,spo3) +#else + q(:,:,:,o3mr) = q(:,:,:,o3mr) + increment_data%tracer_inc(:,:,:,o3mr) +#endif + + end subroutine read_da_inc_cubed_sphere + !============================================================================= !>@brief The subroutine 'remap_coef' calculates the coefficients for horizonal regridding. diff --git a/tools/module_get_cubed_sphere_inc.F90 b/tools/module_get_cubed_sphere_inc.F90 new file mode 100644 index 000000000..0aa8838dd --- /dev/null +++ b/tools/module_get_cubed_sphere_inc.F90 @@ -0,0 +1,84 @@ +#define NC_ERR_STOP(status) \ + if (status /= nf90_noerr) write(0,*) "file: ", __FILE__, " line: ", __LINE__, trim(nf90_strerror(status)); \ + if (status /= nf90_noerr) call ESMF_Finalize(endflag=ESMF_END_ABORT) + +module module_get_cubed_sphere_inc + + use netcdf + use esmf + use mpp_mod, only: mpp_pe, mpp_get_current_pelist + use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names + use field_manager_mod, only: MODEL_ATMOS + use CCPP_data, only: GFS_control + use time_manager_mod, only: time_type, get_time, set_time + use fv_arrays_mod, only: fv_atmos_type + use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfFile_t + use mpi + + implicit none + 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 + + public read_netcdf_inc, iau_internal_data_type + + contains + +!---------------------------------------------------------------------------------------- + + subroutine read_netcdf_inc(filename, increment_data, Atm) + character(*), intent(in) :: filename + type(iau_internal_data_type), intent(inout) :: increment_data + type(fv_atmos_type), intent(in) :: Atm + + type(FmsNetcdfFile_t) :: fileobj + integer :: is, ie, js, je, km, itracer, ntracers, itile + integer, dimension(4) :: corner, edge_lengths + character(len=64) :: tracer_name + + ! Get various dimensions + call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + is = GFS_control%isc + js = GFS_control%jsc + ie = is + GFS_control%nx - 1 + je = js + GFS_Control%ny - 1 + km = GFS_Control%levs + itile = Atm%tile_of_mosaic + + ! Open increment file + if ( open_file(fileobj, trim(filename), "read", pelist=Atm%pelist) ) then + + ! + corner = (/ itile, 1, js, is /) + edge_lengths = (/ 1, km, je-js+1, ie-is+1 /) + + ! Read increments + call read_data(fileobj, 'u_inc', increment_data%ua_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'v_inc', increment_data%va_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'T_inc', increment_data%temp_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'delz_inc', increment_data%delz_inc, corner=corner, edge_lengths=edge_lengths) + + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + if ( variable_exists(fileobj, tracer_name//'_inc') ) then + call read_data(fileobj, tracer_name//'_inc', increment_data%tracer_inc(:,:,:,itracer), & + corner=corner, edge_lengths=edge_lengths) + end if + end do + + ! Close increment file + call close_file(fileobj) + + end if + + end subroutine read_netcdf_inc + +!---------------------------------------------------------------------------------------- +end module module_get_cubed_sphere_inc From 5be98f969984c774fdcddcc9d587739c805d438b Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 30 Apr 2024 16:44:16 +0000 Subject: [PATCH 02/16] Polishing up --- model/fv_control.F90 | 4 +- tools/module_get_cubed_sphere_inc.F90 | 74 ++++++++++++++++----------- 2 files changed, 47 insertions(+), 31 deletions(-) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 31fd1574d..a6fea0287 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -1072,7 +1072,7 @@ subroutine read_namelist_fv_core_nml(Atm) do_schmidt, do_cube_transform, & hord_mt, hord_vt, hord_tm, hord_dp, hord_tr, shift_fac, stretch_fac, target_lat, target_lon, & kord_mt, kord_wz, kord_tm, kord_tr, fv_debug, fv_timers, fv_land, nudge, do_sat_adj, do_inline_mp, do_f3d, & - external_ic, read_increment, gaussian_increment, ncep_ic, nggps_ic, hrrrv3_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, & + external_ic, read_increment, ncep_ic, nggps_ic, hrrrv3_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, & external_eta, res_latlon_dynamics, res_latlon_tracers, scale_z, w_max, z_min, lim_fac, & dddmp, d2_bg, d4_bg, vtdm4, trdm2, d_ext, delt_max, beta, non_ortho, n_sponge, & warm_start, adjust_dry_mass, mountain, d_con, ke_bg, nord, nord_tr, convert_ke, use_old_omega, & @@ -1093,7 +1093,7 @@ subroutine read_namelist_fv_core_nml(Atm) write_coarse_restart_files,& write_coarse_diagnostics,& write_only_coarse_intermediate_restarts, & - write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum + write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum, gaussian_increment ! Read FVCORE namelist diff --git a/tools/module_get_cubed_sphere_inc.F90 b/tools/module_get_cubed_sphere_inc.F90 index 0aa8838dd..3a9e5a3e6 100644 --- a/tools/module_get_cubed_sphere_inc.F90 +++ b/tools/module_get_cubed_sphere_inc.F90 @@ -1,7 +1,3 @@ -#define NC_ERR_STOP(status) \ - if (status /= nf90_noerr) write(0,*) "file: ", __FILE__, " line: ", __LINE__, trim(nf90_strerror(status)); \ - if (status /= nf90_noerr) call ESMF_Finalize(endflag=ESMF_END_ABORT) - module module_get_cubed_sphere_inc use netcdf @@ -17,12 +13,12 @@ module module_get_cubed_sphere_inc implicit none 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(:,:,:,:) + 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 public read_netcdf_inc, iau_internal_data_type @@ -37,38 +33,58 @@ subroutine read_netcdf_inc(filename, increment_data, Atm) type(fv_atmos_type), intent(in) :: Atm type(FmsNetcdfFile_t) :: fileobj - integer :: is, ie, js, je, km, itracer, ntracers, itile + integer :: isd, ied, jsd, jed, isc, iec, jsc, jec, ng, npz, itracer, & + ntracers, itile, i integer, dimension(4) :: corner, edge_lengths character(len=64) :: tracer_name ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - is = GFS_control%isc - js = GFS_control%jsc - ie = is + GFS_control%nx - 1 - je = js + GFS_Control%ny - 1 - km = GFS_Control%levs + isd = Atm%bd%isd + jsd = Atm%bd%jsd + ied = Atm%bd%ied + jed = Atm%bd%jed + isc = Atm%bd%isc + jsc = Atm%bd%jsc + iec = Atm%bd%iec + jec = Atm%bd%jec + ng = Atm%bd%ng + npz = Atm%npz itile = Atm%tile_of_mosaic ! Open increment file - if ( open_file(fileobj, trim(filename), "read", pelist=Atm%pelist) ) then + if ( open_file(fileobj, trim(filename), 'read', pelist=Atm%pelist) ) then + if ( .not. allocated(increment_data%ua_inc) ) then + allocate(increment_data%ua_inc(isd:ied,jsd:jed,npz)) + allocate(increment_data%va_inc(isd:ied,jsd:jed,npz)) + allocate(increment_data%temp_inc(isd:ied,jsd:jed,npz)) + allocate(increment_data%delp_inc(isd:ied,jsd:jed,npz)) + allocate(increment_data%delz_inc(isc:iec,jsc:jec,npz)) + allocate(increment_data%tracer_inc(isd:ied,jsd:jed,npz,ntracers)) + end if + ! - corner = (/ itile, 1, js, is /) - edge_lengths = (/ 1, km, je-js+1, ie-is+1 /) - + corner = (/ isd+ng, jsd+ng, 1, itile /) + edge_lengths = (/ ied-isd-2*ng+1, jed-jsd-2*ng+1, npz, 1 /) + + write(*,*) 'corner ', corner + write(*,*) 'edge_lengths ', edge_lengths + ! Read increments - call read_data(fileobj, 'u_inc', increment_data%ua_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'v_inc', increment_data%va_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'T_inc', increment_data%temp_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'delp_inc', increment_data%delp_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'delz_inc', increment_data%delz_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'u_inc', increment_data%ua_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'v_inc', increment_data%va_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'T_inc', increment_data%temp_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'delz_inc', increment_data%delz_inc(isc:iec,jsc:jec,:), corner=corner, edge_lengths=edge_lengths) ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - if ( variable_exists(fileobj, tracer_name//'_inc') ) then - call read_data(fileobj, tracer_name//'_inc', increment_data%tracer_inc(:,:,:,itracer), & + do i = 1,ntracers + call get_tracer_names(MODEL_ATMOS, i, tracer_name) + itracer = get_tracer_index(MODEL_ATMOS, trim(tracer_name)) + if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then + write(*,*) itracer, ' '//trim(tracer_name)//'_inc ' + call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(isd:ied,jsd:jed,:,itracer), & corner=corner, edge_lengths=edge_lengths) end if end do From a0cca091335b97e62a6556c66671813cd4940203 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Thu, 2 May 2024 15:39:16 +0000 Subject: [PATCH 03/16] Debugging --- tools/fv_restart.F90 | 6 +- tools/fv_treat_da_inc.F90 | 147 +++++++++++++++----------- tools/module_get_cubed_sphere_inc.F90 | 48 +++------ 3 files changed, 107 insertions(+), 94 deletions(-) diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index c1ccbfa08..7012ae8e0 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -371,17 +371,19 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this i = (Atm(n)%bd%isc + Atm(n)%bd%iec)/2 j = (Atm(n)%bd%jsc + Atm(n)%bd%jec)/2 k = Atm(n)%npz/2 - if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k) if ( Atm(n)%flagstruct%gaussian_increment ) then + if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k) call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & isc, jsc, iec, jec ) + if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k) else + if( is_master() ) write(*,*) 'Calling read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k) call read_da_inc_cubed_sphere(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & isc, jsc, iec, jec ) + if( is_master() ) write(*,*) 'Back from read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k) endif - if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k) endif !====== end PJP added DA functionailty====== endif diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index b7cd65d7d..7ad78812b 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -468,7 +468,8 @@ end subroutine apply_inc_on_3d_scalar !--------------------------------------------------------------------------- end subroutine read_da_inc - subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, u, v, q, delp, pt, delz, & + subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & + u, v, q, delp, pt, delz, & is_in, js_in, ie_in, je_in, isc_in, jsc_in, iec_in, jec_in) type(fv_atmos_type), intent(in) :: Atm type(domain2d), intent(inout) :: fv_domain @@ -478,104 +479,126 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, u, v, q, del real, intent(inout) :: v(is_in:ie_in+1, js_in:je_in, npz_in) ! D grid meridional wind (m/s) real, intent(inout) :: delp(is_in:ie_in, js_in:je_in, npz_in) ! pressure thickness (pascal) real, intent(inout) :: pt( is_in:ie_in, js_in:je_in, npz_in) ! temperature (K) - real, intent(inout) :: q( is_in:ie_in, js_in:je_in, npz_in, nq) ! - real, intent(inout) :: delz(isc_in:iec_in, jsc_in:jec_in, npz_in) ! + real, intent(inout) :: q( is_in:ie_in, js_in:je_in, npz_in, nq) ! tracers + real, intent(inout) :: delz(isc_in:iec_in, jsc_in:jec_in, npz_in) ! layer thickness character(len=128) :: fname - integer :: sphum, liq_wat, spo, spo2, spo3, o3mr, i,j, k + integer :: sphum, liq_wat, spo, spo2, spo3, o3mr, i, j, k, itracer type(iau_internal_data_type) :: increment_data type(group_halo_update_type) :: i_pack(2) + real :: ua_inc(is_in:ie_in, js_in:je_in, npz_in) + real :: va_inc(is_in:ie_in, js_in:je_in, npz_in) ! Get increment filename fname = 'INPUT/'//Atm%flagstruct%res_latlon_dynamics ! Ensure file exists if ( .not. file_exists(fname) ) then - call mpp_error(FATAL,'==> Error in read_da_inc: Expected file '& + call mpp_error(FATAL,'==> Error in read_da_inc_cubed_sphere: Expected file '& //trim(fname)//' for DA increment does not exist') endif + ! Allocate increments + allocate(increment_data%ua_inc(isc_in:iec_in,jsc_in:jec_in,npz_in)) + allocate(increment_data%va_inc(isc_in:iec_in,jsc_in:jec_in,npz_in)) + allocate(increment_data%temp_inc(isc_in:iec_in,jsc_in:jec_in,npz_in)) + allocate(increment_data%delp_inc(isc_in:iec_in,jsc_in:jec_in,npz_in)) + if ( .not. Atm%flagstruct%hydrostatic ) then + allocate(increment_data%delz_inc(isc_in:iec_in,jsc_in:jec_in,npz_in)) + endif + allocate(increment_data%tracer_inc(isc_in:iec_in,jsc_in:jec_in,npz_in,nq)) + ! Read increments call read_netcdf_inc(fname, increment_data, Atm) ! Wind increments ! --------------- + ! Put u and v increments on grid that includes halo points + do k = 1,npz_in + do j = jsc_in,jec_in + do i = isc_in,iec_in + ua_inc(i,j,k) = increment_data%ua_inc(i,j,k) + va_inc(i,j,k) = increment_data%va_inc(i,j,k) + enddo + enddo + enddo + ! Start halo update if ( Atm%gridstruct%square_domain ) then - call start_group_halo_update(i_pack(1), increment_data%ua_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.false.) - call start_group_halo_update(i_pack(1), increment_data%va_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.true.) + call start_group_halo_update(i_pack(1), ua_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.false.) + call start_group_halo_update(i_pack(1), va_inc, fv_domain, whalo=1, ehalo=1, shalo=1, nhalo=1, complete=.true.) else - call start_group_halo_update(i_pack(1), increment_data%ua_inc, fv_domain, complete=.false.) - call start_group_halo_update(i_pack(1), increment_data%va_inc, fv_domain, complete=.true.) + call start_group_halo_update(i_pack(1), ua_inc, fv_domain, complete=.false.) + call start_group_halo_update(i_pack(1), va_inc, fv_domain, complete=.true.) endif if ( Atm%flagstruct%dwind_2d ) then ! Apply A-grid wind increments to D-grid winds for case of dwind_2d - call update2d_dwinds_phys(is_in, ie_in, js_in, je_in, bd%isd, bd%ied, bd%jsd, bd%jed, & - 1., increment_data%ua_inc, increment_data%va_inc, u, v, & + call update2d_dwinds_phys(isc_in, iec_in, jsc_in, jec_in, is_in, ie_in, js_in, je_in, & + 1., ua_inc, va_inc, u, v, & Atm%gridstruct, Atm%npx, Atm%npy, npz_in, fv_domain) else ! Complete halo update call complete_group_halo_update(i_pack(1), fv_domain) - ! Treat increment boundary conditions for case of regional grid - if (Atm%gridstruct%regional) then + ! Treat increment boundary conditions for case of regional grid + if ( Atm%gridstruct%regional ) then ! Edges - if (is_in == 1) then - do k=1,npz_in - do j = js_in,je_in - increment_data%ua_inc(is_in-1,j,k) = increment_data%ua_inc(is_in,j,k) - increment_data%va_inc(is_in-1,j,k) = increment_data%va_inc(is_in,j,k) + if ( isc_in == 1 ) then + do k = 1,npz_in + do j = jsc_in,jec_in + ua_inc(isc_in-1,j,k) = ua_inc(isc_in,j,k) + va_inc(isc_in-1,j,k) = va_inc(isc_in,j,k) enddo enddo endif - if (ie_in == Atm%npx) then - do k=1,npz_in - do j = js_in,je_in - increment_data%ua_inc(ie_in+1,j,k) = increment_data%ua_inc(ie_in,j,k) - increment_data%va_inc(ie_in+1,j,k) = increment_data%va_inc(ie_in,j,k) + if ( iec_in == Atm%npx ) then + do k = 1,npz_in + do j = jsc_in,jec_in + ua_inc(iec_in+1,j,k) = ua_inc(iec_in,j,k) + va_inc(iec_in+1,j,k) = va_inc(iec_in,j,k) enddo enddo endif - if (js_in == 1) then - do k=1,npz_in - do i = is_in,ie_in - increment_data%ua_inc(i,js_in-1,k) = increment_data%ua_inc(i,js_in,k) - increment_data%va_inc(i,js_in-1,k) = increment_data%va_inc(i,js_in,k) + if ( jsc_in == 1 ) then + do k = 1,npz_in + do i = isc_in,iec_in + ua_inc(i,jsc_in-1,k) = ua_inc(i,jsc_in,k) + va_inc(i,jsc_in-1,k) = va_inc(i,jsc_in,k) enddo enddo endif - if (je_in == Atm%npy) then - do k=1,npz_in - do i = is_in,ie_in - increment_data%ua_inc(i,je_in+1,k) = increment_data%ua_inc(i,je_in,k) - increment_data%va_inc(i,je_in+1,k) = increment_data%va_inc(i,je_in,k) + if ( jec_in == Atm%npy ) then + do k = 1,npz_in + do i = isc_in,iec_in + ua_inc(i,jec_in+1,k) = ua_inc(i,jec_in,k) + va_inc(i,jec_in+1,k) = va_inc(i,jec_in,k) enddo enddo endif ! corners - do k=1,npz_in - if (is_in == 1 .and. js_in == 1) then - increment_data%ua_inc(is_in-1,js_in-1,k) = increment_data%ua_inc(is_in,js_in,k) - increment_data%va_inc(is_in-1,js_in-1,k) = increment_data%va_inc(is_in,js_in,k) - elseif (is_in == 1 .and. je_in == Atm%npy) then - increment_data%ua_inc(is_in-1,je_in+1,k) = increment_data%ua_inc(is_in,je_in,k) - increment_data%va_inc(is_in-1,je_in+1,k) = increment_data%va_inc(is_in,je_in,k) - elseif (ie_in == Atm%npx .and. js_in == 1) then - increment_data%ua_inc(ie_in+1,js_in-1,k) = increment_data%ua_inc(ie_in,je_in,k) - increment_data%va_inc(ie_in+1,js_in-1,k) = increment_data%va_inc(ie_in,je_in,k) - elseif (ie_in == Atm%npx .and. je_in == Atm%npy) then - increment_data%ua_inc(ie_in+1,je_in+1,k) = increment_data%ua_inc(ie_in,je_in,k) - increment_data%va_inc(ie_in+1,je_in+1,k) = increment_data%va_inc(ie_in,je_in,k) + do k = 1,npz_in + if ( isc_in == 1 .and. jsc_in == 1 ) then + ua_inc(isc_in-1,jsc_in-1,k) = ua_inc(isc_in,jsc_in,k) + va_inc(isc_in-1,jsc_in-1,k) = va_inc(isc_in,jsc_in,k) + elseif ( isc_in == 1 .and. jec_in == Atm%npy ) then + ua_inc(isc_in-1,jec_in+1,k) = ua_inc(isc_in,jec_in,k) + va_inc(isc_in-1,jec_in+1,k) = va_inc(isc_in,jec_in,k) + elseif ( iec_in == Atm%npx .and. jsc_in == 1 ) then + ua_inc(iec_in+1,jsc_in-1,k) = ua_inc(iec_in,jec_in,k) + va_inc(iec_in+1,jsc_in-1,k) = va_inc(iec_in,jec_in,k) + elseif ( iec_in == Atm%npx .and. jec_in == Atm%npy ) then + ua_inc(iec_in+1,jec_in+1,k) = ua_inc(iec_in,jec_in,k) + va_inc(iec_in+1,jec_in+1,k) = va_inc(iec_in,jec_in,k) endif enddo endif ! Apply A-grid wind increments to D-grid winds - call update_dwinds_phys(is_in, ie_in, js_in, je_in, bd%isd, bd%ied, bd%jsd, bd%jed, & - 1., increment_data%ua_inc, increment_data%va_inc, u, v, & + call update_dwinds_phys(isc_in, iec_in, jsc_in, jec_in, is_in, ie_in, js_in, je_in, & + 1., ua_inc, va_inc, u, v, & Atm%gridstruct, Atm%npx, Atm%npy, npz_in, fv_domain) endif @@ -594,20 +617,26 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, u, v, q, del #endif ! Apply increments - pt = pt + increment_data%temp_inc - delp = delp + increment_data%delp_inc - if ( .not. Atm%flagstruct%hydrostatic ) then - delz = delz + increment_data%delz_inc - endif - q(:,:,:,sphum) = q(:,:,:,sphum) + increment_data%tracer_inc(:,:,:,sphum) - q(:,:,:,liq_wat) = q(:,:,:,liq_wat) + increment_data%tracer_inc(:,:,:,liq_wat) + do k = 1,npz_in + do j = jsc_in,jec_in + do i = isc_in,iec_in + pt(i,j,k) = pt(i,j,k) + increment_data%temp_inc(i,j,k) + delp(i,j,k) = delp(i,j,k) + increment_data%delp_inc(i,j,k) + if ( .not. Atm%flagstruct%hydrostatic ) then + delz(i,j,k) = delz(i,j,k) + increment_data%delz_inc(i,j,k) + endif + q(i,j,k,sphum) = q(i,j,k,sphum) + increment_data%tracer_inc(i,j,k,sphum) + q(i,j,k,liq_wat) = q(i,j,k,liq_wat) + increment_data%tracer_inc(i,j,k,liq_wat) #ifdef MULTI_GASES - q(:,:,:,spo) = q(:,:,:,spo) + increment_data%tracer_inc(:,:,:,spo) - q(:,:,:,spo2) = q(:,:,:,spo2) + increment_data%tracer_inc(:,:,:,spo2) - q(:,:,:,spo3) = q(:,:,:,spo3) + increment_data%tracer_inc(:,:,:,spo3) + q(i,j,k,spo) = q(i,j,k,spo) + increment_data%tracer_inc(i,j,k,spo) + q(i,j,k,spo2) = q(i,j,k,spo2) + increment_data%tracer_inc(i,j,k,spo2) + q(i,j,k,spo3) = q(i,j,k,spo3) + increment_data%tracer_inc(i,j,k,spo3) #else - q(:,:,:,o3mr) = q(:,:,:,o3mr) + increment_data%tracer_inc(:,:,:,o3mr) + q(i,j,k,o3mr) = q(i,j,k,o3mr) + increment_data%tracer_inc(i,j,k,o3mr) #endif + enddo + enddo + enddo end subroutine read_da_inc_cubed_sphere diff --git a/tools/module_get_cubed_sphere_inc.F90 b/tools/module_get_cubed_sphere_inc.F90 index 3a9e5a3e6..5d3e6b1e4 100644 --- a/tools/module_get_cubed_sphere_inc.F90 +++ b/tools/module_get_cubed_sphere_inc.F90 @@ -33,58 +33,40 @@ subroutine read_netcdf_inc(filename, increment_data, Atm) type(fv_atmos_type), intent(in) :: Atm type(FmsNetcdfFile_t) :: fileobj - integer :: isd, ied, jsd, jed, isc, iec, jsc, jec, ng, npz, itracer, & - ntracers, itile, i + integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile integer, dimension(4) :: corner, edge_lengths character(len=64) :: tracer_name ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - isd = Atm%bd%isd - jsd = Atm%bd%jsd - ied = Atm%bd%ied - jed = Atm%bd%jed isc = Atm%bd%isc - jsc = Atm%bd%jsc iec = Atm%bd%iec + jsc = Atm%bd%jsc jec = Atm%bd%jec - ng = Atm%bd%ng npz = Atm%npz itile = Atm%tile_of_mosaic - + ! Open increment file if ( open_file(fileobj, trim(filename), 'read', pelist=Atm%pelist) ) then - - if ( .not. allocated(increment_data%ua_inc) ) then - allocate(increment_data%ua_inc(isd:ied,jsd:jed,npz)) - allocate(increment_data%va_inc(isd:ied,jsd:jed,npz)) - allocate(increment_data%temp_inc(isd:ied,jsd:jed,npz)) - allocate(increment_data%delp_inc(isd:ied,jsd:jed,npz)) - allocate(increment_data%delz_inc(isc:iec,jsc:jec,npz)) - allocate(increment_data%tracer_inc(isd:ied,jsd:jed,npz,ntracers)) - end if ! - corner = (/ isd+ng, jsd+ng, 1, itile /) - edge_lengths = (/ ied-isd-2*ng+1, jed-jsd-2*ng+1, npz, 1 /) + corner = (/ isc, jsc, 1, itile /) + edge_lengths = (/ iec-isc+1, jec-jsc+1, npz, 1 /) - write(*,*) 'corner ', corner - write(*,*) 'edge_lengths ', edge_lengths - ! Read increments - call read_data(fileobj, 'u_inc', increment_data%ua_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'v_inc', increment_data%va_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'T_inc', increment_data%temp_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'delp_inc', increment_data%delp_inc(isd:ied,jsd:jed,:), corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'delz_inc', increment_data%delz_inc(isc:iec,jsc:jec,:), corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'u_inc', increment_data%ua_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'v_inc', increment_data%va_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'T_inc', increment_data%temp_inc, corner=corner, edge_lengths=edge_lengths) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc, corner=corner, edge_lengths=edge_lengths) + if ( .not. Atm%flagstruct%hydrostatic ) then + call read_data(fileobj, 'delz_inc', increment_data%delz_inc, corner=corner, edge_lengths=edge_lengths) + endif ! Read tracer increments - do i = 1,ntracers - call get_tracer_names(MODEL_ATMOS, i, tracer_name) - itracer = get_tracer_index(MODEL_ATMOS, trim(tracer_name)) + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then - write(*,*) itracer, ' '//trim(tracer_name)//'_inc ' - call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(isd:ied,jsd:jed,:,itracer), & + call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer), & corner=corner, edge_lengths=edge_lengths) end if end do From 0ffe8b313325f8a1f1e7a3ae73759ef6ff3cccc7 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 14 May 2024 20:07:58 +0000 Subject: [PATCH 04/16] multiple updates --- CMakeLists.txt | 2 +- model/fv_arrays.F90 | 2 +- model/fv_control.F90 | 7 ++-- tools/fv_iau_mod.F90 | 26 ++++++------ tools/fv_restart.F90 | 14 +++---- tools/fv_treat_da_inc.F90 | 8 ++-- tools/module_cubed_sphere_inc.F90 | 69 +++++++++++++++++++++++++++++++ tools/sim_nc_mod.F90 | 62 ++++++++++++++++++++++++++- 8 files changed, 159 insertions(+), 31 deletions(-) create mode 100644 tools/module_cubed_sphere_inc.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index c5b9e058d..478509e6b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -102,7 +102,7 @@ list(APPEND tools_srcs tools/sim_nc_mod.F90 tools/sorted_index.F90 tools/test_cases.F90 - tools/module_get_cubed_sphere_inc.F90) + tools/module_cubed_sphere_inc.F90) list(APPEND tools_srcs_extra tools/fv_iau_mod.F90) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 27c933430..003034727 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -831,7 +831,7 @@ module fv_arrays_mod !< condition file (if nggps_ic or ecwmf_ic are .true.). This overrides the !< hard-coded levels in fv_eta. The default is .false. logical :: read_increment = .false. !< read in analysis increment and add to restart - logical :: gaussian_increment = .true. !< increment is on Gaussian grid + logical :: increment_file_on_native_grid = .false. !< increment is on native cubed sphere grid grid else on Gaussian grid ! following are namelist parameters for Stochastic Energy Baskscatter dissipation estimate logical :: do_skeb = .false. !< save dissipation estimate integer :: skeb_npass = 11 !< Filter dissipation estimate "skeb_npass" times diff --git a/model/fv_control.F90 b/model/fv_control.F90 index a6fea0287..ae7697885 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -373,7 +373,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, logical , pointer :: external_ic logical , pointer :: external_eta logical , pointer :: read_increment - logical , pointer :: gaussian_increment + logical , pointer :: increment_file_on_native_grid logical , pointer :: hydrostatic logical , pointer :: phys_hydrostatic logical , pointer :: use_hydro_pressure @@ -950,8 +950,7 @@ subroutine set_namelist_pointers(Atm) external_ic => Atm%flagstruct%external_ic external_eta => Atm%flagstruct%external_eta read_increment => Atm%flagstruct%read_increment - gaussian_increment => Atm%flagstruct%gaussian_increment - + increment_file_on_native_grid => Atm%flagstruct%increment_file_on_native_grid hydrostatic => Atm%flagstruct%hydrostatic phys_hydrostatic => Atm%flagstruct%phys_hydrostatic use_hydro_pressure => Atm%flagstruct%use_hydro_pressure @@ -1093,7 +1092,7 @@ subroutine read_namelist_fv_core_nml(Atm) write_coarse_restart_files,& write_coarse_diagnostics,& write_only_coarse_intermediate_restarts, & - write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum, gaussian_increment + write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum, increment_file_on_native_grid ! Read FVCORE namelist diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 8e175ad11..a91c91d89 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -66,7 +66,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 + use module_cubed_sphere_inc, only : read_cubed_sphere_inc, increment_data_type implicit none private @@ -96,8 +96,8 @@ module fv_iau_mod logical :: drymassfixer = .false. end type iau_external_data_type type iau_state_type - type(iau_internal_data_type):: inc1 - type(iau_internal_data_type):: inc2 + type(increment_data_type) :: inc1 + type(increment_data_type) :: inc2 real(kind=kind_phys) :: hr1 real(kind=kind_phys) :: hr2 real(kind=kind_phys) :: wt @@ -268,10 +268,10 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) enddo iau_state%wt_normfact = (2*nstep+1)/normfact endif - if ( Atm%flagstruct%gaussian_increment ) then - call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1))) + if ( Atm%flagstruct%increment_file_on_native_grid ) then + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm) else - call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm) + call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1))) 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) @@ -284,10 +284,10 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) 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) - if ( Atm%flagstruct%gaussian_increment ) then - call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2))) + if ( Atm%flagstruct%increment_file_on_native_grid ) then + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)), iau_state%inc2, Atm) else - call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)), iau_state%inc2, Atm) + call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2))) endif endif ! print*,'in IAU init',dt,rdt @@ -374,10 +374,10 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) 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)) - if ( Atm%flagstruct%gaussian_increment ) then - call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext))) + if ( Atm%flagstruct%increment_file_on_native_grid ) then + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)), iau_state%inc2, Atm) else - call read_netcdf_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)), iau_state%inc2, Atm) + call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext))) endif endif call updateiauforcing(IPD_Control,IAU_Data,iau_state%wt) @@ -441,7 +441,7 @@ end subroutine setiauforcing subroutine read_iau_forcing(IPD_Control,increments,fname) type (IPD_control_type), intent(in) :: IPD_Control - type(iau_internal_data_type), intent(inout):: increments + type(increment_data_type), intent(inout):: increments character(len=*), intent(in) :: fname !locals real, dimension(:,:,:), allocatable:: u_inc, v_inc diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index 7012ae8e0..6e289b941 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -371,18 +371,18 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this i = (Atm(n)%bd%isc + Atm(n)%bd%iec)/2 j = (Atm(n)%bd%jsc + Atm(n)%bd%jec)/2 k = Atm(n)%npz/2 - if ( Atm(n)%flagstruct%gaussian_increment ) then - if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k) - call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & - Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & - isc, jsc, iec, jec ) - if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k) - else + if ( Atm(n)%flagstruct%increment_file_on_native_grid ) then if( is_master() ) write(*,*) 'Calling read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k) call read_da_inc_cubed_sphere(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & isc, jsc, iec, jec ) if( is_master() ) write(*,*) 'Back from read_da_inc_cubed_sphere',Atm(n)%pt(i,j,k) + else + if( is_master() ) write(*,*) 'Calling read_da_inc',Atm(n)%pt(i,j,k) + call read_da_inc(Atm(n), Atm(n)%domain, Atm(n)%bd, Atm(n)%npz, Atm(n)%ncnst, & + Atm(n)%u, Atm(n)%v, Atm(n)%q, Atm(n)%delp, Atm(n)%pt, Atm(n)%delz, isd, jsd, ied, jed, & + isc, jsc, iec, jec ) + if( is_master() ) write(*,*) 'Back from read_da_inc',Atm(n)%pt(i,j,k) endif endif !====== end PJP added DA functionailty====== diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 7ad78812b..e7bd2795b 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -145,8 +145,8 @@ module fv_treat_da_inc_mod get_var3_r4, & get_var1_real, & check_var_exists - use module_get_cubed_sphere_inc, only: read_netcdf_inc, & - iau_internal_data_type + use module_cubed_sphere_inc, only: read_cubed_sphere_inc, & + increment_data_type implicit none private @@ -484,7 +484,7 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & character(len=128) :: fname integer :: sphum, liq_wat, spo, spo2, spo3, o3mr, i, j, k, itracer - type(iau_internal_data_type) :: increment_data + type(increment_data_type) :: increment_data type(group_halo_update_type) :: i_pack(2) real :: ua_inc(is_in:ie_in, js_in:je_in, npz_in) real :: va_inc(is_in:ie_in, js_in:je_in, npz_in) @@ -509,7 +509,7 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & allocate(increment_data%tracer_inc(isc_in:iec_in,jsc_in:jec_in,npz_in,nq)) ! Read increments - call read_netcdf_inc(fname, increment_data, Atm) + call read_cubed_sphere_inc(fname, increment_data, Atm) ! Wind increments ! --------------- diff --git a/tools/module_cubed_sphere_inc.F90 b/tools/module_cubed_sphere_inc.F90 new file mode 100644 index 000000000..f44205d0c --- /dev/null +++ b/tools/module_cubed_sphere_inc.F90 @@ -0,0 +1,69 @@ +module module_cubed_sphere_inc + + use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names + use field_manager_mod, only: MODEL_ATMOS + use fv_arrays_mod, only: fv_atmos_type + use sim_nc_mod, only: open_ncfile, get_var5_r4, check_var_exists, close_ncfile + + implicit none + type increment_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 increment_data_type + + public read_cubed_sphere_inc, increment_data_type + + contains + +!---------------------------------------------------------------------------------------- + + subroutine read_cubed_sphere_inc(filename, increment_data, Atm) + character(*), intent(in) :: filename + type(increment_data_type), intent(inout) :: increment_data + type(fv_atmos_type), intent(in) :: Atm + + integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile + character(len=64) :: tracer_name + integer :: ncid, ncerr + + ! Get various dimensions + call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + isc = Atm%bd%isc + iec = Atm%bd%iec + jsc = Atm%bd%jsc + jec = Atm%bd%jec + npz = Atm%npz + itile = Atm%tile_of_mosaic + + ! Open increment file + call open_ncfile( trim(filename), ncid ) + + ! Read increments + call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc ) + call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc ) + call get_var5_r4( ncid, 'temp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) + call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc ) + if ( .not. Atm%flagstruct%hydrostatic ) then + call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc ) + end if + + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + call check_var_exists(ncid, trim(tracer_name)//'_inc', ncerr) + if ( ncerr == 0 ) then + call get_var5_r4( ncid, trim(tracer_name)//'_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%tracer_inc(:,:,:,itracer) ) + end if + end do + + ! Close increment file + call close_ncfile(ncid) + + end subroutine read_cubed_sphere_inc + +!---------------------------------------------------------------------------------------- +end module module_cubed_sphere_inc diff --git a/tools/sim_nc_mod.F90 b/tools/sim_nc_mod.F90 index b461aba45..19ce370c4 100644 --- a/tools/sim_nc_mod.F90 +++ b/tools/sim_nc_mod.F90 @@ -47,7 +47,7 @@ module sim_nc_mod public open_ncfile, close_ncfile, get_ncdim1, get_var1_double, get_var2_double, & get_var3_real, get_var3_double, get_var3_r4, get_var2_real, get_var2_r4, & handle_err, check_var, get_var1_real, get_var_att_double, & - check_var_exists + check_var_exists, get_var5_r4 contains @@ -339,6 +339,66 @@ subroutine get_var4_double( ncid, var4_name, im, jm, km, nt, var4 ) if (status .ne. NF_NOERR) call handle_err('get_var4_double get_vara_double '//var4_name,status) end subroutine get_var4_double + + subroutine get_var5_r4( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) + implicit none +#include + integer, intent(in):: ncid + character*(*), intent(in):: var5_name + integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime + real*4, intent(out):: var5(is:ie,js:je,ks:ke) + integer:: status, var5id + integer:: start(5), icount(5) + + start(1) = is + start(2) = js + start(3) = ks + start(4) = ntile + start(5) = ntime + + icount(1) = ie - is + 1 + icount(2) = je - js + 1 + icount(3) = ke - ks + 1 + icount(4) = 1 ! one tile at a time + icount(5) = 1 ! one time level at a time + + status = nf_inq_varid (ncid, var5_name, var5id) + + status = nf_get_vara_real(ncid, var5id, start, icount, var5) + + if (status .ne. NF_NOERR) call handle_err('get_var5_r4 get_vara_real '//var5_name,status) + + end subroutine get_var5_r4 + + subroutine get_var5_double( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) + implicit none +#include + integer, intent(in):: ncid + character*(*), intent(in):: var5_name + integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime + real*8, intent(out):: var5(is:ie,js:je,ks:ke) + integer:: status, var5id + integer:: start(5), icount(5) + + start(1) = is + start(2) = js + start(3) = ks + start(4) = ntile + start(5) = ntime + + icount(1) = ie - is + 1 + icount(2) = je - js + 1 + icount(3) = ke - ks + 1 + icount(4) = 1 ! one tile at a time + icount(5) = 1 ! one time level at a time + + status = nf_inq_varid (ncid, var5_name, var5id) + + status = nf_get_vara_double(ncid, var5id, start, icount, var5) + + if (status .ne. NF_NOERR) call handle_err('get_var5_real get_vara_real '//var5_name,status) + + end subroutine get_var5_double !------------------------------------------------------------------------ subroutine get_real3( ncid, var4_name, im, jm, nt, var4 ) From 7262f54daa81b92ca405be2f3bcd75fd25a3f371 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Wed, 15 May 2024 15:16:51 +0000 Subject: [PATCH 05/16] Clean up --- CMakeLists.txt | 2 +- tools/cubed_sphere_inc_mod.F90 | 69 ++++++++++++++++++++++++++++++++++ tools/fv_iau_mod.F90 | 2 +- tools/fv_treat_da_inc.F90 | 10 +++-- 4 files changed, 78 insertions(+), 5 deletions(-) create mode 100644 tools/cubed_sphere_inc_mod.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 478509e6b..64de04747 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -102,7 +102,7 @@ list(APPEND tools_srcs tools/sim_nc_mod.F90 tools/sorted_index.F90 tools/test_cases.F90 - tools/module_cubed_sphere_inc.F90) + tools/cubed_sphere_inc_mod.F90) list(APPEND tools_srcs_extra tools/fv_iau_mod.F90) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 new file mode 100644 index 000000000..c0c2edfe8 --- /dev/null +++ b/tools/cubed_sphere_inc_mod.F90 @@ -0,0 +1,69 @@ +module cubed_sphere_inc_mod + + use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names + use field_manager_mod, only: MODEL_ATMOS + use fv_arrays_mod, only: fv_atmos_type + use sim_nc_mod, only: open_ncfile, get_var5_r4, check_var_exists, close_ncfile + + implicit none + type increment_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 increment_data_type + + public read_cubed_sphere_inc, increment_data_type + + contains + +!---------------------------------------------------------------------------------------- + + subroutine read_cubed_sphere_inc(filename, increment_data, Atm) + character(*), intent(in) :: filename + type(increment_data_type), intent(inout) :: increment_data + type(fv_atmos_type), intent(in) :: Atm + + integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile + character(len=64) :: tracer_name + integer :: ncid, ncerr + + ! Get various dimensions + call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + isc = Atm%bd%isc + iec = Atm%bd%iec + jsc = Atm%bd%jsc + jec = Atm%bd%jec + npz = Atm%npz + itile = Atm%tile_of_mosaic + + ! Open increment file + call open_ncfile( trim(filename), ncid ) + + ! Read increments + call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc ) + call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc ) + call get_var5_r4( ncid, 'temp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) + call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc ) + if ( .not. Atm%flagstruct%hydrostatic ) then + call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc ) + end if + + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + call check_var_exists(ncid, trim(tracer_name)//'_inc', ncerr) + if ( ncerr == 0 ) then + call get_var5_r4( ncid, trim(tracer_name)//'_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%tracer_inc(:,:,:,itracer) ) + end if + end do + + ! Close increment file + call close_ncfile(ncid) + + end subroutine read_cubed_sphere_inc + +!---------------------------------------------------------------------------------------- +end module cubed_sphere_inc_mod diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index a91c91d89..19375a373 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -66,7 +66,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_cubed_sphere_inc, only : read_cubed_sphere_inc, increment_data_type + use cubed_sphere_inc_mod, only : read_cubed_sphere_inc, increment_data_type implicit none private diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index e7bd2795b..5e97fcec9 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -103,6 +103,10 @@ module fv_treat_da_inc_mod ! tracer_manager_mod ! get_tracer_names, get_number_tracers, get_tracer_index ! +! +! cubed_sphere_inc_mod +! read_cubed_sphere_inc, increment_data_type +! ! use fms2_io_mod, only: file_exists @@ -145,8 +149,8 @@ module fv_treat_da_inc_mod get_var3_r4, & get_var1_real, & check_var_exists - use module_cubed_sphere_inc, only: read_cubed_sphere_inc, & - increment_data_type + use cubed_sphere_inc_mod, only: read_cubed_sphere_inc, & + increment_data_type implicit none private @@ -601,7 +605,7 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & 1., ua_inc, va_inc, u, v, & Atm%gridstruct, Atm%npx, Atm%npy, npz_in, fv_domain) endif - + ! Remaining increments ! -------------------- From 3ea6348bd316929c84d51e4234723db095c4b04c Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Fri, 17 May 2024 14:07:40 +0000 Subject: [PATCH 06/16] Final updates (I hope) --- tools/cubed_sphere_inc_mod.F90 | 2 +- tools/fv_restart.F90 | 2 +- tools/fv_treat_da_inc.F90 | 7 ++++++- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index c0c2edfe8..1fffa9691 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -45,7 +45,7 @@ subroutine read_cubed_sphere_inc(filename, increment_data, Atm) ! Read increments call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc ) call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc ) - call get_var5_r4( ncid, 'temp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) + call get_var5_r4( ncid, 'T_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc ) if ( .not. Atm%flagstruct%hydrostatic ) then call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc ) diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index 6e289b941..9b6165d3e 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -88,7 +88,7 @@ module fv_restart_mod ! ! ! fv_treat_da_inc_mod -! read_da_inc +! read_da_inc, read_da_inc_cubed_sphere ! ! ! fv_timing_mod diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 5e97fcec9..759b527ca 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -159,7 +159,7 @@ module fv_treat_da_inc_mod contains !============================================================================= !>@brief The subroutine 'read_da_inc' reads the increments of the diagnostic variables - !! from the DA-generated files. + !! from the DA-generated files on a lot-lon grid. !>@details Additional support of prognostic variables such as tracers can be assessed !! and added upon request. !>@author Xi.Chen @@ -472,6 +472,11 @@ end subroutine apply_inc_on_3d_scalar !--------------------------------------------------------------------------- end subroutine read_da_inc + !============================================================================= + !>@brief The subroutine 'read_da_inc_cubed_sphere' reads increments of diagnostic variables + !! from DA-generated files on the native cubed-sphere grid. + !>@author David.New + !>@date 05/16/2024 subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & u, v, q, delp, pt, delz, & is_in, js_in, ie_in, je_in, isc_in, jsc_in, iec_in, jec_in) From 3891c65e55823d53e5d90e2a2b1b9043d1175ebc Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 21 May 2024 20:35:28 +0000 Subject: [PATCH 07/16] Spelling fix --- tools/fv_treat_da_inc.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 759b527ca..3cf1f3e69 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -159,7 +159,7 @@ module fv_treat_da_inc_mod contains !============================================================================= !>@brief The subroutine 'read_da_inc' reads the increments of the diagnostic variables - !! from the DA-generated files on a lot-lon grid. + !! from the DA-generated files on a lat-lon grid. !>@details Additional support of prognostic variables such as tracers can be assessed !! and added upon request. !>@author Xi.Chen From 866ac571f49bcda9baa68bff4a003e9584254209 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Thu, 30 May 2024 15:02:22 +0000 Subject: [PATCH 08/16] Deleted accidentally undeleted files --- tools/module_cubed_sphere_inc.F90 | 69 ---------------------- tools/module_get_cubed_sphere_inc.F90 | 82 --------------------------- 2 files changed, 151 deletions(-) delete mode 100644 tools/module_cubed_sphere_inc.F90 delete mode 100644 tools/module_get_cubed_sphere_inc.F90 diff --git a/tools/module_cubed_sphere_inc.F90 b/tools/module_cubed_sphere_inc.F90 deleted file mode 100644 index f44205d0c..000000000 --- a/tools/module_cubed_sphere_inc.F90 +++ /dev/null @@ -1,69 +0,0 @@ -module module_cubed_sphere_inc - - use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names - use field_manager_mod, only: MODEL_ATMOS - use fv_arrays_mod, only: fv_atmos_type - use sim_nc_mod, only: open_ncfile, get_var5_r4, check_var_exists, close_ncfile - - implicit none - type increment_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 increment_data_type - - public read_cubed_sphere_inc, increment_data_type - - contains - -!---------------------------------------------------------------------------------------- - - subroutine read_cubed_sphere_inc(filename, increment_data, Atm) - character(*), intent(in) :: filename - type(increment_data_type), intent(inout) :: increment_data - type(fv_atmos_type), intent(in) :: Atm - - integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile - character(len=64) :: tracer_name - integer :: ncid, ncerr - - ! Get various dimensions - call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - isc = Atm%bd%isc - iec = Atm%bd%iec - jsc = Atm%bd%jsc - jec = Atm%bd%jec - npz = Atm%npz - itile = Atm%tile_of_mosaic - - ! Open increment file - call open_ncfile( trim(filename), ncid ) - - ! Read increments - call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc ) - call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc ) - call get_var5_r4( ncid, 'temp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) - call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc ) - if ( .not. Atm%flagstruct%hydrostatic ) then - call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc ) - end if - - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - call check_var_exists(ncid, trim(tracer_name)//'_inc', ncerr) - if ( ncerr == 0 ) then - call get_var5_r4( ncid, trim(tracer_name)//'_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%tracer_inc(:,:,:,itracer) ) - end if - end do - - ! Close increment file - call close_ncfile(ncid) - - end subroutine read_cubed_sphere_inc - -!---------------------------------------------------------------------------------------- -end module module_cubed_sphere_inc diff --git a/tools/module_get_cubed_sphere_inc.F90 b/tools/module_get_cubed_sphere_inc.F90 deleted file mode 100644 index 5d3e6b1e4..000000000 --- a/tools/module_get_cubed_sphere_inc.F90 +++ /dev/null @@ -1,82 +0,0 @@ -module module_get_cubed_sphere_inc - - use netcdf - use esmf - use mpp_mod, only: mpp_pe, mpp_get_current_pelist - use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names - use field_manager_mod, only: MODEL_ATMOS - use CCPP_data, only: GFS_control - use time_manager_mod, only: time_type, get_time, set_time - use fv_arrays_mod, only: fv_atmos_type - use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfFile_t - use mpi - - implicit none - 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 - - public read_netcdf_inc, iau_internal_data_type - - contains - -!---------------------------------------------------------------------------------------- - - subroutine read_netcdf_inc(filename, increment_data, Atm) - character(*), intent(in) :: filename - type(iau_internal_data_type), intent(inout) :: increment_data - type(fv_atmos_type), intent(in) :: Atm - - type(FmsNetcdfFile_t) :: fileobj - integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile - integer, dimension(4) :: corner, edge_lengths - character(len=64) :: tracer_name - - ! Get various dimensions - call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - isc = Atm%bd%isc - iec = Atm%bd%iec - jsc = Atm%bd%jsc - jec = Atm%bd%jec - npz = Atm%npz - itile = Atm%tile_of_mosaic - - ! Open increment file - if ( open_file(fileobj, trim(filename), 'read', pelist=Atm%pelist) ) then - - ! - corner = (/ isc, jsc, 1, itile /) - edge_lengths = (/ iec-isc+1, jec-jsc+1, npz, 1 /) - - ! Read increments - call read_data(fileobj, 'u_inc', increment_data%ua_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'v_inc', increment_data%va_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'T_inc', increment_data%temp_inc, corner=corner, edge_lengths=edge_lengths) - call read_data(fileobj, 'delp_inc', increment_data%delp_inc, corner=corner, edge_lengths=edge_lengths) - if ( .not. Atm%flagstruct%hydrostatic ) then - call read_data(fileobj, 'delz_inc', increment_data%delz_inc, corner=corner, edge_lengths=edge_lengths) - endif - - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then - call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer), & - corner=corner, edge_lengths=edge_lengths) - end if - end do - - ! Close increment file - call close_file(fileobj) - - end if - - end subroutine read_netcdf_inc - -!---------------------------------------------------------------------------------------- -end module module_get_cubed_sphere_inc From c9be39460689673f6dfd578d4561355e5affae35 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Thu, 30 May 2024 16:44:02 +0000 Subject: [PATCH 09/16] Fix bug in GCC compilation related to float conversion --- tools/sim_nc_mod.F90 | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/tools/sim_nc_mod.F90 b/tools/sim_nc_mod.F90 index 19ce370c4..0ffd0f10a 100644 --- a/tools/sim_nc_mod.F90 +++ b/tools/sim_nc_mod.F90 @@ -49,6 +49,11 @@ module sim_nc_mod handle_err, check_var, get_var1_real, get_var_att_double, & check_var_exists, get_var5_r4 + interface get_var5_r4 + module procedure get_var5_r4_to_r4 + module procedure get_var5_r4_to_double + end interface get_var5_r4 + contains subroutine open_ncfile( iflnm, ncid ) @@ -340,7 +345,7 @@ subroutine get_var4_double( ncid, var4_name, im, jm, km, nt, var4 ) end subroutine get_var4_double - subroutine get_var5_r4( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) + subroutine get_var5_r4_to_r4( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) implicit none #include integer, intent(in):: ncid @@ -368,7 +373,37 @@ subroutine get_var5_r4( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 if (status .ne. NF_NOERR) call handle_err('get_var5_r4 get_vara_real '//var5_name,status) - end subroutine get_var5_r4 + end subroutine get_var5_r4_to_r4 + + subroutine get_var5_r4_to_double( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) + implicit none +#include + integer, intent(in):: ncid + character*(*), intent(in):: var5_name + integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime + real*8, intent(out):: var5(is:ie,js:je,ks:ke) + integer:: status, var5id + integer:: start(5), icount(5) + + start(1) = is + start(2) = js + start(3) = ks + start(4) = ntile + start(5) = ntime + + icount(1) = ie - is + 1 + icount(2) = je - js + 1 + icount(3) = ke - ks + 1 + icount(4) = 1 ! one tile at a time + icount(5) = 1 ! one time level at a time + + status = nf_inq_varid (ncid, var5_name, var5id) + + status = nf_get_vara_real(ncid, var5id, start, icount, var5) + + if (status .ne. NF_NOERR) call handle_err('get_var5_r4 get_vara_real '//var5_name,status) + + end subroutine get_var5_r4_to_double subroutine get_var5_double( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) implicit none From ff741fd7afabf9b9e9d0079224e0c3de6853864b Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 2 Jul 2024 09:28:23 -0500 Subject: [PATCH 10/16] Initial commit --- tools/cubed_sphere_inc_mod.F90 | 61 ++++++++++++++++------------------ 1 file changed, 28 insertions(+), 33 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 1fffa9691..4657a69df 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -3,8 +3,8 @@ module cubed_sphere_inc_mod use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names use field_manager_mod, only: MODEL_ATMOS use fv_arrays_mod, only: fv_atmos_type - use sim_nc_mod, only: open_ncfile, get_var5_r4, check_var_exists, close_ncfile - + use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfFile_t + implicit none type increment_data_type real, allocatable :: ua_inc(:,:,:) @@ -21,47 +21,42 @@ module cubed_sphere_inc_mod !---------------------------------------------------------------------------------------- - subroutine read_cubed_sphere_inc(filename, increment_data, Atm) - character(*), intent(in) :: filename + subroutine read_cubed_sphere_inc(fname_prefix, increment_data, Atm) + character(*), intent(in) :: fname_prefix type(increment_data_type), intent(inout) :: increment_data type(fv_atmos_type), intent(in) :: Atm - integer :: isc, iec, jsc, jec, npz, itracer, ntracers, itile - character(len=64) :: tracer_name - integer :: ncid, ncerr + type(FmsNetcdfFile_t) :: fileobj + integer :: itracer, ntracers, itile + character(len=64) :: tracer_name + character(len=1) :: itile_str ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - isc = Atm%bd%isc - iec = Atm%bd%iec - jsc = Atm%bd%jsc - jec = Atm%bd%jec - npz = Atm%npz itile = Atm%tile_of_mosaic + write(itile_str, '(I0)') itile - ! Open increment file - call open_ncfile( trim(filename), ncid ) - - ! Read increments - call get_var5_r4( ncid, 'u_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%ua_inc ) - call get_var5_r4( ncid, 'v_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%va_inc ) - call get_var5_r4( ncid, 'T_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%temp_inc ) - call get_var5_r4( ncid, 'delp_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delp_inc ) - if ( .not. Atm%flagstruct%hydrostatic ) then - call get_var5_r4( ncid, 'delz_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%delz_inc ) - end if - - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - call check_var_exists(ncid, trim(tracer_name)//'_inc', ncerr) - if ( ncerr == 0 ) then - call get_var5_r4( ncid, trim(tracer_name)//'_inc', isc,iec, jsc,jec, 1,npz, itile, 1, increment_data%tracer_inc(:,:,:,itracer) ) + ! Open file + if (open_file(fileobj, trim(fname_prefix) // '.tile' // itile_str // '.nc', "read", Atm%domain)) then + ! Read increments + call read_data(fileobj, 'u_inc', increment_data%ua_inc) + call read_data(fileobj, 'v_inc', increment_data%va_inc) + call read_data(fileobj, 'T_inc', increment_data%T_inc) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc) + if ( .not. Atm%flagstruct%hydrostatic ) then + call read_data(fileobj, 'delz_inc', increment_data%delz_inc) end if - end do + + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then + call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) + end if + end do - ! Close increment file - call close_ncfile(ncid) + call close_file(fileobj) + end if end subroutine read_cubed_sphere_inc From 2335f8ef39ab6ae262fe03ea397892698cecc6c9 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 2 Jul 2024 09:30:04 -0500 Subject: [PATCH 11/16] revert sim_nc_mod --- tools/sim_nc_mod.F90 | 97 +------------------------------------------- 1 file changed, 1 insertion(+), 96 deletions(-) diff --git a/tools/sim_nc_mod.F90 b/tools/sim_nc_mod.F90 index 0ffd0f10a..b461aba45 100644 --- a/tools/sim_nc_mod.F90 +++ b/tools/sim_nc_mod.F90 @@ -47,13 +47,8 @@ module sim_nc_mod public open_ncfile, close_ncfile, get_ncdim1, get_var1_double, get_var2_double, & get_var3_real, get_var3_double, get_var3_r4, get_var2_real, get_var2_r4, & handle_err, check_var, get_var1_real, get_var_att_double, & - check_var_exists, get_var5_r4 + check_var_exists - interface get_var5_r4 - module procedure get_var5_r4_to_r4 - module procedure get_var5_r4_to_double - end interface get_var5_r4 - contains subroutine open_ncfile( iflnm, ncid ) @@ -344,96 +339,6 @@ subroutine get_var4_double( ncid, var4_name, im, jm, km, nt, var4 ) if (status .ne. NF_NOERR) call handle_err('get_var4_double get_vara_double '//var4_name,status) end subroutine get_var4_double - - subroutine get_var5_r4_to_r4( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) - implicit none -#include - integer, intent(in):: ncid - character*(*), intent(in):: var5_name - integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime - real*4, intent(out):: var5(is:ie,js:je,ks:ke) - integer:: status, var5id - integer:: start(5), icount(5) - - start(1) = is - start(2) = js - start(3) = ks - start(4) = ntile - start(5) = ntime - - icount(1) = ie - is + 1 - icount(2) = je - js + 1 - icount(3) = ke - ks + 1 - icount(4) = 1 ! one tile at a time - icount(5) = 1 ! one time level at a time - - status = nf_inq_varid (ncid, var5_name, var5id) - - status = nf_get_vara_real(ncid, var5id, start, icount, var5) - - if (status .ne. NF_NOERR) call handle_err('get_var5_r4 get_vara_real '//var5_name,status) - - end subroutine get_var5_r4_to_r4 - - subroutine get_var5_r4_to_double( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) - implicit none -#include - integer, intent(in):: ncid - character*(*), intent(in):: var5_name - integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime - real*8, intent(out):: var5(is:ie,js:je,ks:ke) - integer:: status, var5id - integer:: start(5), icount(5) - - start(1) = is - start(2) = js - start(3) = ks - start(4) = ntile - start(5) = ntime - - icount(1) = ie - is + 1 - icount(2) = je - js + 1 - icount(3) = ke - ks + 1 - icount(4) = 1 ! one tile at a time - icount(5) = 1 ! one time level at a time - - status = nf_inq_varid (ncid, var5_name, var5id) - - status = nf_get_vara_real(ncid, var5id, start, icount, var5) - - if (status .ne. NF_NOERR) call handle_err('get_var5_r4 get_vara_real '//var5_name,status) - - end subroutine get_var5_r4_to_double - - subroutine get_var5_double( ncid, var5_name, is,ie, js,je, ks,ke, ntile, ntime, var5 ) - implicit none -#include - integer, intent(in):: ncid - character*(*), intent(in):: var5_name - integer, intent(in):: is, ie, js, je, ks, ke, ntile, ntime - real*8, intent(out):: var5(is:ie,js:je,ks:ke) - integer:: status, var5id - integer:: start(5), icount(5) - - start(1) = is - start(2) = js - start(3) = ks - start(4) = ntile - start(5) = ntime - - icount(1) = ie - is + 1 - icount(2) = je - js + 1 - icount(3) = ke - ks + 1 - icount(4) = 1 ! one tile at a time - icount(5) = 1 ! one time level at a time - - status = nf_inq_varid (ncid, var5_name, var5id) - - status = nf_get_vara_double(ncid, var5id, start, icount, var5) - - if (status .ne. NF_NOERR) call handle_err('get_var5_real get_vara_real '//var5_name,status) - - end subroutine get_var5_double !------------------------------------------------------------------------ subroutine get_real3( ncid, var4_name, im, jm, nt, var4 ) From 39508026528eea2d8df01cf9238fa2ab9a4a40dd Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Wed, 3 Jul 2024 21:39:38 +0000 Subject: [PATCH 12/16] Implement FMS increment reads --- tools/cubed_sphere_inc_mod.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 4657a69df..433e16622 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -3,7 +3,7 @@ module cubed_sphere_inc_mod use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names use field_manager_mod, only: MODEL_ATMOS use fv_arrays_mod, only: fv_atmos_type - use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfFile_t + use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfDomainFile_t implicit none type increment_data_type @@ -26,10 +26,10 @@ subroutine read_cubed_sphere_inc(fname_prefix, increment_data, Atm) type(increment_data_type), intent(inout) :: increment_data type(fv_atmos_type), intent(in) :: Atm - type(FmsNetcdfFile_t) :: fileobj - integer :: itracer, ntracers, itile - character(len=64) :: tracer_name - character(len=1) :: itile_str + type(FmsNetcdfDomainFile_t) :: fileobj + integer :: itracer, ntracers, itile + character(len=64) :: tracer_name + character(len=1) :: itile_str ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) @@ -37,11 +37,11 @@ subroutine read_cubed_sphere_inc(fname_prefix, increment_data, Atm) write(itile_str, '(I0)') itile ! Open file - if (open_file(fileobj, trim(fname_prefix) // '.tile' // itile_str // '.nc', "read", Atm%domain)) then + if ( open_file(fileobj, trim(fname_prefix) // '.tile' // itile_str // '.nc', "read", Atm%domain) ) then ! Read increments call read_data(fileobj, 'u_inc', increment_data%ua_inc) call read_data(fileobj, 'v_inc', increment_data%va_inc) - call read_data(fileobj, 'T_inc', increment_data%T_inc) + call read_data(fileobj, 'T_inc', increment_data%temp_inc) call read_data(fileobj, 'delp_inc', increment_data%delp_inc) if ( .not. Atm%flagstruct%hydrostatic ) then call read_data(fileobj, 'delz_inc', increment_data%delz_inc) From 4d6e8d13402c9525aa862a1667c1ebbc123ecb77 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Thu, 18 Jul 2024 12:58:56 +0000 Subject: [PATCH 13/16] Fix filename issues --- tools/cubed_sphere_inc_mod.F90 | 11 ++++------- tools/fv_treat_da_inc.F90 | 8 ++++++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 433e16622..bc6382542 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -21,23 +21,20 @@ module cubed_sphere_inc_mod !---------------------------------------------------------------------------------------- - subroutine read_cubed_sphere_inc(fname_prefix, increment_data, Atm) - character(*), intent(in) :: fname_prefix + subroutine read_cubed_sphere_inc(fname, increment_data, Atm) + character(*), intent(in) :: fname type(increment_data_type), intent(inout) :: increment_data type(fv_atmos_type), intent(in) :: Atm type(FmsNetcdfDomainFile_t) :: fileobj - integer :: itracer, ntracers, itile + integer :: itracer, ntracers character(len=64) :: tracer_name - character(len=1) :: itile_str ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - itile = Atm%tile_of_mosaic - write(itile_str, '(I0)') itile ! Open file - if ( open_file(fileobj, trim(fname_prefix) // '.tile' // itile_str // '.nc', "read", Atm%domain) ) then + if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then ! Read increments call read_data(fileobj, 'u_inc', increment_data%ua_inc) call read_data(fileobj, 'v_inc', increment_data%va_inc) diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 3cf1f3e69..7e985fb13 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -491,17 +491,20 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & real, intent(inout) :: q( is_in:ie_in, js_in:je_in, npz_in, nq) ! tracers real, intent(inout) :: delz(isc_in:iec_in, jsc_in:jec_in, npz_in) ! layer thickness - character(len=128) :: fname + character(len=128) :: fname_prefix, fname integer :: sphum, liq_wat, spo, spo2, spo3, o3mr, i, j, k, itracer type(increment_data_type) :: increment_data type(group_halo_update_type) :: i_pack(2) real :: ua_inc(is_in:ie_in, js_in:je_in, npz_in) real :: va_inc(is_in:ie_in, js_in:je_in, npz_in) + character(len=1) :: itile_str ! Get increment filename - fname = 'INPUT/'//Atm%flagstruct%res_latlon_dynamics + fname_prefix = 'INPUT/'//Atm%flagstruct%res_latlon_dynamics ! Ensure file exists + write(itile_str, '(I0)') Atm%tile_of_mosaic + fname = trim(fname_prefix) // '.tile' // itile_str // '.nc' if ( .not. file_exists(fname) ) then call mpp_error(FATAL,'==> Error in read_da_inc_cubed_sphere: Expected file '& //trim(fname)//' for DA increment does not exist') @@ -518,6 +521,7 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & allocate(increment_data%tracer_inc(isc_in:iec_in,jsc_in:jec_in,npz_in,nq)) ! Read increments + fname = trim(fname_prefix) // '.nc' call read_cubed_sphere_inc(fname, increment_data, Atm) ! Wind increments From 409eca8ea893fa6fbc5937c612a30df38d8f2437 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Tue, 30 Jul 2024 19:04:52 +0000 Subject: [PATCH 14/16] Make sure to register axes --- tools/cubed_sphere_inc_mod.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index bc6382542..03eb467c2 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -3,7 +3,8 @@ module cubed_sphere_inc_mod use tracer_manager_mod, only: get_tracer_index, get_number_tracers, get_tracer_names use field_manager_mod, only: MODEL_ATMOS use fv_arrays_mod, only: fv_atmos_type - use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, FmsNetcdfDomainFile_t + use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, & + FmsNetcdfDomainFile_t, register_axis implicit none type increment_data_type @@ -35,6 +36,10 @@ subroutine read_cubed_sphere_inc(fname, increment_data, Atm) ! Open file if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then + ! Register axes + call register_axis(fileobj, 'grid_xt', 'x') + call register_axis(fileobj, 'grid_yt', 'y') + ! Read increments call read_data(fileobj, 'u_inc', increment_data%ua_inc) call read_data(fileobj, 'v_inc', increment_data%va_inc) From b7c4d907cb3fe074e604573b76aa3eef8411c9a8 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA <134300700+DavidNew-NOAA@users.noreply.github.com> Date: Mon, 5 Aug 2024 17:02:18 -0400 Subject: [PATCH 15/16] Update x and y-axis names for increment to match restart axis names --- tools/cubed_sphere_inc_mod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 03eb467c2..d64d7077a 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -37,8 +37,8 @@ subroutine read_cubed_sphere_inc(fname, increment_data, Atm) ! Open file if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then ! Register axes - call register_axis(fileobj, 'grid_xt', 'x') - call register_axis(fileobj, 'grid_yt', 'y') + call register_axis(fileobj, 'xaxis_1', 'x') + call register_axis(fileobj, 'yaxis_2', 'y') ! Read increments call read_data(fileobj, 'u_inc', increment_data%ua_inc) From 04b336f3da7fde10f1c63ca82699fd5c1378fbf3 Mon Sep 17 00:00:00 2001 From: DavidNew-NOAA Date: Thu, 15 Aug 2024 13:42:08 +0000 Subject: [PATCH 16/16] Last minute update of y-axis name --- tools/cubed_sphere_inc_mod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index d64d7077a..2b169b6e7 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -38,7 +38,7 @@ subroutine read_cubed_sphere_inc(fname, increment_data, Atm) if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then ! Register axes call register_axis(fileobj, 'xaxis_1', 'x') - call register_axis(fileobj, 'yaxis_2', 'y') + call register_axis(fileobj, 'yaxis_1', 'y') ! Read increments call read_data(fileobj, 'u_inc', increment_data%ua_inc)