From 205cc762231f6b558c9bc2253debf41aab6cd89b Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 22 Jun 2020 06:06:08 -0600 Subject: [PATCH] Fixes for GNU compilation issues (#22) * Fixes for GNU compilation issues (https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/32) * Fixes for mesh generation in init_grid (https://github.com/NOAA-GFDL/GFDL_atmos_cubed_sphere/pull/39) * Remove trailing whitespace and any tabs * Add default values for nest_*offsets in fv_control * Remove the fix for mesh generation and it will be fixed later * fix GNU out-of-bound error in driver/fvGFS/atmosphere.F90 * Correct Z dimensions in driver/fvGFS/atmosphere.F90 Co-authored-by: Xiaqiong Zhou --- driver/fvGFS/atmosphere.F90 | 4 +- model/boundary.F90 | 8 ++-- model/fv_control.F90 | 9 ++-- model/fv_nesting.F90 | 2 +- model/tp_core.F90 | 4 +- tools/external_ic.F90 | 4 +- tools/fv_diagnostics.F90 | 2 +- tools/fv_restart.F90 | 4 +- tools/test_cases.F90 | 84 ++++++++++++++++++------------------- 9 files changed, 60 insertions(+), 61 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 2a4122c84..d94c5ccc5 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -944,14 +944,14 @@ subroutine atmosphere_hgt (hgt, position, relative, flip) !--- if needed, flip the indexing during this step if (flip) then if (.not. relative) then - z(:,:,1) = Atm(mygrid)%phis(:,:)/grav + z(1:iec-isc+1,1:jec-jsc+1,1) = Atm(mygrid)%phis(isc:iec,jsc:jec)/grav endif do k = 2,npz+1 z(:,:,k) = z(:,:,k-1) - dz(:,:,npz+2-k) enddo else if (.not. relative) then - z(:,:,npz+1) = Atm(mygrid)%phis(:,:)/grav + z(1:iec-isc+1,1:jec-jsc+1,npz+1) = Atm(mygrid)%phis(isc:iec,jsc:jec)/grav endif do k = npz,1,-1 z(:,:,k) = z(:,:,k+1) - dz(:,:,k) diff --git a/model/boundary.F90 b/model/boundary.F90 index 0645ab2e4..5c087c483 100644 --- a/model/boundary.F90 +++ b/model/boundary.F90 @@ -2365,7 +2365,7 @@ subroutine update_coarse_grid_mpp(var_coarse, var_nest, nest_domain, dx, dy, are position = CENTER end if - !Note that *_c does not have values on the parent_proc. + !Note that *_c does not have values on the parent_proc. !Must use isu, etc. to get bounds of update region on parent. call mpp_get_F2C_index(nest_domain, is_c, ie_c, js_c, je_c, is_f, ie_f, js_f, je_f, nest_level=nest_level, position=position) if (child_proc) then @@ -2536,7 +2536,7 @@ subroutine fill_var_coarse(var_coarse, coarse_dat_recv, isd_p, ied_p, jsd_p, jed select case (nestupdate) case (1,2,6,7,8) ! 1 = Conserving update on all variables; 2 = conserving update for cell-centered values; 6 = conserving remap-update -!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) +!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) do k=1,npz do j=jsu,jeu do i=isu,ieu @@ -2559,7 +2559,7 @@ subroutine fill_var_coarse(var_coarse, coarse_dat_recv, isd_p, ied_p, jsd_p, jed select case (nestupdate) case (1,6,7,8) -!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) +!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) do k=1,npz do j=jsu,jeu+1 do i=isu,ieu @@ -2579,7 +2579,7 @@ subroutine fill_var_coarse(var_coarse, coarse_dat_recv, isd_p, ied_p, jsd_p, jed select case (nestupdate) case (1,6,7,8) !averaging update; in-line average for face-averaged values instead of areal average -!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) +!$OMP parallel do default(none) shared(npz,jsu,jeu,isu,ieu,coarse_dat_recv,parent_grid,var_coarse) do k=1,npz do j=jsu,jeu do i=isu,ieu+1 diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 200be4418..1d01c5288 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -143,7 +143,7 @@ module fv_control_mod use fv_mp_mod, only: mp_start, domain_decomp, mp_assign_gid, global_nest_domain use fv_mp_mod, only: broadcast_domains, mp_barrier, is_master, setup_master, grids_master_procs, tile_fine use fv_mp_mod, only: MAX_NNEST, MAX_NTILE - !use test_cases_mod, only: test_case, bubble_do, alpha, nsolitons, soliton_Umax, soliton_size + use test_cases_mod, only: read_namelist_test_case_nml use fv_timing_mod, only: timing_on, timing_off, timing_init, timing_prt use mpp_domains_mod, only: domain2D use mpp_domains_mod, only: mpp_define_nest_domains, nest_domain_type, mpp_get_global_domain @@ -200,7 +200,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) integer, dimension(MAX_NNEST) :: grid_pes = 0 integer, dimension(MAX_NNEST) :: grid_coarse = -1 integer, dimension(MAX_NNEST) :: nest_refine = 3 - integer, dimension(MAX_NNEST) :: nest_ioffsets, nest_joffsets + integer, dimension(MAX_NNEST) :: nest_ioffsets = -999, nest_joffsets = -999 integer, dimension(MAX_NNEST) :: all_npx = 0 integer, dimension(MAX_NNEST) :: all_npy = 0 integer, dimension(MAX_NNEST) :: all_npz = 0 @@ -537,7 +537,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) #endif call read_namelist_fv_grid_nml call read_namelist_fv_core_nml(Atm(this_grid)) ! do options processing here too? - !TODO test_case_nml moved to test_cases + call read_namelist_test_case_nml(Atm(this_grid)%nml_filename) call mpp_get_current_pelist(Atm(this_grid)%pelist, commID=commID) ! for commID call mp_start(commID,halo_update_type) @@ -679,7 +679,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) endif - allocate(Atm(this_grid)%neststruct%child_grids(ngrids)) + allocate(Atm(this_grid)%neststruct%child_grids(ngrids)) do n=1,ngrids Atm(this_grid)%neststruct%child_grids(n) = (grid_coarse(n) == this_grid) allocate(Atm(n)%neststruct%do_remap_bc(ngrids)) @@ -1218,7 +1218,6 @@ subroutine setup_update_regions upoff = Atm(this_grid)%neststruct%upoff do n=2,ngrids - write(*,'(I, A, 4I)') mpp_pe(), 'SETUP_UPDATE_REGIONS 0: ', mpp_pe(), tile_coarse(n), Atm(this_grid)%global_tile if (tile_coarse(n) == Atm(this_grid)%global_tile) then isu = nest_ioffsets(n) diff --git a/model/fv_nesting.F90 b/model/fv_nesting.F90 index 2c4f50efd..01748aff7 100644 --- a/model/fv_nesting.F90 +++ b/model/fv_nesting.F90 @@ -1870,7 +1870,7 @@ end subroutine set_BCs_t0 subroutine d2c_setup(u, v, & ua, va, & - uc, vc, dord4, & + uc, vc, dord4, & isd,ied,jsd,jed, is,ie,js,je, npx,npy, & grid_type, bounded_domain, & se_corner, sw_corner, ne_corner, nw_corner, & diff --git a/model/tp_core.F90 b/model/tp_core.F90 index dac334f83..c5cc3b29e 100644 --- a/model/tp_core.F90 +++ b/model/tp_core.F90 @@ -159,7 +159,7 @@ subroutine fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, & ord_ou = hord if (.not. gridstruct%bounded_domain) & - call copy_corners(q, npx, npy, 2, gridstruct%bounded_domain, bd, & + call copy_corners(q, npx, npy, 2, gridstruct%bounded_domain, bd, & gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) call yppm(fy2, q, cry, ord_in, isd,ied,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dya, gridstruct%bounded_domain, gridstruct%grid_type, lim_fac) @@ -178,7 +178,7 @@ subroutine fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, & call xppm(fx, q_i, crx(is,js), ord_ou, is,ie,isd,ied, js,je,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%bounded_domain, gridstruct%grid_type, lim_fac) if (.not. gridstruct%bounded_domain) & - call copy_corners(q, npx, npy, 1, gridstruct%bounded_domain, bd, & + call copy_corners(q, npx, npy, 1, gridstruct%bounded_domain, bd, & gridstruct%sw_corner, gridstruct%se_corner, gridstruct%nw_corner, gridstruct%ne_corner) call xppm(fx2, q, crx, ord_in, is,ie,isd,ied, jsd,jed,jsd,jed, npx,npy, gridstruct%dxa, gridstruct%bounded_domain, gridstruct%grid_type, lim_fac) diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 0eb166c41..d122f2dd5 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -810,7 +810,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos ) Atm%gridstruct%dxc, Atm%gridstruct%dyc, Atm%gridstruct%sin_sg, & Atm%flagstruct%n_zs_filter, cnst_0p20*Atm%gridstruct%da_min, & .false., oro_g, Atm%gridstruct%bounded_domain, & - Atm%domain, Atm%bd) + Atm%domain, Atm%bd) if ( is_master() ) write(*,*) 'Warning !!! del-2 terrain filter has been applied ', & Atm%flagstruct%n_zs_filter, ' times' else if( Atm%flagstruct%nord_zs_filter == 4 ) then @@ -818,7 +818,7 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos ) Atm%gridstruct%dx, Atm%gridstruct%dy, & Atm%gridstruct%dxc, Atm%gridstruct%dyc, Atm%gridstruct%sin_sg, & Atm%flagstruct%n_zs_filter, .false., oro_g, & - Atm%gridstruct%bounded_domain, & + Atm%gridstruct%bounded_domain, & Atm%domain, Atm%bd) if ( is_master() ) write(*,*) 'Warning !!! del-4 terrain filter has been applied ', & Atm%flagstruct%n_zs_filter, ' times' diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90 index 938b74760..cbdac9c79 100644 --- a/tools/fv_diagnostics.F90 +++ b/tools/fv_diagnostics.F90 @@ -2036,7 +2036,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq) - if( idiag%id_slp>0 .or. idiag%id_tm>0 .or. idiag%id_any_hght>0 .or. idiag%id_hght3d .or. idiag%id_c15>0 .or. idiag%id_ctz ) then + if( idiag%id_slp>0 .or. idiag%id_tm>0 .or. idiag%id_any_hght>0 .or. idiag%id_hght3d>0 .or. idiag%id_c15>0 .or. idiag%id_ctz>0 ) then allocate ( wz(isc:iec,jsc:jec,npz+1) ) call get_height_field(isc, iec, jsc, jec, ngc, npz, Atm(n)%flagstruct%hydrostatic, Atm(n)%delz, & diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index 2496ce828..c4722b33d 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -541,7 +541,7 @@ subroutine fv_restart(fv_domain, Atm, dt_atmos, seconds, days, cold_start, grid_ do n = ntileMe,1,-1 - if (new_nest_topo(n)) then + if (new_nest_topo(n) > 0 ) then call twoway_topo_update(Atm(n), n==this_grid) endif end do @@ -566,7 +566,7 @@ subroutine fv_restart(fv_domain, Atm, dt_atmos, seconds, days, cold_start, grid_ ntdiag = size(Atm(n)%qdiag,4) - if (.not. ideal_test_case(n)) then + if ( ideal_test_case(n) == 0 ) then #ifdef SW_DYNAMICS Atm(n)%pt(:,:,:)=1. #else diff --git a/tools/test_cases.F90 b/tools/test_cases.F90 index 935d4f26d..ae3019a5a 100644 --- a/tools/test_cases.F90 +++ b/tools/test_cases.F90 @@ -230,7 +230,7 @@ module test_cases_mod integer, parameter :: interpOrder = 1 public :: pz0, zz0 - public :: test_case, bubble_do, alpha, tracer_test, wind_field, nsolitons, soliton_Umax, soliton_size + public :: read_namelist_test_case_nml, alpha public :: init_case #ifdef NCDF_OUTPUT public :: output, output_ncdf @@ -6360,26 +6360,26 @@ subroutine init_double_periodic(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, & delp, ak, bk, pt, delz, area, ng, .false., hydrostatic, hybrid_z, domain) - ! *** Add Initial perturbation *** - if (bubble_do) then - r0 = 100.*sqrt(dx_const**2 + dy_const**2) - icenter = npx/2 - jcenter = npy/2 - - do j=js,je - do i=is,ie - dist = (i-icenter)*dx_const*(i-icenter)*dx_const & - +(j-jcenter)*dy_const*(j-jcenter)*dy_const - dist = min(r0, sqrt(dist)) - do k=1,npz - prf = ak(k) + ps(i,j)*bk(k) - if ( prf > 100.E2 ) then - pt(i,j,k) = pt(i,j,k) + 0.01*(1. - (dist/r0)) * prf/ps(i,j) - endif - enddo - enddo - enddo - endif + ! *** Add Initial perturbation *** + if (bubble_do) then + r0 = 100.*sqrt(dx_const**2 + dy_const**2) + icenter = npx/2 + jcenter = npy/2 + + do j=js,je + do i=is,ie + dist = (i-icenter)*dx_const*(i-icenter)*dx_const & + +(j-jcenter)*dy_const*(j-jcenter)*dy_const + dist = min(r0, sqrt(dist)) + do k=1,npz + prf = ak(k) + ps(i,j)*bk(k) + if ( prf > 100.E2 ) then + pt(i,j,k) = pt(i,j,k) + 0.01*(1. - (dist/r0)) * prf/ps(i,j) + endif + enddo + enddo + enddo + endif if ( hydrostatic ) then call p_var(npz, is, ie, js, je, ptop, ptop_min, delp, delz, pt, ps, & pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., & @@ -6734,26 +6734,26 @@ subroutine init_double_periodic(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, .true., hydrostatic, nwat, domain, flagstruct%adiabatic) ! *** Add Initial perturbation *** - if (bubble_do) then - r0 = 10.e3 - zc = 1.4e3 ! center of bubble from surface - icenter = (npx-1)/2 + 1 - jcenter = (npy-1)/2 + 1 - do k=1, npz - zm = 0.5*(ze1(k)+ze1(k+1)) - ptmp = ( (zm-zc)/zc ) **2 - if ( ptmp < 1. ) then - do j=js,je - do i=is,ie - dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2 - if ( dist < 1. ) then - pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist)) - endif - enddo - enddo - endif - enddo - endif + if (bubble_do) then + r0 = 10.e3 + zc = 1.4e3 ! center of bubble from surface + icenter = (npx-1)/2 + 1 + jcenter = (npy-1)/2 + 1 + do k=1, npz + zm = 0.5*(ze1(k)+ze1(k+1)) + ptmp = ( (zm-zc)/zc ) **2 + if ( ptmp < 1. ) then + do j=js,je + do i=is,ie + dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2 + if ( dist < 1. ) then + pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist)) + endif + enddo + enddo + endif + enddo + endif case ( 101 ) @@ -6904,6 +6904,7 @@ subroutine read_namelist_test_case_nml(nml_filename) integer :: ierr, f_unit, unit, ios #include + namelist /test_case_nml/test_case, bubble_do, alpha, nsolitons,soliton_Umax, soliton_size unit = stdlog() @@ -6911,7 +6912,6 @@ subroutine read_namelist_test_case_nml(nml_filename) alpha = 0. bubble_do = .false. test_case = 11 ! (USGS terrain) - namelist /test_case_nml/test_case, bubble_do, alpha, nsolitons, soliton_Umax, soliton_size #ifdef INTERNAL_FILE_NML ! Read Test_Case namelist