diff --git a/CMakeLists.txt b/CMakeLists.txt index 37b0066fa..415a0068f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,7 +54,8 @@ endif() if(NOT FMS_FOUND) find_package(FMS REQUIRED COMPONENTS ${kind}) - add_library(fms ALIAS FMS::fms_${kind}) + string(TOLOWER ${kind} kind_lower) + add_library(fms ALIAS FMS::fms_${kind_lower}) endif() list(APPEND moving_srcs @@ -189,7 +190,7 @@ set_property(SOURCE model/fv_mapz.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "${F set_property(SOURCE tools/fv_diagnostics.F90 APPEND_STRING PROPERTY COMPILE_FLAGS "-O0") set_target_properties(fv3 PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include/fv3) -target_include_directories(fv3 INTERFACE $ +target_include_directories(fv3 PUBLIC $ $) target_include_directories(fv3 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 7af45f664..d4ba3e672 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -265,7 +265,7 @@ module atmosphere_mod logical :: cold_start = .false. ! used in initial condition integer, dimension(:), allocatable :: id_tracerdt_dyn - integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, cld_amt ! condensate species tracer indices + integer :: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, cld_amt ! condensate species tracer indices integer :: mygrid = 1 integer :: p_split = 1 @@ -324,10 +324,12 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) logical :: dycore_only = .false. logical :: debug = .false. logical :: sync = .false. + logical :: ignore_rst_cksum = .false. integer, parameter :: maxhr = 4096 real, dimension(maxhr) :: fdiag = 0. real :: fhmax=384.0, fhmaxhf=120.0, fhout=3.0, fhouthf=1.0,avg_max_length=3600. - namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, avg_max_length + namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, & + & avg_max_length, ignore_rst_cksum ! *DH 20210326 !For regional @@ -357,7 +359,7 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) #ifdef MOVING_NEST call fv_tracker_init(size(Atm)) if (mygrid .eq. 2) call allocate_tracker(mygrid, Atm(mygrid)%bd%isc, Atm(mygrid)%bd%iec, Atm(mygrid)%bd%jsc, Atm(mygrid)%bd%jec) -#endif +#endif Atm(mygrid)%Time_init = Time_init @@ -403,9 +405,10 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat' ) snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat' ) graupel = get_tracer_index (MODEL_ATMOS, 'graupel' ) + hailwat = get_tracer_index (MODEL_ATMOS, 'hailwat' ) cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt') - if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel) > Atm(mygrid)%flagstruct%nwat) then + if (max(sphum,liq_wat,ice_wat,rainwat,snowwat,graupel,hailwat) > Atm(mygrid)%flagstruct%nwat) then call mpp_error (FATAL,' atmosphere_init: condensate species are not first in the list of & &tracers defined in the field_table') endif @@ -449,6 +452,15 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) !--- allocate pref allocate(pref(npz+1,2), dum1d(npz+1)) + ! DH* 20210326 + ! First, read atmos_model_nml namelist section - this is a workaround to avoid + ! unnecessary additional changes to the input namelists, in anticipation of the + ! implementation of a generic interface for GFDL and CCPP fast physics soon + read(input_nml_file, nml=atmos_model_nml, iostat=io) + ierr = check_nml_error(io, 'atmos_model_nml') + !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'" + ! *DH 20210326 + call fv_restart(Atm(mygrid)%domain, Atm, seconds, days, cold_start, Atm(mygrid)%gridstruct%grid_type, mygrid) fv_time = Time @@ -490,15 +502,6 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) ! Do CCPP fast physics initialization before call to adiabatic_init (since this calls fv_dynamics) - ! DH* 20210326 - ! First, read atmos_model_nml namelist section - this is a workaround to avoid - ! unnecessary additional changes to the input namelists, in anticipation of the - ! implementation of a generic interface for GFDL and CCPP fast physics soon - read(input_nml_file, nml=atmos_model_nml, iostat=io) - ierr = check_nml_error(io, 'atmos_model_nml') - !write(0,'(a)') "It's me, and my physics suite is '" // trim(ccpp_suite) // "'" - ! *DH 20210326 - ! For fast physics running over the entire domain, block ! and thread number are not used; set to safe values cdata%blk_no = 1 @@ -953,10 +956,10 @@ end subroutine get_nth_domain_info !! the "domain2d" variable associated with the coupling grid and the !! decomposition for the current cubed-sphere tile. !>@detail Coupling is done using the mass/temperature grid with no halos. - subroutine atmosphere_domain ( fv_domain, layout, regional, nested, & + subroutine atmosphere_domain ( fv_domain, rd_domain, layout, regional, nested, & moving_nest_parent, is_moving_nest, & ngrids_atmos, mygrid_atmos, pelist ) - type(domain2d), intent(out) :: fv_domain + type(domain2d), intent(out) :: fv_domain, rd_domain integer, intent(out) :: layout(2) logical, intent(out) :: regional logical, intent(out) :: nested @@ -969,6 +972,7 @@ subroutine atmosphere_domain ( fv_domain, layout, regional, nested, & integer :: n fv_domain = Atm(mygrid)%domain_for_coupler + rd_domain = Atm(mygrid)%domain_for_read layout(1:2) = Atm(mygrid)%layout(1:2) regional = Atm(mygrid)%flagstruct%regional nested = ngrids > 1 diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 7806df52e..f858d8804 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -1248,7 +1248,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, call timing_on('COMM_TOTAL') #ifndef ROT3 - if ( it/=n_split) & + if ( .not. flagstruct%regional .and. it/=n_split) & call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE) #endif call timing_off('COMM_TOTAL') @@ -1351,7 +1351,10 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, isd, ied, jsd, jed, & reg_bc_update_time,it ) - call mpp_update_domains(u, v, domain, gridtype=DGRID_NE) +#ifndef ROT3 + if (it/=n_split) & + call start_group_halo_update(i_pack(8), u, v, domain, gridtype=DGRID_NE) +#endif endif diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index b89dcefa1..4a0539509 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -516,6 +516,7 @@ module fv_arrays_mod !----------------------------------------------------------------------------------------------- logical :: reset_eta = .false. + logical :: ignore_rst_cksum = .false. !< enfore (.false.) or override (.true.) data integrity restart checksums real :: p_fac = 0.05 !< Safety factor for minimum nonhydrostatic pressures, which !< will be limited so the full pressure is no less than p_fac !< times the hydrostatic pressure. This is only of concern in mid-top @@ -1303,6 +1304,7 @@ module fv_arrays_mod #if defined(SPMD) type(domain2D) :: domain_for_coupler !< domain used in coupled model with halo = 1. + type(domain2D) :: domain_for_read !< domain used for reads to increase performance when io_layout=(1,1) !global tile and tile_of_mosaic only have a meaning for the CURRENT pe integer :: num_contact, npes_per_tile, global_tile, tile_of_mosaic, npes_this_grid diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 5b24ccf56..6951a15a5 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -198,7 +198,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, real, intent(in) :: dt_atmos integer, intent(OUT) :: this_grid logical, allocatable, intent(OUT) :: grids_on_this_pe(:) - character(len=32), optional, intent(in) :: nml_filename_in ! alternate nml + character(len=32), optional, intent(in) :: nml_filename_in ! alternate nml logical, optional, intent(in) :: skip_nml_read_in ! use previously loaded nml integer, intent(INOUT) :: p_split @@ -291,6 +291,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, real(kind=R_GRID) , pointer :: target_lon logical , pointer :: reset_eta + logical , pointer :: ignore_rst_cksum real , pointer :: p_fac real , pointer :: a_imp integer , pointer :: n_split @@ -676,7 +677,8 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, Atm(this_grid)%flagstruct%grid_type,Atm(this_grid)%neststruct%nested, & Atm(this_grid)%layout,Atm(this_grid)%io_layout,Atm(this_grid)%bd,Atm(this_grid)%tile_of_mosaic, & Atm(this_grid)%gridstruct%square_domain,Atm(this_grid)%npes_per_tile,Atm(this_grid)%domain, & - Atm(this_grid)%domain_for_coupler,Atm(this_grid)%num_contact,Atm(this_grid)%pelist) + Atm(this_grid)%domain_for_coupler,Atm(this_grid)%domain_for_read,Atm(this_grid)%num_contact, & + Atm(this_grid)%pelist) call broadcast_domains(Atm,Atm(this_grid)%pelist,size(Atm(this_grid)%pelist)) do n=1,ngrids tile_id = mpp_get_tile_id(Atm(n)%domain) @@ -857,6 +859,7 @@ subroutine set_namelist_pointers(Atm) regional_bcs_from_gsi => Atm%flagstruct%regional_bcs_from_gsi write_restart_with_bcs => Atm%flagstruct%write_restart_with_bcs reset_eta => Atm%flagstruct%reset_eta + ignore_rst_cksum => Atm%flagstruct%ignore_rst_cksum p_fac => Atm%flagstruct%p_fac a_imp => Atm%flagstruct%a_imp n_split => Atm%flagstruct%n_split @@ -1075,7 +1078,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 + write_coarse_agrid_vel_rst, write_coarse_dgrid_vel_rst, ignore_rst_cksum ! Read FVCORE namelist diff --git a/model/fv_fill.F90 b/model/fv_fill.F90 index faab16fdf..5579c7221 100644 --- a/model/fv_fill.F90 +++ b/model/fv_fill.F90 @@ -37,7 +37,7 @@ module fv_fill_mod ! use mpp_domains_mod, only: mpp_update_domains, domain2D - use platform_mod, only: kind_phys => r8_kind + use GFS_typedefs, only: kind_phys implicit none public fillz diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 06a5158b4..4443d04ca 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -235,12 +235,42 @@ module fv_regional_mod character(len=100) :: grid_data='grid.tile7.halo4.nc' & ,oro_data ='oro_data.tile7.halo4.nc' +! +!----------------------------------------------------------------------- +!*** All arrays and variables that contain boundary data should be +!*** initialized to real_snan or dbl_snan (as appropriate) because that +!*** is the only way to detect errors in the boundary generation. +!*** Usual approaches like detecting out-of-bounds accesses won't work +!*** because the arrays may be larger than the boundaries. +!*** +!*** Regions of the arrays that should have valid boundary data should +!*** never contain a NaN. If they do, then there is an error that has +!*** not yet been identified. In particular, the blending code has an +!*** unidentified problem. +!*** +!*** Presently, the blending code is asked to blend data in regions +!*** that contain NaN. This means either: +!*** +!*** 1) not all boundary areas are filled in, OR +!*** 2) the blending code is accessing beyond the boundaries +!*** +!*** The reasons for this are unknown, but the workaround does work: +!*** detect the NaNs and don't use those points in blending. +!*** +!*** Also, the wind remapping code needs surface pressure beyond the +!*** wind boundary regions, but the surface pressure boundary data +!*** doesn't extend that far. +!----------------------------------------------------------------------- +! #ifdef OVERLOAD_R4 - real, parameter:: real_snan=real(Z'FFBFFFFF') + real, parameter:: real_snan=real(Z'FFBFFFFF') ! <-- One of many ways to express 32-bit signalling NaN + real, parameter:: real_qnan=real(Z'FFFFFFFF') ! <-- One of many ways to express 32-bit quiet NaN #else - real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') + real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') ! <-- One of many ways to express 64-bit signalling NaN + real, parameter:: real_qnan=real(Z'FFFFFFFFFFFFFFFF') ! <-- One of many ways to express 64-bit quiet NaN #endif - real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) + real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit signalling NaN + real(kind=R_GRID), parameter:: dbl_qnan=real(Z'FFFFFFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit quiet NaN interface dump_field module procedure dump_field_3d @@ -254,6 +284,68 @@ module fv_regional_mod logical :: data_source_fv3gfs contains +!----------------------------------------------------------------------- +! + logical function is_not_finite(val) +! +!----------------------------------------------------------------------- +!*** This routine is equivalent to ".not. ieee_is_finite(val)" +!*** which returns .true. for infinite and Not a Number (NaN), or +!*** .false. otherwise. It's here as a workaround for this gfortran bug: +!*** +!*** https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82207 +!*** +!*** The compiler must use IEEE-standard floating point for this to work +!----------------------------------------------------------------------- +! + use, intrinsic :: iso_c_binding, only: c_int32_t, c_int64_t + implicit none +! +!----------------------------------------------------------------------- +! Portability note: shiftr() is part of Fortran 2008, but it is widely +! supported in older compilers. +!----------------------------------------------------------------------- +! + intrinsic shiftr, transfer, iand ! <--- for the benefit of older compilers +! +!----------------------------------------------------------------------- +! Use value-based argument passing instead of reference-based to avoid +! signaling a NaN on conversion to addressable storage. +!----------------------------------------------------------------------- +! + real, value :: val ! <-- bit pattern to test for infinity or NaN +! +!----------------------------------------------------------------------- +! Bit manipulation constants for testing 32-bit floating-point +! non-finite values. +!----------------------------------------------------------------------- +! +#ifdef OVERLOAD_R4 + integer(c_int32_t), parameter :: check = 255 ! <-- all bits on, size of exponent (8 bits) + integer, parameter :: shift = 23 ! <-- number of mantissa bits except sign +! +!----------------------------------------------------------------------- +! Bit manipulation constants for testing 64-bit floating-point +! non-finite values. +!----------------------------------------------------------------------- +! +#else + integer(c_int64_t), parameter :: check = 2047 ! <-- all bits on, size of exponent (11 bits) + integer, parameter :: shift = 52 ! <-- number of mantissa bits except sign +#endif +! +!----------------------------------------------------------------------- +! For IEEE standard floating point numbers, non-finite values follow +! a mandatory bit pattern. They have the mantissa sign bit on, and all +! exponent bits on, except the exponent sign which can be on or off. +!----------------------------------------------------------------------- +! + is_not_finite = iand(shiftr(transfer(val,check),shift),check)==check +! + end function is_not_finite +! +!----------------------------------------------------------------------- + !----------------------------------------------------------------------- ! subroutine setup_regional_BC(Atm & @@ -410,6 +502,17 @@ subroutine setup_regional_BC(Atm & else nrows_blend=nrows_blend_in_data !<-- # of blending rows in the BC files. endif + + IF ( north_bc .or. south_bc ) THEN + IF ( nrows_blend_user > jed - nhalo_model - (jsd + nhalo_model) + 1 ) THEN + call mpp_error(FATAL,'Number of blending rows is greater than the north-south tile size!') + ENDIF + ENDIF + IF ( west_bc .or. east_bc ) THEN + IF ( nrows_blend_user > ied - nhalo_model - (isd + nhalo_model) + 1 ) THEN + call mpp_error(FATAL,'Number of blending rows is greater than the east-west tile size!') + ENDIF + ENDIF ! call check(nf90_close(ncid)) !<-- Close the BC file for now. ! @@ -754,8 +857,8 @@ subroutine setup_regional_BC(Atm & !*** reference pressure profile. Compute it now. !----------------------------------------------------------------------- ! - allocate(pref(npz+1)) - allocate(dum1d(npz+1)) + allocate(pref(npz+1)) ; pref=real_snan + allocate(dum1d(npz+1)) ; dum1d=real_snan ! ps1=101325. pref(npz+1)=ps1 @@ -1296,8 +1399,8 @@ subroutine start_regional_cold_start(Atm, ak, bk, levp & enddo enddo ! - allocate (ak_in(1:levp+1)) !<-- Save the input vertical structure for - allocate (bk_in(1:levp+1)) ! remapping BC updates during the forecast. + allocate (ak_in(1:levp+1)) ; ak_in=real_snan !<-- Save the input vertical structure for + allocate (bk_in(1:levp+1)) ; bk_in=real_snan ! remapping BC updates during the forecast. do k=1,levp+1 ak_in(k)=ak(k) bk_in(k)=bk(k) @@ -1391,9 +1494,9 @@ subroutine start_regional_restart(Atm & ,isd, ied, jsd, jed & ,Atm%npx, Atm%npy ) ! - allocate (wk2(levp+1,2)) - allocate (ak_in(levp+1)) !<-- Save the input vertical structure for - allocate (bk_in(levp+1)) ! remapping BC updates during the forecast. + allocate (wk2(levp+1,2)) ; wk2=real_snan + allocate (ak_in(levp+1)) ; ak_in=real_snan !<-- Save the input vertical structure for + allocate (bk_in(levp+1)) ; bk_in=real_snan ! remapping BC updates during the forecast. if (Atm%flagstruct%hrrrv3_ic) then if (open_file(Grid_input, 'INPUT/hrrr_ctrl.nc', "read", pelist=pes)) then call read_data(Grid_input,'vcoord',wk2) @@ -1897,7 +2000,8 @@ subroutine regional_bc_data(Atm,bc_hour & !*** the integration levels. !----------------------------------------------------------------------- ! - allocate(ps_reg(is_input:ie_input,js_input:je_input)) ; ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed + allocate(ps_reg(is_input:ie_input,js_input:je_input)) !<-- Sfc pressure in domain's boundary region derived from BC files + ps_reg=real_snan !<-- detect access of uninitialized pressures ! !----------------------------------------------------------------------- !*** We have the boundary variables from the BC file on the levels @@ -1931,16 +2035,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==2)then - if(south_bc)then - call_remap=.true. - side='south' - bc_side_t1=>BC_t1%south - bc_side_t0=>BC_t0%south - endif - endif -! - if(nside==3)then + if(nside==2)then if(east_bc)then call_remap=.true. side='east ' @@ -1949,7 +2044,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_remap=.true. side='west ' @@ -1957,6 +2052,15 @@ subroutine regional_bc_data(Atm,bc_hour & bc_side_t0=>BC_t0%west endif endif +! + if(nside==4)then + if(south_bc)then + call_remap=.true. + side='south' + bc_side_t1=>BC_t1%south + bc_side_t0=>BC_t0%south + endif + endif ! if(call_remap)then call remap_scalar_nggps_regional_bc(Atm & @@ -2034,57 +2138,9 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==2)then - if(south_bc)then - if (is == 1) then - istart = 1 - else - istart = isd - endif - if (ie == npx-1) then - iend = npx-1 - else - iend = ied - endif - - do k=1,npz - do j=npy,jed - do i=istart,iend - delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - - ! North, south include all corners - if (is == 1) then - do k=1,npz - do j=npy,jed - do i=isd,0 - delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - - - if (ie == npx-1) then - do k=1,npz - do j=npy,jed - do i=npx,ied - delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - endif - endif - ! - if(nside==3)then + if(nside==2)then if(east_bc)then if (js == 1) then jstart = 1 @@ -2109,7 +2165,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==4)then + if(nside==3)then if(west_bc)then if (js == 1) then jstart = 1 @@ -2134,11 +2190,74 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif + if(nside==4)then + if(south_bc)then + if (is == 1) then + istart = 1 + else + istart = isd + endif + if (ie == npx-1) then + iend = npx-1 + else + iend = ied + endif + + do k=1,npz + do j=npy,jed + do i=istart,iend + delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + + ! North, south include all corners + if (is == 1) then + do k=1,npz + do j=npy,jed + do i=isd,0 + delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + + + if (ie == npx-1) then + do k=1,npz + do j=npy,jed + do i=npx,ied + delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + endif + endif + endif ! !----------------------------------------------------------------------- enddo sides_scalars !----------------------------------------------------------------------- + +! +!----------------------------------------------------------------------- +!*** The caller may eventually print or write Atm%ps to disk, so it +!*** cannot contain signalling NaN. Here, we replace all non-finite +!*** data with quiet NaN instead. +!----------------------------------------------------------------------- +! + do j=lbound(Atm%ps,2),ubound(Atm%ps,2) + do i=lbound(Atm%ps,1),ubound(Atm%ps,1) + if(is_not_finite(Atm%ps(i,j))) then + Atm%ps(i,j) = real_qnan + endif + enddo + enddo ! !----------------------------------------------------------------------- !*** Now that we have the pressure throughout the boundary region @@ -2214,23 +2333,6 @@ subroutine regional_bc_data(Atm,bc_hour & endif ! if(nside==2)then - if(south_bc)then - call_remap=.true. - bc_side_t1=>BC_t1%south -! - is_u=Atm%regional_bc_bounds%is_south_uvs - ie_u=Atm%regional_bc_bounds%ie_south_uvs - js_u=Atm%regional_bc_bounds%js_south_uvs - je_u=Atm%regional_bc_bounds%je_south_uvs -! - is_v=Atm%regional_bc_bounds%is_south_uvw - ie_v=Atm%regional_bc_bounds%ie_south_uvw - js_v=Atm%regional_bc_bounds%js_south_uvw - je_v=Atm%regional_bc_bounds%je_south_uvw - endif - endif -! - if(nside==3)then if(east_bc)then call_remap=.true. bc_side_t1=>BC_t1%east @@ -2247,7 +2349,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_remap=.true. bc_side_t1=>BC_t1%west @@ -2263,6 +2365,23 @@ subroutine regional_bc_data(Atm,bc_hour & je_v=Atm%regional_bc_bounds%je_west_uvw endif endif +! + if(nside==4)then + if(south_bc)then + call_remap=.true. + bc_side_t1=>BC_t1%south +! + is_u=Atm%regional_bc_bounds%is_south_uvs + ie_u=Atm%regional_bc_bounds%ie_south_uvs + js_u=Atm%regional_bc_bounds%js_south_uvs + je_u=Atm%regional_bc_bounds%je_south_uvs +! + is_v=Atm%regional_bc_bounds%is_south_uvw + ie_v=Atm%regional_bc_bounds%ie_south_uvw + js_v=Atm%regional_bc_bounds%js_south_uvw + je_v=Atm%regional_bc_bounds%je_south_uvw + endif + endif ! if(call_remap)then ! @@ -2756,17 +2875,6 @@ subroutine fill_divgd_BC endif ! if(nside==2)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%divgd_BC,1) - ie0=ubound(BC_t1%south%divgd_BC,1) - js0=lbound(BC_t1%south%divgd_BC,2) - je0=ubound(BC_t1%south%divgd_BC,2) - endif - endif -! - if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2777,7 +2885,7 @@ subroutine fill_divgd_BC endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2787,6 +2895,17 @@ subroutine fill_divgd_BC je0=ubound(BC_t1%west%divgd_BC,2) endif endif +! + if(nside==4)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%divgd_BC,1) + ie0=ubound(BC_t1%south%divgd_BC,1) + js0=lbound(BC_t1%south%divgd_BC,2) + je0=ubound(BC_t1%south%divgd_BC,2) + endif + endif ! if(call_set)then do k=1,klev_out @@ -2849,17 +2968,6 @@ subroutine fill_q_con_BC endif ! if(nside==2)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%q_con_BC,1) - ie0=ubound(BC_t1%south%q_con_BC,1) - js0=lbound(BC_t1%south%q_con_BC,2) - je0=ubound(BC_t1%south%q_con_BC,2) - endif - endif -! - if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2870,7 +2978,7 @@ subroutine fill_q_con_BC endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2880,6 +2988,17 @@ subroutine fill_q_con_BC je0=ubound(BC_t1%west%q_con_BC,2) endif endif +! + if(nside==4)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%q_con_BC,1) + ie0=ubound(BC_t1%south%q_con_BC,1) + js0=lbound(BC_t1%south%q_con_BC,2) + je0=ubound(BC_t1%south%q_con_BC,2) + endif + endif ! if(call_set)then do k=1,klev_out @@ -2937,23 +3056,23 @@ subroutine fill_cappa_BC endif ! if(nside==2)then - if(south_bc)then + if(east_bc)then call_compute=.true. - bc_side_t1=>BC_t1%south + bc_side_t1=>BC_t1%east endif endif ! if(nside==3)then - if(east_bc)then + if(west_bc)then call_compute=.true. - bc_side_t1=>BC_t1%east + bc_side_t1=>BC_t1%west endif endif ! if(nside==4)then - if(west_bc)then + if(south_bc)then call_compute=.true. - bc_side_t1=>BC_t1%west + bc_side_t1=>BC_t1%south endif endif ! @@ -3168,34 +3287,11 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==2)then - if(south_bc)then - call_get_var=.true. -! - var_name=trim(var_name_root)//"_top" -! - i_start_array=is_input - i_end_array =ie_input - j_start_array=je_input-nhalo_data-nrows_blend+1 - j_end_array =je_input -! - i_start_data=i_start_array+nhalo_data - i_count=i_end_array-i_start_array+1 - j_start_data=1 - j_count=j_end_array-j_start_array+1 -! - endif - endif -! !---------- !*** East !---------- ! - if(nside==3)then + if(nside==2)then if(east_bc)then call_get_var=.true. ! @@ -3239,7 +3335,7 @@ subroutine read_regional_bc_file(is_input,ie_input & !*** West !---------- ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_get_var=.true. ! @@ -3275,6 +3371,29 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! +!----------- +!*** South +!----------- +! + if(nside==4)then + if(south_bc)then + call_get_var=.true. +! + var_name=trim(var_name_root)//"_top" +! + i_start_array=is_input + i_end_array =ie_input + j_start_array=je_input-nhalo_data-nrows_blend+1 + j_end_array =je_input +! + i_start_data=i_start_array+nhalo_data + i_count=i_end_array-i_start_array+1 + j_start_data=1 + j_count=j_end_array-j_start_array+1 +! + endif + endif +! !----------------------------------------------------------------------- !*** Fill this task's subset of boundary data for this 3-D !*** or 4-D variable. This includes the data in the domain's @@ -3452,7 +3571,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & ,npz & !<-- # of levels in final 3-D integration variables ,ncnst !<-- # of tracer variables real, intent(in):: ak0(km+1), bk0(km+1) - real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc):: psc + real, intent(inout), dimension(is_bc:ie_bc,js_bc:je_bc):: psc real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: t_in real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: omga real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km,ncnst):: qa @@ -3480,7 +3599,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & real(kind=R_GRID), dimension(is_bc:ie_bc,km+1):: pn0 real(kind=R_GRID):: pst !!! High-precision - integer i,ie,is,j,je,js,k,l,m, k2,iq + integer i,ie,is,j,je,js,k,l,m, k2,iq, iss,ies,jss,jes, isv,iev,jsv,jev integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, cld_amt ! !--------------------------------------------------------------------------------- @@ -3514,25 +3633,102 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !*** This is needed to obtain pressures that will surround the wind points. !--------------------------------------------------------------------------------- ! - is=is_bc + isv=is_bc + iss=isv if(side=='west')then - is=ie_bc-nhalo_data-nrows_blend+1 + isv=ie_bc-nhalo_data-nrows_blend + iss=isv+1 endif ! - ie=ie_bc + iev=ie_bc + ies=iev if(side=='east')then - ie=is_bc+nhalo_data+nrows_blend-1 + iev=is_bc+nhalo_data+nrows_blend + ies=iev-1 endif ! - js=js_bc + jsv=js_bc + jss=jsv if(side=='south')then - js=je_bc-nhalo_data-nrows_blend+1 + jsv=je_bc-nhalo_data-nrows_blend + jss=jsv+1 endif ! - je=je_bc + jev=je_bc + jes=jev if(side=='north')then - je=js_bc+nhalo_data+nrows_blend-1 + jev=js_bc+nhalo_data+nrows_blend + jes=jev-1 + endif +! +!--------------------------------------------------------------------------------- +!*** Workaround for incomplete 'ps' in boundary condition files. Copy from +!*** scalar halo regions (s) to velocity halo regions (v) beyond them. +!--------------------------------------------------------------------------------- +! + if(side=='west')then + do j=jss,jes + psc(isv,j) = psc(iss,j) + enddo + endif +! + if(side=='east')then + do j=jss,jes + psc(iev,j) = psc(ies,j) + enddo + endif +! + if(side=='south')then + ! Special handling of corners + if(west_bc) then + is=iss + psc(isv,jsv) = psc(iss,jss) + else + is=isv + endif + if(east_bc) then + ie=ies + psc(iev,jsv) = psc(ies,jss) + else + ie=iev + endif + ! The rest of the south boundary + do i=is,ie + psc(i,jsv) = psc(i,jss) + enddo endif +! + if(side=='north')then + ! Special handling of corners + if(west_bc) then + is=iss + psc(isv,jev) = psc(iss,jes) + else + is=isv + endif + if(east_bc) then + ie=ies + psc(iev,jev) = psc(ies,jes) + else + ie=iev + endif + ! The rest of the north boundary + do i=is,ie + psc(i,jev) = psc(i,jes) + enddo + endif +! + is = iss + js = jss + ie = ies + je = jes +! +! Ensure uninitialized memory isn't used + pn0 = real_snan + pn1 = real_snan + gz_fv = real_snan + gz = real_snan + pn = real_snan ! allocate(pe0(is:ie,km+1)) ; pe0=real_snan allocate(qn1(is:ie,npz)) ; qn1=real_snan @@ -3564,29 +3760,87 @@ subroutine remap_scalar_nggps_regional_bc(Atm & pn(k) = 2.*pn(km+1) - pn(l) enddo + pst = real_snan do k=km+k2-1, 2, -1 if( phis_reg(i,j).le.gz(k) .and. phis_reg(i,j).ge.gz(k+1) ) then pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-phis_reg(i,j))/(gz(k)-gz(k+1)) - go to 123 + exit endif enddo - 123 ps(i,j) = exp(pst) + ps(i,j) = exp(pst) enddo ! i-loop !--------------------------------------------------------------------------------- enddo jloop1 !--------------------------------------------------------------------------------- +! +!--------------------------------------------------------------------------------- +!*** Workaround for incomplete 'ps' in boundary condition files. Copy from +!*** scalar halo regions (s) to velocity halo regions (v) beyond them. +!--------------------------------------------------------------------------------- +! + if(side=='west')then + do j=jss,jes + ps(isv,j) = ps(iss,j) + enddo + endif +! + if(side=='east')then + do j=jss,jes + ps(iev,j) = ps(ies,j) + enddo + endif +! + if(side=='south')then + ! Special handling of corners + if(west_bc) then + is=iss + ps(isv,jsv) = ps(iss,jss) + else + is=isv + endif + if(east_bc) then + ie=ies + ps(iev,jsv) = ps(ies,jss) + else + ie=iev + endif + ! The rest of the south boundary + do i=is,ie + ps(i,jsv) = ps(i,jss) + enddo + endif +! + if(side=='north')then + ! Special handling of corners + if(west_bc) then + is=iss + ps(isv,jev) = ps(iss,jes) + else + is=isv + endif + if(east_bc) then + ie=ies + ps(iev,jev) = ps(ies,jes) + else + ie=iev + endif + ! The rest of the north boundary + do i=is,ie + ps(i,jev) = ps(i,jes) + enddo + endif !--------------------------------------------------------------------------------- !*** Transfer values from the expanded boundary array for sfc pressure into !*** the Atm object. !--------------------------------------------------------------------------------- ! - is=lbound(Atm%ps,1) - ie=ubound(Atm%ps,1) - js=lbound(Atm%ps,2) - je=ubound(Atm%ps,2) + is=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),isv)) + ie=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),iev)) + js=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jsv)) + je=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jev)) ! do j=js,je do i=is,ie @@ -3923,12 +4177,22 @@ subroutine remap_dwinds_regional_bc(Atm & !------ do k=1,km+1 do i=is_u,ie_u - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) + if(is_not_finite(psc(i,j-1))) then + ! Workaround for bug: PS not available in some velocity areas + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + else + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) + endif enddo enddo do k=1,npz+1 do i=is_u,ie_u - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) + if(is_not_finite(psc(i,j-1))) then + ! Workaround for bug: PS not available in some velocity areas + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + else + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) + endif enddo enddo call mappm(km, pe0(is_u:ie_u,1:km+1), ud(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), & @@ -3964,12 +4228,22 @@ subroutine remap_dwinds_regional_bc(Atm & do k=1,km+1 do i=is_v,ie_v - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) + if(is_not_finite(psc(i-1,j))) then + ! Workaround for bug: PS not available in some velocity areas + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + else + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) + endif enddo enddo do k=1,npz+1 do i=is_v,ie_v - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) + if(is_not_finite(psc(i-1,j))) then + ! Workaround for bug: PS not available in some velocity areas + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + else + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) + endif enddo enddo call mappm(km, pe0(is_v:ie_v,1:km+1), vd(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), & @@ -4352,7 +4626,7 @@ subroutine regional_boundary_update(array & ! real,dimension(:,:,:),pointer :: bc_t0,bc_t1 !<-- Boundary data at the two bracketing times. ! - logical :: blend,call_interp + logical :: blend,call_interp,blendtmp ! !--------------------------------------------------------------------- !********************************************************************* @@ -4396,13 +4670,21 @@ subroutine regional_boundary_update(array & i2=ied+1 endif ! - j1=jsd - j2=js-1 + j1=jsd ! -2 -- outermost boundary ghost zone + j2=js-1 ! 0 -- first boundary ghost zone ! + IF ( east_bc ) THEN i1_blend=is + ELSE + i1_blend=isd !is-nhalo_model + ENDIF + IF ( west_bc ) THEN i2_blend=ie + ELSE + i2_blend=ied ! ie+nhalo_model + ENDIF if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then - i2_blend=ie+1 + i2_blend=i2_blend+1 ! ie+1 endif j1_blend=js j2_blend=js+nrows_blend_user-1 @@ -4412,50 +4694,11 @@ subroutine regional_boundary_update(array & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==2)then - if(south_bc)then - call_interp=.true. - bc_side_t0=>bc_south_t0 - bc_side_t1=>bc_south_t1 -! - i1=isd - i2=ied - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then - i2=ied+1 - endif -! - j1=je+1 - j2=jed - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j1=je+2 - j2=jed+1 - endif -! - i1_blend=is - i2_blend=ie - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then - i2_blend=ie+1 - endif - j2_blend=je - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j2_blend=je+1 - endif - j1_blend=j2_blend-nrows_blend_user+1 - i_bc=-9e9 - j_bc=j1 -! - endif - endif -! !---------- !*** East !---------- ! - if(nside==3)then + if(nside==2)then if(east_bc)then call_interp=.true. bc_side_t0=>bc_east_t0 @@ -4483,16 +4726,21 @@ subroutine regional_boundary_update(array & endif endif ! - i1_blend=is - i2_blend=is+nrows_blend_user-1 +! Note: Original code checked for corner region and avoided overlap, but changed this to blend corners from both boundaries + i1_blend=is + i2_blend=is+nrows_blend_user-1 + + IF ( north_bc ) THEN j1_blend=js + ELSE + j1_blend=jsd !js-nhalo_model + ENDIF + IF ( south_bc ) THEN j2_blend=je - if(north_bc)then - j1_blend=js+nrows_blend_user !<-- North BC already handles nrows_blend_user blending rows - endif - if(south_bc)then - j2_blend=je-nrows_blend_user !<-- South BC already handles nrows_blend_user blending rows - endif + ELSE + j2_blend=jed ! ie+nhalo_model + ENDIF + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then j2_blend=j2_blend+1 endif @@ -4506,7 +4754,7 @@ subroutine regional_boundary_update(array & !*** West !---------- ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_interp=.true. bc_side_t0=>bc_west_t0 @@ -4538,16 +4786,20 @@ subroutine regional_boundary_update(array & endif endif ! +! Note: Original code checked for corner region and avoided overlap, but changed this to blend corners from both boundaries i1_blend=i1-nrows_blend_user i2_blend=i1-1 + + IF ( north_bc ) THEN j1_blend=js + ELSE + j1_blend=jsd !is-nhalo_model + ENDIF + IF ( south_bc ) THEN j2_blend=je - if(north_bc)then - j1_blend=js+nrows_blend_user !<-- North BC already handled nrows_blend_user blending rows. - endif - if(south_bc)then - j2_blend=je-nrows_blend_user !<-- South BC already handled nrows_blend_user blending rows. - endif + ELSE + j2_blend=jed ! ie+nhalo_model + ENDIF if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then j2_blend=j2_blend+1 endif @@ -4557,12 +4809,64 @@ subroutine regional_boundary_update(array & endif endif ! +!----------- +!*** South +!----------- +! + if(nside==4)then + if(south_bc)then + call_interp=.true. + bc_side_t0=>bc_south_t0 + bc_side_t1=>bc_south_t1 +! + i1=isd + i2=ied + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then + i2=ied+1 + endif +! + j1=je+1 + j2=jed + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j1=je+2 + j2=jed+1 + endif +! + i1_blend=is + i2_blend=ie + IF ( east_bc ) THEN + i1_blend=is + ELSE + i1_blend=isd !is-nhalo_model + ENDIF + IF ( west_bc ) THEN + i2_blend=ie + ELSE + i2_blend=ied ! ie+nhalo_model + ENDIF + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then +! i2_blend=ie+1 + i2_blend=i2_blend+1 + endif + j2_blend=je + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j2_blend=je+1 + endif + j1_blend=j2_blend-nrows_blend_user+1 + i_bc=-9e9 + j_bc=j1 +! + endif + endif +! !--------------------------------------------------------------------- !*** Get the pointers pointing at the boundary arrays holding the !*** two time levels of the given prognostic array's boundary region !*** then update the boundary points. !--------------------------------------------------------------------- ! + + if(call_interp)then ! call retrieve_bc_variable_data(bc_vbl_name & @@ -4787,7 +5091,7 @@ subroutine bc_time_interpolation(array & ! !--------------------------------------------------------------------- -! +! Set values in the boundary points only do k=1,ubnd_z do j=j1,j2 do i=i1,i2 @@ -4818,11 +5122,14 @@ subroutine bc_time_interpolation(array & !----------- ! if(nside==1.and.north_bc)then - rdenom=1./real(j2_blend-j_bc-1) + rdenom=1./real(Max(1,j2_blend-j_bc-1)) do k=1,ubnd_z do j=j1_blend,j2_blend factor_dist=exp(-(blend_exp1+blend_exp2*(j-j_bc-1)*rdenom)) !<-- Exponential falloff of blending weights. do i=i1_blend,i2_blend + if(is_not_finite(array(i,j,k))) then + cycle ! Outside boundary + endif blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1. ! @@ -4832,34 +5139,19 @@ subroutine bc_time_interpolation(array & enddo endif ! -!----------- -!*** South -!----------- -! - if(nside==2.and.south_bc)then - rdenom=1./real(j_bc-j1_blend-1) - do k=1,ubnd_z - do j=j1_blend,j2_blend - factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom)) !<-- Exponential falloff of blending weights. - do i=i1_blend,i2_blend - blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated - +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1. - array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value - enddo - enddo - enddo - endif -! !---------- !*** East !---------- ! - if(nside==3.and.east_bc)then - rdenom=1./real(i2_blend-i_bc-1) + if(nside==2.and.east_bc)then + rdenom=1./real(Max(1,i2_blend-i_bc-1)) do k=1,ubnd_z do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! + if(is_not_finite(array(i,j,k))) then + cycle ! Outside boundary + endif blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1. ! @@ -4875,12 +5167,15 @@ subroutine bc_time_interpolation(array & !*** West !---------- ! - if(nside==4.and.west_bc)then - rdenom=1./real(i_bc-i1_blend-1) + if(nside==3.and.west_bc)then + rdenom=1./real(Max(1, i_bc-i1_blend-1)) do k=1,ubnd_z do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! + if(is_not_finite(array(i,j,k))) then + cycle ! Outside boundary + endif blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1. ! @@ -4892,6 +5187,27 @@ subroutine bc_time_interpolation(array & enddo endif ! +!----------- +!*** South +!----------- +! + if(nside==4.and.south_bc)then + rdenom=1./real(Max(1,j_bc-j1_blend-1)) + do k=1,ubnd_z + do j=j1_blend,j2_blend + factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom)) !<-- Exponential falloff of blending weights. + do i=i1_blend,i2_blend + if(is_not_finite(array(i,j,k))) then + cycle ! Outside boundary + endif + blend_value=bc_t0(i,j,k) & !<-- Blend data interpolated + +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1. + array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value + enddo + enddo + enddo + endif +! !--------------------------------------------------------------------- ! end subroutine bc_time_interpolation @@ -4970,29 +5286,6 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif ! if(nside==2)then - if(south_bc)then - move=.true. - bc_side_t0=>BC_t0%south - bc_side_t1=>BC_t1%south -! - is_c=bnds%is_south !<-- - ie_c=bnds%ie_south ! South BC index limits - js_c=bnds%js_south ! for centers of grid cells. - je_c=bnds%je_south !<-- -! - is_s=bnds%is_south_uvs !<-- - ie_s=bnds%ie_south_uvs ! South BC index limits - js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. - je_s=bnds%je_south_uvs !<-- -! - is_w=bnds%is_south_uvw !<-- - ie_w=bnds%ie_south_uvw ! South BC index limits - js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. - je_w=bnds%je_south_uvw !<-- - endif - endif -! - if(nside==3)then if(east_bc)then move=.true. bc_side_t0=>BC_t0%east @@ -5015,7 +5308,7 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then move=.true. bc_side_t0=>BC_t0%west @@ -5037,6 +5330,29 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & je_w=bnds%je_west_uvw !<-- endif endif +! + if(nside==4)then + if(south_bc)then + move=.true. + bc_side_t0=>BC_t0%south + bc_side_t1=>BC_t1%south +! + is_c=bnds%is_south !<-- + ie_c=bnds%ie_south ! South BC index limits + js_c=bnds%js_south ! for centers of grid cells. + je_c=bnds%je_south !<-- +! + is_s=bnds%is_south_uvs !<-- + ie_s=bnds%ie_south_uvs ! South BC index limits + js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. + je_s=bnds%je_south_uvs !<-- +! + is_w=bnds%is_south_uvw !<-- + ie_w=bnds%ie_south_uvw ! South BC index limits + js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. + je_w=bnds%je_south_uvw !<-- + endif + endif ! if(move)then do k=1,nlev @@ -5686,7 +6002,7 @@ subroutine dump_field_3d (domain, name, field, isd, ied, jsd, jed, nlev, stag) nyg = npy + 2*halo + jext nz = size(field,dim=3) - allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) ) + allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) ) ; glob_field=real_snan isection_s = is isection_e = ie @@ -5805,7 +6121,7 @@ subroutine dump_field_2d (domain, name, field, isd, ied, jsd, jed, stag) nxg = npx + 2*halo + iext nyg = npy + 2*halo + jext - allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) ) + allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) ) ; glob_field=real_snan isection_s = is isection_e = ie @@ -6113,6 +6429,7 @@ subroutine write_full_fields(Atm) nz=size(fields_core(nv)%ptr,3) ! allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) ) + global_field=real_snan ! !----------------------------------------------------------------------- !*** What is the local extent of the variable on the task subdomain? @@ -6217,6 +6534,7 @@ subroutine write_full_fields(Atm) nz=size(fields_tracers(1)%ptr,3) ! allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) ) + global_field=real_snan ! !----------------------------------------------------------------------- !*** What is the local extent of the variable on the task subdomain? @@ -6573,10 +6891,10 @@ subroutine exch_uv(domain, bd, npz, u, v) ! buf1 and buf4 must be of the same size (sim. for buf2 and buf3). ! Changes to the code below should be tested with debug flags ! enabled (out-of-bounds reads/writes). - allocate(buf1(1:24*npz)) - allocate(buf2(1:36*npz)) - allocate(buf3(1:36*npz)) - allocate(buf4(1:24*npz)) + allocate(buf1(1:24*npz)) ; buf1=real_snan + allocate(buf2(1:36*npz)) ; buf2=real_snan + allocate(buf3(1:36*npz)) ; buf3=real_snan + allocate(buf4(1:24*npz)) ; buf4=real_snan ! FIXME: MPI_COMM_WORLD diff --git a/model/fv_sg.F90 b/model/fv_sg.F90 index 2446b1144..e60693203 100644 --- a/model/fv_sg.F90 +++ b/model/fv_sg.F90 @@ -404,12 +404,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat) enddo enddo - else + elseif ( nwat==6 ) then do k=1,kbot do i=is,ie qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel) enddo enddo + elseif ( nwat==7 ) then + do k=1,kbot + do i=is,ie + qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)+q0(i,k,hailwat) + enddo + enddo endif do k=kbot, 2, -1 @@ -985,12 +991,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & do i=is,ie qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat) enddo - else + elseif ( nwat==6 ) then do k=1,kbot do i=is,ie qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel) enddo enddo + elseif ( nwat==7 ) then + do k=1,kbot + do i=is,ie + qcon(i,k) = q0(i,k,liq_wat)+q0(i,k,ice_wat)+q0(i,k,snowwat)+q0(i,k,rainwat)+q0(i,k,graupel)+q0(i,k,hailwat) + enddo + enddo endif do k=kbot, 2, -1 diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 56cd554f7..6be85a45f 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -202,14 +202,13 @@ module external_ic_mod ! Include variable "version" to be written to log file. #include - public get_external_ic, get_cubed_sphere_terrain + public get_external_ic, get_cubed_sphere_terrain, remap_scalar, remap_dwinds contains - subroutine get_external_ic( Atm, fv_domain, cold_start ) + subroutine get_external_ic( Atm, cold_start ) type(fv_atmos_type), intent(inout), target :: Atm - type(domain2d), intent(inout) :: fv_domain logical, intent(IN) :: cold_start real:: alpha = 0. real rdg @@ -260,14 +259,14 @@ subroutine get_external_ic( Atm, fv_domain, cold_start ) enddo enddo - call mpp_update_domains( f0, fv_domain ) + call mpp_update_domains( f0, Atm%domain ) if ( Atm%gridstruct%cubed_sphere .and. (.not. Atm%gridstruct%bounded_domain))then call fill_corners(f0, Atm%npx, Atm%npy, YDir) endif ! Read in cubed_sphere terrain if ( Atm%flagstruct%mountain ) then - call get_cubed_sphere_terrain(Atm, fv_domain) + call get_cubed_sphere_terrain(Atm) else if (.not. Atm%neststruct%nested) Atm%phis = 0. !TODO: Not sure about this line --- lmh 30 may 18 endif @@ -276,32 +275,32 @@ subroutine get_external_ic( Atm, fv_domain, cold_start ) if ( Atm%flagstruct%ncep_ic ) then nq = 1 call timing_on('NCEP_IC') - call get_ncep_ic( Atm, fv_domain, nq ) + call get_ncep_ic( Atm, nq ) call timing_off('NCEP_IC') #ifdef FV_TRACERS if (.not. cold_start) then - call fv_io_read_tracers( fv_domain, Atm ) + call fv_io_read_tracers( Atm ) if(is_master()) write(*,*) 'All tracers except sphum replaced by FV IC' endif #endif elseif ( Atm%flagstruct%nggps_ic ) then call timing_on('NGGPS_IC') - call get_nggps_ic( Atm, fv_domain ) + call get_nggps_ic( Atm ) call timing_off('NGGPS_IC') elseif ( Atm%flagstruct%hrrrv3_ic ) then call timing_on('HRRR_IC') - call get_hrrr_ic( Atm, fv_domain ) + call get_hrrr_ic( Atm ) call timing_off('HRRR_IC') elseif ( Atm%flagstruct%ecmwf_ic ) then if( is_master() ) write(*,*) 'Calling get_ecmwf_ic' call timing_on('ECMWF_IC') - call get_ecmwf_ic( Atm, fv_domain ) + call get_ecmwf_ic( Atm ) call timing_off('ECMWF_IC') else ! The following is to read in legacy lat-lon FV core restart file ! is Atm%q defined in all cases? nq = size(Atm%q,4) - call get_fv_ic( Atm, fv_domain, nq ) + call get_fv_ic( Atm, nq ) endif call prt_maxmin('PS', Atm%ps, is, ie, js, je, ng, 1, 0.01) @@ -368,9 +367,8 @@ end subroutine get_external_ic !------------------------------------------------------------------ - subroutine get_cubed_sphere_terrain( Atm, fv_domain ) + subroutine get_cubed_sphere_terrain( Atm ) type(fv_atmos_type), intent(inout), target :: Atm - type(domain2d), intent(inout) :: fv_domain type(FmsNetcdfDomainFile_t) :: Fv_core integer :: tile_id(1) character(len=64) :: fname @@ -393,13 +391,13 @@ subroutine get_cubed_sphere_terrain( Atm, fv_domain ) jed = Atm%bd%jed ng = Atm%bd%ng - tile_id = mpp_get_tile_id( fv_domain ) + tile_id = mpp_get_tile_id( Atm%domain ) fname = 'INPUT/fv_core.res.nc' call mpp_error(NOTE, 'external_ic: looking for '//fname) - if( open_file(Fv_core, fname, "read", fv_domain, is_restart=.true.) ) then + if( open_file(Fv_core, fname, "read", Atm%domain_for_read, is_restart=.true.) ) then call read_data(Fv_core, 'phis', Atm%phis(is:ie,js:je)) call close_file(Fv_core) else @@ -428,7 +426,7 @@ end subroutine get_cubed_sphere_terrain !>@brief The subroutine 'get_nggps_ic' reads in data after it has been preprocessed with !! NCEP/EMC orography maker and 'global_chgres', and has been horiztontally !! interpolated to the current cubed-sphere grid - subroutine get_nggps_ic (Atm, fv_domain) + subroutine get_nggps_ic (Atm) !>variables read in from 'gfs_ctrl.nc' !> VCOORD - level information @@ -456,7 +454,6 @@ subroutine get_nggps_ic (Atm, fv_domain) #endif type(fv_atmos_type), intent(inout) :: Atm - type(domain2d), intent(inout) :: fv_domain ! local: real, dimension(:), allocatable:: ak, bk real, dimension(:,:), allocatable:: wk2, ps, oro_g @@ -585,7 +582,7 @@ subroutine get_nggps_ic (Atm, fv_domain) !--- read in surface temperature (k) and land-frac ! surface skin temperature - if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then naxis_dims = get_variable_num_dimensions(SFC_restart, 'tsea') allocate (dim_names_alloc(naxis_dims)) call get_variable_dimension_names(SFC_restart, 'tsea', dim_names_alloc) @@ -605,7 +602,7 @@ subroutine get_nggps_ic (Atm, fv_domain) dim_names_2d(2) = "lon" ! terrain surface height -- (needs to be transformed into phis = zs*grav) - if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(ORO_restart, "lat", "y") call register_axis(ORO_restart, "lon", "x") if (filtered_terrain) then @@ -915,7 +912,7 @@ subroutine read_gfs_ic() dim_names_3d4(1) = "levp" ! surface pressure (Pa) - if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(GFS_restart, "lat", "y") call register_axis(GFS_restart, "lon", "x") call register_axis(GFS_restart, "lonp", "x", domain_position=east) @@ -945,7 +942,7 @@ subroutine read_gfs_ic() do nt = 1, ntracers q(:,:,:,nt) = -999.99 - + call get_tracer_names(MODEL_ATMOS, nt, tracer_name) call register_restart_field(GFS_restart, trim(tracer_name), q(:,:,:,nt), dim_names_3d3, is_optional=.true.) enddo @@ -964,7 +961,7 @@ subroutine read_gfs_ic() end subroutine get_nggps_ic !------------------------------------------------------------------ !------------------------------------------------------------------ - subroutine get_hrrr_ic (Atm, fv_domain) + subroutine get_hrrr_ic (Atm) ! read in data after it has been preprocessed with ! NCEP/EMC orography maker ! @@ -990,7 +987,6 @@ subroutine get_hrrr_ic (Atm, fv_domain) type(fv_atmos_type), intent(inout) :: Atm - type(domain2d), intent(inout) :: fv_domain ! local: real, dimension(:), allocatable:: ak, bk real, dimension(:,:), allocatable:: wk2, ps, oro_g @@ -1104,7 +1100,7 @@ subroutine get_hrrr_ic (Atm, fv_domain) !--- read in surface temperature (k) and land-frac ! surface skin temperature - if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(SFC_restart, fn_sfc_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call get_variable_dimension_names(SFC_restart, 'tsea', dim_names_2d) call register_axis(SFC_restart, dim_names_2d(2), "y") call register_axis(SFC_restart, dim_names_2d(1), "x") @@ -1121,7 +1117,7 @@ subroutine get_hrrr_ic (Atm, fv_domain) dim_names_2d(2) = "lon" ! terrain surface height -- (needs to be transformed into phis = zs*grav) - if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(ORO_restart, "lat", "y") call register_axis(ORO_restart, "lon", "x") if (filtered_terrain) then @@ -1164,7 +1160,7 @@ subroutine get_hrrr_ic (Atm, fv_domain) dim_names_3d4(1) = "levp" ! edge pressure (Pa) - if( open_file(HRRR_restart, fn_hrr_ics, "read", Atm%domain,is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(HRRR_restart, fn_hrr_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(HRRR_restart, "lat", "y") call register_axis(HRRR_restart, "lon", "x") call register_axis(HRRR_restart, "lonp", "x", domain_position=east) @@ -1360,9 +1356,8 @@ end subroutine get_hrrr_ic !------------------------------------------------------------------ !------------------------------------------------------------------ !>@brief The subroutine 'get_ncep_ic' reads in the specified NCEP analysis or reanalysis dataset - subroutine get_ncep_ic( Atm, fv_domain, nq ) + subroutine get_ncep_ic( Atm, nq ) type(fv_atmos_type), intent(inout) :: Atm - type(domain2d), intent(inout) :: fv_domain integer, intent(in):: nq ! local: #ifdef HIWPP_ETA @@ -1818,7 +1813,7 @@ end subroutine get_ncep_ic !>@brief The subroutine 'get_ecmwf_ic' reads in initial conditions from ECMWF analyses !! (EXPERIMENTAL: contact Jan-Huey Chen jan-huey.chen@noaa.gov for support) !>@authors Jan-Huey Chen, Xi Chen, Shian-Jiann Lin - subroutine get_ecmwf_ic( Atm, fv_domain ) + subroutine get_ecmwf_ic( Atm ) #ifdef __PGI use GFS_restart, only : GFS_restart_type @@ -1827,7 +1822,6 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) #endif type(fv_atmos_type), intent(inout) :: Atm - type(domain2d), intent(inout) :: fv_domain ! local: real :: ak_ec(138), bk_ec(138) data ak_ec/ 0.000000, 2.000365, 3.102241, 4.666084, 6.827977, 9.746966, & @@ -2062,7 +2056,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) dim_names_3d4(1) = "levp" !! Read in model terrain from oro_data.tile?.nc - if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(ORO_restart, fn_oro_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(ORO_restart, "lat", "y") call register_axis(ORO_restart, "lon", "x") if (filtered_terrain) then @@ -2088,7 +2082,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) allocate (ps_gfs(is:ie,js:je)) allocate (zh_gfs(is:ie,js:je,levp_gfs+1)) - if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain, is_restart=.true., dont_add_res_to_filename=.true.) ) then + if( open_file(GFS_restart, fn_gfs_ics, "read", Atm%domain_for_read, is_restart=.true., dont_add_res_to_filename=.true.) ) then call register_axis(GFS_restart, "lat", "y") call register_axis(GFS_restart, "lon", "x") call register_axis(GFS_restart, "levp", size(zh_gfs,3)) @@ -2622,9 +2616,8 @@ subroutine get_ecmwf_ic( Atm, fv_domain ) end subroutine get_ecmwf_ic !------------------------------------------------------------------ !------------------------------------------------------------------ - subroutine get_fv_ic( Atm, fv_domain, nq ) + subroutine get_fv_ic( Atm, nq ) type(fv_atmos_type), intent(inout) :: Atm - type(domain2d), intent(inout) :: fv_domain integer, intent(in):: nq type(FmsNetcdfFile_t) :: Latlon_dyn, Latlon_tra diff --git a/tools/fv_io.F90 b/tools/fv_io.F90 index e892b6811..fd5294b32 100644 --- a/tools/fv_io.F90 +++ b/tools/fv_io.F90 @@ -419,7 +419,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) Atm(1)%Fv_restart_tile_is_open = open_file(Atm(1)%Fv_restart_tile, fname, "read", fv_domain, is_restart=.true.) if (Atm(1)%Fv_restart_tile_is_open) then call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Fv_restart_tile) + call read_restart(Atm(1)%Fv_restart_tile, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Fv_restart_tile) Atm(1)%Fv_restart_tile_is_open = .false. endif @@ -429,7 +429,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) Atm(1)%Tra_restart_is_open = open_file(Atm(1)%Tra_restart, fname, "read", fv_domain, is_restart=.true.) if (Atm(1)%Tra_restart_is_open) then call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Tra_restart) + call read_restart(Atm(1)%Tra_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Tra_restart) Atm(1)%Tra_restart_is_open = .false. else @@ -442,7 +442,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) if (Atm(1)%Rsf_restart_is_open) then Atm(1)%flagstruct%srf_init = .true. call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Rsf_restart) + call read_restart(Atm(1)%Rsf_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Rsf_restart) Atm(1)%Rsf_restart_is_open = .false. else @@ -456,7 +456,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) Atm(1)%Mg_restart_is_open = open_file(Atm(1)%Mg_restart, fname, "read", fv_domain, is_restart=.true.) if (Atm(1)%Mg_restart_is_open) then call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Mg_restart) + call read_restart(Atm(1)%Mg_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Mg_restart) Atm(1)%Mg_restart_is_open = .false. else @@ -467,7 +467,7 @@ subroutine fv_io_read_restart(fv_domain,Atm) Atm(1)%Lnd_restart_is_open = open_file(Atm(1)%Lnd_restart, fname, "read", fv_domain, is_restart=.true.) if (Atm(1)%Lnd_restart_is_open) then call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Lnd_restart) + call read_restart(Atm(1)%Lnd_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Lnd_restart) Atm(1)%Lnd_restart_is_open = .false. else @@ -483,8 +483,7 @@ end subroutine fv_io_read_restart !>@details This subroutine is useful when initializing a cycled or nudged model !! from an analysis that does not have a whole set of microphysical, aerosol, or !! chemical tracers - subroutine fv_io_read_tracers(fv_domain,Atm) - type(domain2d), intent(inout) :: fv_domain + subroutine fv_io_read_tracers(Atm) type(fv_atmos_type), intent(inout) :: Atm(:) integer :: ntracers, ntprog, nt, isc, iec, jsc, jec character(len=6) :: stile_name @@ -499,7 +498,7 @@ subroutine fv_io_read_tracers(fv_domain,Atm) call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers, num_prog=ntprog) ! fix for single tile runs where you need fv_core.res.nc and fv_core.res.tile1.nc - ntiles = mpp_get_ntile_count(fv_domain) + ntiles = mpp_get_ntile_count(Atm(1)%domain_for_read) if(ntiles == 1 .and. .not. Atm(1)%neststruct%nested) then stile_name = '.tile1' else @@ -508,7 +507,7 @@ subroutine fv_io_read_tracers(fv_domain,Atm) fname = 'INPUT/fv_tracer.res'//trim(stile_name)//'.nc' - if (open_file(Tra_restart_r,fname,"read",fv_domain, is_restart=.true.)) then + if (open_file(Tra_restart_r,fname,"read",Atm(1)%domain_for_read, is_restart=.true.)) then do nt = 2, ntprog call get_tracer_names(MODEL_ATMOS, nt, tracer_name) call set_tracer_profile (MODEL_ATMOS, nt, Atm(1)%q(isc:iec,jsc:jec,:,nt) ) @@ -534,10 +533,9 @@ end subroutine fv_io_read_tracers !>@brief The subroutine 'remap_restart' remaps the model state from remap files !! to a new set of Eulerian coordinates. !>@details Use if npz (run time z-dimension) /= npz_rst (restart z-dimension) - subroutine remap_restart(fv_domain,Atm) + subroutine remap_restart(Atm) use fv_mapz_mod, only: rst_remap - type(domain2d), intent(inout) :: fv_domain type(fv_atmos_type), intent(inout) :: Atm(:) character(len=64) :: fname, tracer_name @@ -546,6 +544,7 @@ subroutine remap_restart(fv_domain,Atm) integer :: isd, ied, jsd, jed integer :: ntiles + type(domain2d) :: fv_domain type(FmsNetcdfDomainFile_t) :: FV_tile_restart_r, Tra_restart_r type(FmsNetcdfFile_t) :: Fv_restart_r integer, allocatable, dimension(:) :: pes !< Array of the pes in the current pelist @@ -560,6 +559,7 @@ subroutine remap_restart(fv_domain,Atm) integer npz, npz_rst, ng integer i,j,k + fv_domain = Atm(1)%domain_for_read npz = Atm(1)%npz ! run time z dimension npz_rst = Atm(1)%flagstruct%npz_rst ! restart z dimension isc = Atm(1)%bd%isc; iec = Atm(1)%bd%iec; jsc = Atm(1)%bd%jsc; jec = Atm(1)%bd%jec @@ -640,7 +640,7 @@ subroutine remap_restart(fv_domain,Atm) if (Atm(1)%Rsf_restart_is_open) then Atm%flagstruct%srf_init = .true. call fv_io_register_restart(Atm(1)) - call read_restart(Atm(1)%Rsf_restart) + call read_restart(Atm(1)%Rsf_restart, ignore_checksum=Atm(1)%flagstruct%ignore_rst_cksum) call close_file(Atm(1)%Rsf_restart) Atm(1)%Rsf_restart_is_open = .false. call mpp_error(NOTE,'==> Warning from remap_restart: Expected file '//trim(fname)//' does not exist') @@ -1279,12 +1279,12 @@ subroutine fv_io_read_BCs(Atm) call fv_io_register_restart_BCs(Atm) if (Atm%neststruct%BCfile_sw_is_open) then - call read_restart_bc(Atm%neststruct%BCfile_sw) + call read_restart_bc(Atm%neststruct%BCfile_sw, ignore_checksum=Atm%flagstruct%ignore_rst_cksum) call close_file(Atm%neststruct%BCfile_sw) endif if (Atm%neststruct%BCfile_ne_is_open) then - call read_restart_bc(Atm%neststruct%BCfile_ne) + call read_restart_bc(Atm%neststruct%BCfile_ne, ignore_checksum=Atm%flagstruct%ignore_rst_cksum) call close_file(Atm%neststruct%BCfile_ne) endif diff --git a/tools/fv_mp_mod.F90 b/tools/fv_mp_mod.F90 index b6f279e31..4f2c0f2c2 100644 --- a/tools/fv_mp_mod.F90 +++ b/tools/fv_mp_mod.F90 @@ -91,7 +91,7 @@ module fv_mp_mod use mpp_domains_mod, only : mpp_group_update_initialized, mpp_do_group_update use mpp_domains_mod, only : mpp_create_group_update,mpp_reset_group_update_field use mpp_domains_mod, only : group_halo_update_type => mpp_group_update_type - use mpp_domains_mod, only: nest_domain_type + use mpp_domains_mod, only : nest_domain_type, mpp_get_io_domain_layout, mpp_get_layout, mpp_copy_domain use mpp_parameter_mod, only : WUPDATE, EUPDATE, SUPDATE, NUPDATE, XUPDATE, YUPDATE use fv_arrays_mod, only: fv_atmos_type, fv_grid_bounds_type use mpp_mod, only : mpp_get_current_pelist, mpp_set_current_pelist @@ -367,7 +367,7 @@ end subroutine mp_stop ! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ! !>@brief The subroutine 'domain_decomp' sets up the domain decomposition. subroutine domain_decomp(grid_num,npx,npy,nregions,grid_type,nested,layout,io_layout,bd,tile,square_domain,& - npes_per_tile,domain,domain_for_coupler,num_contact,pelist) + npes_per_tile,domain,domain_for_coupler,domain_for_read,num_contact,pelist) integer, intent(IN) :: grid_num integer, intent(IN) :: npx,npy,grid_type @@ -389,8 +389,9 @@ subroutine domain_decomp(grid_num,npx,npy,nregions,grid_type,nested,layout,io_la integer, intent(INOUT) :: pelist(:) integer, intent(OUT) :: num_contact, npes_per_tile logical, intent(OUT) :: square_domain - type(domain2D), intent(OUT) :: domain, domain_for_coupler + type(domain2D), intent(OUT) :: domain, domain_for_coupler, domain_for_read type(fv_grid_bounds_type), intent(INOUT) :: bd + integer :: l_layout(2) nx = npx-1 ny = npy-1 @@ -652,6 +653,17 @@ subroutine domain_decomp(grid_num,npx,npy,nregions,grid_type,nested,layout,io_la call mpp_define_io_domain(domain, io_layout) call mpp_define_io_domain(domain_for_coupler, io_layout) + !--- create a read domain that can be used to improve read performance + !--- if io_layout=(1,1) then read io_layout=layout (all read) + !--- if io_layout\=(1,1) then read io_layout=io_layout (no change) + l_layout = mpp_get_io_domain_layout(domain) + call mpp_copy_domain(domain, domain_for_read) + if (ALL(l_layout == 1)) then + call mpp_get_layout(domain, l_layout) + call mpp_define_io_domain(domain_for_read, l_layout) + else + call mpp_define_io_domain(domain_for_read, l_layout) + endif endif deallocate(pe_start,pe_end) diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index bcd6ea13b..313ef7355 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -343,7 +343,7 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this !3. External_ic if (Atm(n)%flagstruct%external_ic) then if( is_master() ) write(*,*) 'Calling get_external_ic' - call get_external_ic(Atm(n), Atm(n)%domain, .not. do_read_restart) + call get_external_ic(Atm(n), .not. do_read_restart) if( is_master() ) write(*,*) 'IC generated from the specified external source' !4. Restart @@ -358,11 +358,11 @@ subroutine fv_restart(fv_domain, Atm, seconds, days, cold_start, grid_type, this write(*,*) '***** End Note from FV core **************************' write(*,*) ' ' endif - call remap_restart( Atm(n)%domain, Atm(n:n) ) + call remap_restart( Atm(n:n) ) if( is_master() ) write(*,*) 'Done remapping dynamical IC' else if( is_master() ) write(*,*) 'Warm starting, calling fv_io_restart' - call fv_io_read_restart(Atm(n)%domain,Atm(n:n)) + call fv_io_read_restart(Atm(n)%domain_for_read,Atm(n:n)) !====== PJP added DA functionality ====== if (Atm(n)%flagstruct%read_increment) then ! print point in middle of domain for a sanity check