diff --git a/CMakeLists.txt b/CMakeLists.txt index def27a8e6..0d588415f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -102,7 +102,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/cubed_sphere_inc_mod.F90) list(APPEND tools_srcs_extra tools/fv_iau_mod.F90) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index cddf08cd6..d05167078 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 :: 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 ba1f88e9e..2ead0274c 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 :: increment_file_on_native_grid logical , pointer :: hydrostatic logical , pointer :: phys_hydrostatic logical , pointer :: use_hydro_pressure @@ -950,7 +951,7 @@ subroutine set_namelist_pointers(Atm) external_ic => Atm%flagstruct%external_ic external_eta => Atm%flagstruct%external_eta read_increment => Atm%flagstruct%read_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,10 +1094,9 @@ 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, & + write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, increment_file_on_native_grid, & pass_full_omega_to_physics_in_non_hydrostatic_mode, ignore_rst_cksum - ! Read FVCORE namelist read (input_nml_file,fv_core_nml,iostat=ios) ierr = check_nml_error(ios,'fv_core_nml') diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 new file mode 100644 index 000000000..2b169b6e7 --- /dev/null +++ b/tools/cubed_sphere_inc_mod.F90 @@ -0,0 +1,66 @@ +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, register_axis + + 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(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 + character(len=64) :: tracer_name + + ! Get various dimensions + call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + + ! Open file + if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then + ! Register axes + call register_axis(fileobj, 'xaxis_1', 'x') + call register_axis(fileobj, 'yaxis_1', 'y') + + ! 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%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) + end if + + ! 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 + + call close_file(fileobj) + end if + + 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 e4ca56500..19375a373 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 cubed_sphere_inc_mod, only : read_cubed_sphere_inc, increment_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(:,:,:) @@ -103,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 @@ -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%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_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) 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%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_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2))) + 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%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_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) endif @@ -434,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 acc765686..dae3d0496 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 @@ -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 @@ -471,11 +471,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) - 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) + 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====== endif !n==this_grid diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 4fbee842d..7e985fb13 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 @@ -126,12 +130,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,15 +149,17 @@ module fv_treat_da_inc_mod get_var3_r4, & get_var1_real, & check_var_exists + use cubed_sphere_inc_mod, only: read_cubed_sphere_inc, & + increment_data_type implicit none private - public :: read_da_inc,remap_coef + public :: read_da_inc, read_da_inc_cubed_sphere, remap_coef 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 lat-lon grid. !>@details Additional support of prognostic variables such as tracers can be assessed !! and added upon request. !>@author Xi.Chen @@ -461,6 +472,187 @@ 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) + 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) ! tracers + real, intent(inout) :: delz(isc_in:iec_in, jsc_in:jec_in, npz_in) ! layer thickness + + 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_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') + 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 + fname = trim(fname_prefix) // '.nc' + call read_cubed_sphere_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), 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), 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(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 + ! Edges + 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 ( 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 ( 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 ( 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 ( 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(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 + + ! 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 + 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(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(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 + !============================================================================= !>@brief The subroutine 'remap_coef' calculates the coefficients for horizonal regridding.