diff --git a/.gitmodules b/.gitmodules
index 6bb663df1..047dc19f3 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -8,8 +8,8 @@
branch = main
[submodule "ccpp/physics"]
path = ccpp/physics
- url = https://github.com/NCAR/ccpp-physics
- branch = main
+ url = https://github.com/dustinswales/ccpp-physics
+ branch = ufs-dev-PR30
[submodule "upp"]
path = upp
url = https://github.com/NOAA-EMC/UPP
diff --git a/atmos_cubed_sphere b/atmos_cubed_sphere
index a839395d5..3f7b57a28 160000
--- a/atmos_cubed_sphere
+++ b/atmos_cubed_sphere
@@ -1 +1 @@
-Subproject commit a839395d542ab860ae0afa7601e35b37d68d773e
+Subproject commit 3f7b57a28b334c7f2485b26c29f1d5b9f3c1cd7d
diff --git a/atmos_model.F90 b/atmos_model.F90
index 81a5272c3..2a5091bd6 100644
--- a/atmos_model.F90
+++ b/atmos_model.F90
@@ -96,8 +96,7 @@ module atmos_model_mod
DIAG_SIZE
use fv_iau_mod, only: iau_external_data_type,getiauforcing,iau_initialize
use module_fv3_config, only: output_1st_tstep_rst, first_kdt, nsout, &
- restart_endfcst, output_fh, fcst_mpi_comm, &
- fcst_ntasks
+ output_fh, fcst_mpi_comm, fcst_ntasks
use module_block_data, only: block_atmos_copy, block_data_copy, &
block_data_copy_or_fill, &
block_data_combine_fractions
@@ -178,6 +177,7 @@ module atmos_model_mod
logical :: debug = .false.
!logical :: debug = .true.
logical :: sync = .false.
+logical :: restart_endfcst = .false.
real :: avg_max_length=3600.
logical :: ignore_rst_cksum = .false.
namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, ccpp_suite, avg_max_length, &
@@ -1477,6 +1477,10 @@ subroutine update_atmos_chemistry(state, rc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__, rcToReturn=rc)) return
+ call cplFieldGet(state,'inst_pres_interface', farrayPtr3d=prsi, rc=localrc)
+ if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, file=__FILE__, rcToReturn=rc)) return
+
if (GFS_Control%cplaqm) then
call cplFieldGet(state,'canopy_moisture_storage', farrayPtr2d=canopy, rc=localrc)
@@ -1541,10 +1545,6 @@ subroutine update_atmos_chemistry(state, rc)
else
- call cplFieldGet(state,'inst_pres_interface', farrayPtr3d=prsi, rc=localrc)
- if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
- line=__LINE__, file=__FILE__, rcToReturn=rc)) return
-
call cplFieldGet(state,'inst_liq_nonconv_tendency_levels', &
farrayPtr3d=pflls, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -1592,6 +1592,7 @@ subroutine update_atmos_chemistry(state, rc)
ix = Atm_block%ixp(ib,jb)
!--- interface values
phii(i,j,k) = GFS_data(nb)%Statein%phii(ix,k)
+ prsi(i,j,k) = GFS_data(nb)%Statein%prsi(ix,k)
!--- layer values
prsl(i,j,k) = GFS_Data(nb)%Statein%prsl(ix,k)
phil(i,j,k) = GFS_Data(nb)%Statein%phil(ix,k)
@@ -1600,8 +1601,6 @@ subroutine update_atmos_chemistry(state, rc)
va (i,j,k) = GFS_Data(nb)%Stateout%gv0(ix,k)
cldfra(i,j,k) = GFS_Data(nb)%IntDiag%cldfra(ix,k)
if (.not.GFS_Control%cplaqm) then
- !--- interface values
- prsi(i,j,k) = GFS_data(nb)%Statein%prsi(ix,k)
!--- layer values
pfils (i,j,k) = GFS_Data(nb)%Coupling%pfi_lsan(ix,k)
pflls (i,j,k) = GFS_Data(nb)%Coupling%pfl_lsan(ix,k)
@@ -1620,8 +1619,7 @@ subroutine update_atmos_chemistry(state, rc)
nb = Atm_block%blkno(ib,jb)
ix = Atm_block%ixp(ib,jb)
phii(i,j,k) = GFS_data(nb)%Statein%phii(ix,k)
- if (.not.GFS_Control%cplaqm) &
- prsi(i,j,k) = GFS_data(nb)%Statein%prsi(ix,k)
+ prsi(i,j,k) = GFS_data(nb)%Statein%prsi(ix,k)
enddo
enddo
@@ -1728,6 +1726,7 @@ subroutine update_atmos_chemistry(state, rc)
if (GFS_control%debug) then
! -- diagnostics
+ write(6,'("update_atmos: prsi - min/max/avg",3g16.6)') minval(prsi), maxval(prsi), sum(prsi)/size(prsi)
write(6,'("update_atmos: phii - min/max/avg",3g16.6)') minval(phii), maxval(phii), sum(phii)/size(phii)
write(6,'("update_atmos: prsl - min/max/avg",3g16.6)') minval(prsl), maxval(prsl), sum(prsl)/size(prsl)
write(6,'("update_atmos: phil - min/max/avg",3g16.6)') minval(phil), maxval(phil), sum(phil)/size(phil)
@@ -1766,7 +1765,6 @@ subroutine update_atmos_chemistry(state, rc)
write(6,'("update_atmos: xlai - min/max/avg",3g16.6)') minval(xlai), maxval(xlai), sum(xlai)/size(xlai)
write(6,'("update_atmos: stype - min/max/avg",3g16.6)') minval(stype), maxval(stype), sum(stype)/size(stype)
else
- write(6,'("update_atmos: prsi - min/max/avg",3g16.6)') minval(prsi), maxval(prsi), sum(prsi)/size(prsi)
write(6,'("update_atmos: flake - min/max/avg",3g16.6)') minval(flake), maxval(flake), sum(flake)/size(flake)
write(6,'("update_atmos: focn - min/max/avg",3g16.6)') minval(focn), maxval(focn), sum(focn)/size(focn)
write(6,'("update_atmos: shfsfc - min/max/avg",3g16.6)') minval(shfsfc), maxval(shfsfc), sum(shfsfc)/size(shfsfc)
diff --git a/ccpp/data/CCPP_typedefs.F90 b/ccpp/data/CCPP_typedefs.F90
index 1c6e4316f..b886beb3e 100644
--- a/ccpp/data/CCPP_typedefs.F90
+++ b/ccpp/data/CCPP_typedefs.F90
@@ -249,7 +249,7 @@ module CCPP_typedefs
real (kind=kind_phys), pointer :: qss_ice(:) => null() !<
real (kind=kind_phys), pointer :: qss_land(:) => null() !<
real (kind=kind_phys), pointer :: qss_water(:) => null() !<
- logical :: radar_reset !<
+ logical :: fullradar_diag !<
real (kind=kind_phys) :: raddt !<
real (kind=kind_phys), pointer :: rainmp(:) => null() !<
real (kind=kind_phys), pointer :: raincd(:) => null() !<
@@ -1478,11 +1478,11 @@ subroutine gfs_interstitial_phys_reset (Interstitial, Model)
! Use same logic in UFS to reset Thompson extended diagnostics
Interstitial%ext_diag_thompson_reset = Interstitial%max_hourly_reset
!
- ! Set flag for resetting radar reflectivity calculation
- if (Model%nsradar_reset<0) then
- Interstitial%radar_reset = .true.
+ ! Frequency flag for computing the full radar reflectivity (water coated ice)
+ if (Model%nsfullradar_diag<0) then
+ Interstitial%fullradar_diag = .true.
else
- Interstitial%radar_reset = mod(Model%kdt-1, nint(Model%nsradar_reset/Model%dtp)) == 0
+ Interstitial%fullradar_diag = (Model%kdt == 1 .or. mod(Model%kdt, nint(Model%nsfullradar_diag/Model%dtp)) == 0)
end if
!
end subroutine gfs_interstitial_phys_reset
diff --git a/ccpp/data/CCPP_typedefs.meta b/ccpp/data/CCPP_typedefs.meta
index c5e7fb5be..c7c5b4521 100644
--- a/ccpp/data/CCPP_typedefs.meta
+++ b/ccpp/data/CCPP_typedefs.meta
@@ -1699,9 +1699,9 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
-[radar_reset]
- standard_name = flag_for_resetting_radar_reflectivity_calculation
- long_name = flag for resetting radar reflectivity calculation
+[fullradar_diag]
+ standard_name = do_full_radar_reflectivity
+ long_name = flag for computing full radar reflectivity
units = flag
dimensions = ()
type = logical
diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90
index 38df00ddc..0dac224ff 100644
--- a/ccpp/data/GFS_typedefs.F90
+++ b/ccpp/data/GFS_typedefs.F90
@@ -372,7 +372,6 @@ module GFS_typedefs
real (kind=kind_phys), pointer :: clw_surf_ice(:) => null() !< RUC LSM: moist cloud water mixing ratio at surface over ice
real (kind=kind_phys), pointer :: qwv_surf_land(:) => null() !< RUC LSM: water vapor mixing ratio at surface over land
real (kind=kind_phys), pointer :: qwv_surf_ice(:) => null() !< RUC LSM: water vapor mixing ratio at surface over ice
- real (kind=kind_phys), pointer :: rhofr(:) => null() !< RUC LSM: density of frozen precipitation
real (kind=kind_phys), pointer :: tsnow_land(:) => null() !< RUC LSM: snow temperature at the bottom of the first snow layer over land
real (kind=kind_phys), pointer :: tsnow_ice(:) => null() !< RUC LSM: snow temperature at the bottom of the first snow layer over ice
real (kind=kind_phys), pointer :: snowfallac_land(:) => null() !< ruc lsm diagnostics over land
@@ -920,7 +919,7 @@ module GFS_typedefs
logical :: ltaerosol !< flag for aerosol version
logical :: mraerosol !< flag for merra2_aerosol_aware
logical :: lradar !< flag for radar reflectivity
- real(kind=kind_phys) :: nsradar_reset !< seconds between resetting radar reflectivity calculation
+ real(kind=kind_phys) :: nsfullradar_diag!< seconds between resetting radar reflectivity calculation
real(kind=kind_phys) :: ttendlim !< temperature tendency limiter per time step in K/s
logical :: ext_diag_thompson !< flag for extended diagnostic output from Thompson
integer :: thompson_ext_ndiag3d=37 !< number of 3d arrays for extended diagnostic output from Thompson
@@ -954,6 +953,7 @@ module GFS_typedefs
integer :: lsnow_lsm !< maximum number of snow layers internal to land surface model
integer :: lsnow_lsm_lbound!< lower bound for snow arrays, depending on lsnow_lsm
integer :: lsnow_lsm_ubound!< upper bound for snow arrays, depending on lsnow_lsm
+ logical :: exticeden !< flag for calculating frozen precip ice density outside of the LSM
real(kind=kind_phys), pointer :: zs(:) => null() !< depth of soil levels for land surface model
real(kind=kind_phys), pointer :: dzs(:) => null() !< thickness of soil levels for land surface model
real(kind=kind_phys), pointer :: pores(:) => null() !< max soil moisture for a given soil type for land surface model
@@ -1758,6 +1758,14 @@ module GFS_typedefs
real (kind=kind_phys), pointer :: toticeb(:) => null() !< accumulated ice precipitation in bucket (kg/m2)
real (kind=kind_phys), pointer :: totsnwb(:) => null() !< accumulated snow precipitation in bucket (kg/m2)
real (kind=kind_phys), pointer :: totgrpb(:) => null() !< accumulated graupel precipitation in bucket (kg/m2)
+ real (kind=kind_phys), pointer :: frzr (:) => null() !< accumulated surface freezing rain (m)
+ real (kind=kind_phys), pointer :: frzrb (:) => null() !< accumulated surface freezing rain in bucket (m)
+ real (kind=kind_phys), pointer :: frozr (:) => null() !< accumulated surface graupel (m)
+ real (kind=kind_phys), pointer :: frozrb (:) => null() !< accumulated surface graupel in bucket (m)
+ real (kind=kind_phys), pointer :: tsnowp (:) => null() !< accumulated surface snowfall (m)
+ real (kind=kind_phys), pointer :: tsnowpb(:) => null() !< accumulated surface snowfall in bucket (m)
+ real (kind=kind_phys), pointer :: rhonewsn1(:) => null() !< precipitation ice density outside RUC LSM (kg/m3)
+ real (kind=kind_phys), pointer :: rhosnf(:) => null() !< precipitation ice density inside RUC LSM (kg/m3)
!--- MYNN variables
real (kind=kind_phys), pointer :: edmf_a (:,:) => null() !
@@ -2449,7 +2457,6 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
allocate (Sfcprop%clw_surf_ice (IM))
allocate (Sfcprop%qwv_surf_land (IM))
allocate (Sfcprop%qwv_surf_ice (IM))
- allocate (Sfcprop%rhofr (IM))
allocate (Sfcprop%tsnow_land (IM))
allocate (Sfcprop%tsnow_ice (IM))
allocate (Sfcprop%snowfallac_land (IM))
@@ -2465,7 +2472,6 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
Sfcprop%qwv_surf_land = clear_val
Sfcprop%qwv_surf_ice = clear_val
Sfcprop%flag_frsoil = clear_val
- Sfcprop%rhofr = clear_val
Sfcprop%tsnow_land = clear_val
Sfcprop%tsnow_ice = clear_val
Sfcprop%snowfallac_land = clear_val
@@ -3164,7 +3170,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
logical :: ltaerosol = .false. !< flag for aerosol version
logical :: mraerosol = .false. !< flag for merra2_aerosol_aware
logical :: lradar = .false. !< flag for radar reflectivity
- real(kind=kind_phys) :: nsradar_reset = -999.0 !< seconds between resetting radar reflectivity calculation, set to <0 for every time step
+ real(kind=kind_phys) :: nsfullradar_diag = -999.0 !< seconds between resetting radar reflectivity calculation, set to <0 for every time step
real(kind=kind_phys) :: ttendlim = -999.0 !< temperature tendency limiter, set to <0 to deactivate
logical :: ext_diag_thompson = .false. !< flag for extended diagnostic output from Thompson
real(kind=kind_phys) :: dt_inner = -999.0 !< time step for the inner loop
@@ -3182,6 +3188,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
integer :: lsoil = 4 !< number of soil layers
integer :: lsoil_lsm = -1 !< number of soil layers internal to land surface model; -1 use lsoil
integer :: lsnow_lsm = 3 !< maximum number of snow layers internal to land surface model
+ logical :: exticeden = .false. !< Use variable precip ice density for NOAH LSM if true or original formulation
logical :: rdlai = .false. !< read LAI from input file (for RUC LSM or NOAH LSM WRFv4)
logical :: ua_phys = .false. !< flag for using University of Arizona? extension to NOAH LSM WRFv4
logical :: usemonalb = .true. !< flag to read surface diffused shortwave albedo from input file for NOAH LSM WRFv4
@@ -3576,7 +3583,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
mg_do_graupel, mg_do_hail, mg_nccons, mg_nicons, mg_ngcons, &
mg_ncnst, mg_ninst, mg_ngnst, sed_supersat, do_sb_physics, &
mg_alf, mg_qcmin, mg_do_ice_gmao, mg_do_liq_liu, &
- ltaerosol, lradar, nsradar_reset, lrefres, ttendlim, &
+ ltaerosol, lradar, nsfullradar_diag, lrefres, ttendlim, &
ext_diag_thompson, dt_inner, lgfdlmprad, &
sedi_semi, decfl, &
nssl_cccn, nssl_alphah, nssl_alphahl, &
@@ -3586,7 +3593,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
!--- land/surface model control
lsm, lsoil, lsoil_lsm, lsnow_lsm, kice, rdlai, &
nmtvr, ivegsrc, use_ufo, iopt_thcnd, ua_phys, usemonalb, &
- aoasis, fasdas, &
+ aoasis, fasdas,exticeden, &
! Noah MP options
iopt_dveg,iopt_crs,iopt_btr,iopt_run,iopt_sfc, iopt_frz, &
iopt_inf, iopt_rad,iopt_alb,iopt_snf,iopt_tbot,iopt_stc, &
@@ -4166,7 +4173,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
stop
end if
Model%lradar = lradar
- Model%nsradar_reset = nsradar_reset
+ Model%nsfullradar_diag = nsfullradar_diag
Model%ttendlim = ttendlim
Model%ext_diag_thompson= ext_diag_thompson
if (dt_inner>0) then
@@ -4274,7 +4281,15 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%ivegsrc = ivegsrc
Model%isot = isot
Model%use_ufo = use_ufo
-
+ Model%exticeden = exticeden
+ if (Model%exticeden .and. &
+ (Model%imp_physics /= Model%imp_physics_gfdl .and. Model%imp_physics /= Model%imp_physics_thompson .and. &
+ Model%imp_physics /= Model%imp_physics_nssl )) then
+ !see GFS_MP_generic_post.F90; exticeden is only compatible with GFDL,
+ !Thompson, or NSSL MP
+ print *,' Using exticeden = T is only valid when using GFDL, Thompson, or NSSL microphysics.'
+ stop
+ end if
! GFDL surface layer options
Model%lcurr_sf = lcurr_sf
Model%pert_cd = pert_cd
@@ -4295,7 +4310,11 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%iopt_inf = iopt_inf
Model%iopt_rad = iopt_rad
Model%iopt_alb = iopt_alb
- Model%iopt_snf = iopt_snf
+ if (Model%lsm==Model%lsm_noahmp .and. Model%exticeden .and. iopt_snf == 4) then
+ Model%iopt_snf = 5
+ else
+ Model%iopt_snf = iopt_snf
+ end if
Model%iopt_tbot = iopt_tbot
Model%iopt_stc = iopt_stc
Model%iopt_trs = iopt_trs
@@ -5175,6 +5194,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
if (Model%me == Model%master) then
if (Model%lsm == 1) then
print *,' NOAH Land Surface Model used'
+
elseif (Model%lsm == 0) then
print *,' OSU no longer supported - job aborted'
stop
@@ -5471,7 +5491,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
' decfl=',decfl, &
' effr_in =',Model%effr_in, &
' lradar =',Model%lradar, &
- ' nsradar_reset =',Model%nsradar_reset, &
+ ' nsfullradar_diag =',Model%nsfullradar_diag, &
' num_p3d =',Model%num_p3d, &
' num_p2d =',Model%num_p2d
@@ -5985,7 +6005,7 @@ subroutine control_print(Model)
print *, ' ltaerosol : ', Model%ltaerosol
print *, ' mraerosol : ', Model%mraerosol
print *, ' lradar : ', Model%lradar
- print *, ' nsradar_reset : ', Model%nsradar_reset
+ print *, ' nsfullradar_diag : ', Model%nsfullradar_diag
print *, ' lrefres : ', Model%lrefres
print *, ' ttendlim : ', Model%ttendlim
print *, ' ext_diag_thompson : ', Model%ext_diag_thompson
@@ -7009,6 +7029,14 @@ subroutine diag_create (Diag, IM, Model)
allocate (Diag%epi (IM))
allocate (Diag%smcwlt2 (IM))
allocate (Diag%smcref2 (IM))
+ allocate (Diag%rhonewsn1 (IM))
+ allocate (Diag%rhosnf (IM))
+ allocate (Diag%frzr (IM))
+ allocate (Diag%frzrb (IM))
+ allocate (Diag%frozr (IM))
+ allocate (Diag%frozrb (IM))
+ allocate (Diag%tsnowp (IM))
+ allocate (Diag%tsnowpb (IM))
if (.not. Model%lsm == Model%lsm_ruc) then
allocate (Diag%wet1 (IM))
end if
@@ -7344,6 +7372,9 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
Diag%toticeb = zero
Diag%totsnwb = zero
Diag%totgrpb = zero
+ Diag%frzrb = zero
+ Diag%frozrb = zero
+ Diag%tsnowpb = zero
!--- MYNN variables:
if (Model%do_mynnedmf) then
@@ -7464,7 +7495,9 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
Diag%t02min = 999.
Diag%rh02max = -999.
Diag%rh02min = 999.
- Diag%pratemax = 0.
+ Diag%pratemax = 0.
+ Diag%rhonewsn1 = 200.
+ Diag%rhosnf = -1.e3
set_totprcp = .false.
if (present(linit) ) set_totprcp = linit
if (present(iauwindow_center) ) set_totprcp = iauwindow_center
@@ -7474,6 +7507,9 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
Diag%totice = zero
Diag%totsnw = zero
Diag%totgrp = zero
+ Diag%frzr = zero
+ Diag%frozr = zero
+ Diag%tsnowp = zero
endif
end subroutine diag_phys_zero
diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta
index 439d37cc8..edbb1dbaa 100644
--- a/ccpp/data/GFS_typedefs.meta
+++ b/ccpp/data/GFS_typedefs.meta
@@ -1627,14 +1627,6 @@
type = real
kind = kind_phys
active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme)
-[rhofr]
- standard_name = frozen_precipitation_density
- long_name = density of frozen precipitation
- units = kg m-3
- dimensions = (horizontal_loop_extent)
- type = real
- kind = kind_phys
- active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme)
[tsnow_land]
standard_name = temperature_in_surface_snow_at_surface_adjacent_layer_over_land
long_name = snow temperature at the bottom of the first snow layer over land
@@ -4279,6 +4271,12 @@
units = count
dimensions = ()
type = integer
+[exticeden]
+ standard_name = do_external_surface_frozen_precipitation_density
+ long_name = flag for calculating frozen precip ice density outside of the LSM
+ units = flag
+ dimensions = ()
+ type = logical
[zs]
standard_name = depth_of_soil_layers
long_name = depth of soil levels for land surface model
@@ -7614,6 +7612,62 @@
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
+[frzr]
+ standard_name = cumulative_lwe_thickness_of_surface_freezing_rain_amount
+ long_name = accumulated surface freezing rain
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[frzrb]
+ standard_name = cumulative_lwe_thickness_of_surface_freezing_rain_amount_in_bucket
+ long_name = accumulated surface freezing rain in bucket
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[frozr]
+ standard_name = cumulative_lwe_thickness_of_surface_graupel_amount
+ long_name = accumulated surface graupel
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[frozrb]
+ standard_name = cumulative_lwe_thickness_of_surface_graupel_amount_in_bucket
+ long_name = accumulated surface graupel in bucket
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[tsnowp]
+ standard_name = cumulative_lwe_thickness_of_surface_snow_amount
+ long_name = accumulated surface snow
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[tsnowpb]
+ standard_name = cumulative_lwe_thickness_of_surface_snow_amount_in_bucket
+ long_name = accumulated surface snow in bucket
+ units = m
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[rhonewsn1]
+ standard_name = surface_frozen_precipitation_density
+ long_name = density of precipitation ice
+ units = kg m-3
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
+[rhosnf]
+ standard_name = lsm_internal_surface_frozen_precipitation_density
+ long_name = density of frozen precipitation
+ units = kg m-3
+ dimensions = (horizontal_loop_extent)
+ type = real
+ kind = kind_phys
[train]
standard_name = accumulated_change_of_air_temperature_due_to_FA_scheme
long_name = accumulated change of air temperature due to FA MP scheme
diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90
index ff01c6630..9617b3fc2 100644
--- a/ccpp/driver/GFS_diagnostics.F90
+++ b/ccpp/driver/GFS_diagnostics.F90
@@ -1642,6 +1642,89 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%pratemax(:)
enddo
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'frzr'
+ ExtDiag(idx)%desc = 'accumulated surface freezing rain'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%frzr(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'frzrb'
+ ExtDiag(idx)%desc = 'accumulated surface freezing rain in bucket'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%frzrb(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'frozr'
+ ExtDiag(idx)%desc = 'accumulated surface graupel'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%frozr(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'frozrb'
+ ExtDiag(idx)%desc = 'accumulated surface graupel in bucket'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%frozrb(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'tsnowp'
+ ExtDiag(idx)%desc = 'accumulated surface snow'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%tsnowp(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'tsnowpb'
+ ExtDiag(idx)%desc = 'accumulated surface snow in bucket'
+ ExtDiag(idx)%unit = 'kg/m**2'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ ExtDiag(idx)%cnvfac = cn_th
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%tsnowpb(:)
+ enddo
+
+ idx = idx + 1
+ ExtDiag(idx)%axes = 2
+ ExtDiag(idx)%name = 'rhonewsn'
+ ExtDiag(idx)%desc = 'precipitation ice density'
+ ExtDiag(idx)%unit = 'kg m^-3'
+ ExtDiag(idx)%mod_name = 'gfs_phys'
+ allocate (ExtDiag(idx)%data(nblks))
+ do nb = 1,nblks
+ ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%rhonewsn1(:)
+ enddo
+
idx = idx + 1
ExtDiag(idx)%axes = 2
ExtDiag(idx)%name = 'rain'
diff --git a/ccpp/driver/GFS_restart.F90 b/ccpp/driver/GFS_restart.F90
index 4774ff299..c3930d752 100644
--- a/ccpp/driver/GFS_restart.F90
+++ b/ccpp/driver/GFS_restart.F90
@@ -88,6 +88,15 @@ subroutine GFS_restart_populate (Restart, Model, Statein, Stateout, Sfcprop, &
else if( trim(ExtDiag(idx)%name) == 'totgrp_ave') then
ndiag_rst = ndiag_rst +1
ndiag_idx(ndiag_rst) = idx
+ else if( trim(ExtDiag(idx)%name) == 'tsnowp') then
+ ndiag_rst = ndiag_rst +1
+ ndiag_idx(ndiag_rst) = idx
+ else if( trim(ExtDiag(idx)%name) == 'frozr') then
+ ndiag_rst = ndiag_rst +1
+ ndiag_idx(ndiag_rst) = idx
+ else if( trim(ExtDiag(idx)%name) == 'frzr') then
+ ndiag_rst = ndiag_rst +1
+ ndiag_idx(ndiag_rst) = idx
endif
endif
enddo
diff --git a/ccpp/physics b/ccpp/physics
index fb0a90fc5..cb84c3eb0 160000
--- a/ccpp/physics
+++ b/ccpp/physics
@@ -1 +1 @@
-Subproject commit fb0a90fc54cc9ce9e31e085953ed9971c4282e25
+Subproject commit cb84c3eb026f939580fecbc365d4975bef3b828e
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml
index 1760afe97..3dc3c8d54 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml
@@ -62,8 +62,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml
index 90693645c..254df77e0 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml
@@ -60,8 +60,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml
index 0e7127bb4..746c8ceb7 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml
@@ -57,8 +57,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml
index 8d3bab273..2022ef4d6 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml
@@ -57,8 +57,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml
index b61422021..f66543c80 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml
@@ -55,8 +55,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml
index 0cb393317..665a0b7e2 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml
@@ -55,8 +55,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml
index e68c9e687..a5db1110b 100644
--- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml
+++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml
@@ -57,8 +57,8 @@
satmedmfvdifq
GFS_PBL_generic_post
GFS_GWD_generic_pre
- cires_ugwp
- cires_ugwp_post
+ unified_ugwp
+ unified_ugwp_post
GFS_GWD_generic_post
GFS_suite_stateout_update
ozphys_2015
diff --git a/io/FV3GFS_io.F90 b/io/FV3GFS_io.F90
index 775d7630d..6ee558fcb 100644
--- a/io/FV3GFS_io.F90
+++ b/io/FV3GFS_io.F90
@@ -3784,7 +3784,7 @@ subroutine add_field_to_phybundle(var_name,long_name,units,cell_methods, axes,ph
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
- name='missing_value',value=missing_value,rc=rc)
+ name='missing_value',value=real(missing_value,kind=4),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
@@ -3794,7 +3794,7 @@ subroutine add_field_to_phybundle(var_name,long_name,units,cell_methods, axes,ph
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
- name='_FillValue',value=missing_value,rc=rc)
+ name='_FillValue',value=real(missing_value,kind=4),rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90
index ff142b559..c3b76701a 100644
--- a/io/module_wrt_grid_comp.F90
+++ b/io/module_wrt_grid_comp.F90
@@ -754,12 +754,14 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock,
allocate( wrt_int_state%out_grid_info(n)%latPtr(wrt_int_state%out_grid_info(n)%i_start:wrt_int_state%out_grid_info(n)%i_end, &
wrt_int_state%out_grid_info(n)%j_start:wrt_int_state%out_grid_info(n)%j_end) )
- do j=wrt_int_state%out_grid_info(n)%j_start, wrt_int_state%out_grid_info(n)%j_end
- do i=wrt_int_state%out_grid_info(n)%i_start, wrt_int_state%out_grid_info(n)%i_end
- wrt_int_state%out_grid_info(n)%latPtr(i,j) = latPtr(i,j)
- wrt_int_state%out_grid_info(n)%lonPtr(i,j) = lonPtr(i,j)
- enddo
- enddo
+ if ( trim(output_grid(n)) /= 'regional_latlon_moving' .and. trim(output_grid(n)) /= 'rotated_latlon_moving' ) then
+ do j=wrt_int_state%out_grid_info(n)%j_start, wrt_int_state%out_grid_info(n)%j_end
+ do i=wrt_int_state%out_grid_info(n)%i_start, wrt_int_state%out_grid_info(n)%i_end
+ wrt_int_state%out_grid_info(n)%latPtr(i,j) = latPtr(i,j)
+ wrt_int_state%out_grid_info(n)%lonPtr(i,j) = lonPtr(i,j)
+ enddo
+ enddo
+ endif
wrt_int_state%out_grid_info(n)%im = imo(n)
wrt_int_state%out_grid_info(n)%jm = jmo(n)
@@ -1076,19 +1078,20 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock,
name="grid", value="latlon", rc=rc)
call ESMF_AttributeAdd(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
attrList=(/"lon1","lat1","lon2","lat2","dlon","dlat"/), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lon1", value=lon1(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lat1", value=lat1(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lon2", value=lon2(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lat2", value=lat2(grid_id), rc=rc)
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
name="dlon", value=dlon(grid_id), rc=rc)
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
name="dlat", value=dlat(grid_id), rc=rc)
-
+ if (trim(output_grid(grid_id)) /= 'regional_latlon_moving') then
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lon1", value=lon1(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lat1", value=lat1(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lon2", value=lon2(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lat2", value=lat2(grid_id), rc=rc)
+ endif
else if (trim(output_grid(grid_id)) == 'rotated_latlon' &
.or. trim(output_grid(grid_id)) == 'rotated_latlon_moving') then
@@ -1108,19 +1111,20 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock,
name="cen_lon", value=cen_lon(grid_id), rc=rc)
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
name="cen_lat", value=cen_lat(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lon1", value=lon1(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lat1", value=lat1(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lon2", value=lon2(grid_id), rc=rc)
- call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
- name="lat2", value=lat2(grid_id), rc=rc)
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
name="dlon", value=dlon(grid_id), rc=rc)
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
name="dlat", value=dlat(grid_id), rc=rc)
-
+ if (trim(output_grid(grid_id)) /= 'rotated_latlon_moving') then
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lon1", value=lon1(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lat1", value=lat1(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lon2", value=lon2(grid_id), rc=rc)
+ call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
+ name="lat2", value=lat2(grid_id), rc=rc)
+ endif
else if (trim(output_grid(grid_id)) == 'lambert_conformal') then
call ESMF_AttributeSet(wrt_int_state%wrtFB(i), convention="NetCDF", purpose="FV3", &
@@ -1393,7 +1397,7 @@ subroutine wrt_initialize_p1(wrt_comp, imp_state_write, exp_state_write, clock,
deallocate(attNameList, attNameList2, typekindList)
!
-! write_init_tim = MPI_Wtime() - btim0
+! write_init_tim = MPI_Wtime() - btim0
!
!-----------------------------------------------------------------------
!
@@ -1931,22 +1935,31 @@ subroutine wrt_run(wrt_comp, imp_state_write, exp_state_write,clock,rc)
#ifdef INLINE_POST
wbeg = MPI_Wtime()
do n=1,ngrids
- if (trim(output_grid(n)) == 'regional_latlon' .or. &
- trim(output_grid(n)) == 'regional_latlon_moving' .or. &
- trim(output_grid(n)) == 'rotated_latlon' .or. &
- trim(output_grid(n)) == 'rotated_latlon_moving' .or. &
- trim(output_grid(n)) == 'lambert_conformal') then
-
- !mask fields according to sfc pressure
- do nbdl=1, wrt_int_state%FBCount
- call mask_fields(wrt_int_state%wrtFB(nbdl),rc)
- if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
- enddo
- lmask_fields = .true.
- endif
- call post_run_fv3(wrt_int_state, n, mype, wrt_mpi_comm, lead_write_task, &
- itasks, jtasks, nf_hours, nf_minutes, nf_seconds)
+ if (trim(output_grid(n)) /= 'cubed_sphere_grid') then
+
+ if (trim(output_grid(n)) == 'regional_latlon' .or. &
+ trim(output_grid(n)) == 'regional_latlon_moving' .or. &
+ trim(output_grid(n)) == 'rotated_latlon' .or. &
+ trim(output_grid(n)) == 'rotated_latlon_moving' .or. &
+ trim(output_grid(n)) == 'lambert_conformal') then
+
+ !mask fields according to sfc pressure
+ do nbdl=1, wrt_int_state%FBCount
+ call mask_fields(wrt_int_state%wrtFB(nbdl),rc)
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+ enddo
+ lmask_fields = .true.
+ endif
+
+ call post_run_fv3(wrt_int_state, n, mype, wrt_mpi_comm, lead_write_task, &
+ itasks, jtasks, nf_hours, nf_minutes, nf_seconds)
+ else
+ rc = ESMF_RC_NOT_IMPL
+ print *,'Inline post not available for cubed_sphere_grid'
+ if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+ line=__LINE__, file=__FILE__)) return
+ endif
enddo
wend = MPI_Wtime()
if (lprnt) then
diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 6b08c4785..58852ce41 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -527,7 +527,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
rel_vort_max, rel_vort_maxhy1, refd_max, &
refdm10c_max, u10max, v10max, wspd10max, sfcuxi, &
sfcvxi, t10m, t10avg, psfcavg, akhsavg, akmsavg, &
- albedo, tg, prate_max, pwat
+ albedo, tg, prate_max, pwat, snow_acm, snow_bkt, &
+ acgraup, graup_bucket, acfrain, frzrn_bucket
use soil, only: sldpth, sh2o, smc, stc
use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
use ctlblk_mod, only: im, jm, lm, lp1, jsta, jend, jsta_2l, jend_2u, jsta_m,jend_m, &
@@ -1191,6 +1192,72 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
enddo
endif
+ !Accumulated snowfall
+ if(trim(fieldname)=='tsnowp') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snow_acm,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ snow_acm(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) snow_acm(i,j) = spval
+ enddo
+ enddo
+ endif
+
+ !Snowfall bucket
+ if(trim(fieldname)=='tsnowpb') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,snow_bkt,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ snow_bkt(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) snow_bkt(i,j) = spval
+ enddo
+ enddo
+ endif
+
+ !Accumulated graupel
+ if(trim(fieldname)=='frozr') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,acgraup,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ acgraup(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) acgraup(i,j) = spval
+ enddo
+ enddo
+ endif
+
+ !Graupel bucket
+ if(trim(fieldname)=='frozrb') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,graup_bucket,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ graup_bucket(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) graup_bucket(i,j) = spval
+ enddo
+ enddo
+ endif
+
+ !Accumulated freezing rain
+ if(trim(fieldname)=='frzr') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,acfrain,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ acfrain(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) acfrain(i,j) = spval
+ enddo
+ enddo
+ endif
+
+ !Freezing rain bucket
+ if(trim(fieldname)=='frzrb') then
+ !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,frzrn_bucket,arrayr42d,sm,fillValue)
+ do j=jsta,jend
+ do i=ista, iend
+ frzrn_bucket(i,j) = arrayr42d(i,j)
+ if (abs(arrayr42d(i,j)-fillValue) < small) frzrn_bucket(i,j) = spval
+ enddo
+ enddo
+ endif
+
! max hourly surface precipitation rate
if(trim(fieldname)=='pratemax') then
!$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,spval,prate_max,arrayr42d,sm,fillValue)
diff --git a/module_fcst_grid_comp.F90 b/module_fcst_grid_comp.F90
index 0fc3fff42..c8df6acc0 100644
--- a/module_fcst_grid_comp.F90
+++ b/module_fcst_grid_comp.F90
@@ -65,7 +65,7 @@ module module_fcst_grid_comp
nbdlphys, iau_offset
use module_fv3_config, only: dt_atmos, fcst_mpi_comm, fcst_ntasks, &
quilting, calendar, cpl_grid_id, &
- cplprint_flag, restart_endfcst
+ cplprint_flag
use get_stochy_pattern_mod, only: write_stoch_restart_atm
use module_cplfields, only: nExportFields, exportFields, exportFieldsInfo, &
@@ -73,7 +73,6 @@ module module_fcst_grid_comp
use module_cplfields, only: realizeConnectedCplFields
use atmos_model_mod, only: setup_exportdata
- use CCPP_data, only: GFS_control
!
!-----------------------------------------------------------------------
!
@@ -90,7 +89,7 @@ module module_fcst_grid_comp
type(ESMF_GridComp),dimension(:),allocatable :: fcstGridComp
integer :: ngrids, mygrid
- integer :: intrm_rst, n_atmsteps
+ integer :: n_atmsteps
!----- coupled model data -----
@@ -744,20 +743,7 @@ subroutine fcst_initialize(fcst_comp, importState, exportState, clock, rc)
endif
endif
! if to write out restart at the end of forecast
- restart_endfcst = .false.
- if ( ANY(frestart(:) == total_inttime) ) restart_endfcst = .true.
-! frestart only contains intermediate restart
- do i=1,size(frestart)
- if(frestart(i) == total_inttime) then
- frestart(i) = 0
- exit
- endif
- enddo
- if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'restart_endfcst=',restart_endfcst, &
- 'total_inttime=',total_inttime
-! if there is restart writing during integration
- intrm_rst = 0
- if (frestart(1)>0) intrm_rst = 1
+ if (mype == 0) print *,'frestart=',frestart(1:10)/3600, 'total_inttime=',total_inttime
!------ initialize component models ------
@@ -1255,9 +1241,8 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
!--- intermediate restart
- if (intrm_rst>0) then
- call get_time(Atmos%Time - Atmos%Time_init, seconds)
- if (ANY(frestart(:) == seconds)) then
+ call get_time(Atmos%Time - Atmos%Time_init, seconds)
+ if (ANY(frestart(:) == seconds)) then
if (mype == 0) write(*,*)'write out restart at n_atmsteps=',n_atmsteps,' seconds=',seconds, &
'integration length=',n_atmsteps*dt_atmos/3600.
@@ -1279,7 +1264,6 @@ subroutine fcst_run_phase_2(fcst_comp, importState, exportState,clock,rc)
'Current model time: year, month, day, hour, minute, second'
close( unit )
endif
- endif
endif
if (mype == 0) write(*,'(A,I16,A,F16.6)')'PASS: fcstRUN phase 2, n_atmsteps = ', &
@@ -1323,24 +1307,6 @@ subroutine fcst_finalize(fcst_comp, importState, exportState,clock,rc)
call atmos_model_end (Atmos)
-!*** write restart file
- if( restart_endfcst ) then
- call get_date (Atmos%Time, date(1), date(2), date(3), &
- date(4), date(5), date(6))
- call mpp_set_current_pelist()
- if (mpp_pe() == mpp_root_pe())then
- open( newunit=unit, file='RESTART/coupler.res' )
- write( unit, '(i6,8x,a)' )calendar_type, &
- '(Calendar: no_calendar=0, thirty_day_months=1, julian=2, gregorian=3, noleap=4)'
-
- write( unit, '(6i6,8x,a)' )date_init, &
- 'Model start time: year, month, day, hour, minute, second'
- write( unit, '(6i6,8x,a)' )date, &
- 'Current model time: year, month, day, hour, minute, second'
- close( unit )
- endif
- endif
-
call diag_manager_end (Atmos%Time)
call fms_end
diff --git a/module_fv3_config.F90 b/module_fv3_config.F90
index bb3546772..a62800b34 100644
--- a/module_fv3_config.F90
+++ b/module_fv3_config.F90
@@ -20,7 +20,6 @@ module module_fv3_config
integer :: cpl_grid_id
logical :: cplprint_flag
logical :: quilting, output_1st_tstep_rst
- logical :: restart_endfcst
!
real,dimension(:),allocatable :: output_fh
character(esmf_maxstr),dimension(:),allocatable :: filename_base
diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90
index 4e2ab8b59..009f04b87 100644
--- a/moving_nest/fv_moving_nest.F90
+++ b/moving_nest/fv_moving_nest.F90
@@ -54,7 +54,7 @@ module fv_moving_nest_mod
use block_control_mod, only : block_control_type
use fms_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE, clock_flag_default, CLOCK_SUBCOMPONENT
- use mpp_mod, only : mpp_pe, mpp_sync, mpp_sync_self, mpp_send, mpp_error, NOTE, FATAL
+ use mpp_mod, only : mpp_pe, mpp_sync, mpp_sync_self, mpp_send, mpp_error, NOTE, FATAL, WARNING
use mpp_domains_mod, only : mpp_update_domains, mpp_get_data_domain, mpp_get_global_domain
use mpp_domains_mod, only : mpp_define_nest_domains, mpp_shift_nest_domains, nest_domain_type, domain2d
use mpp_domains_mod, only : mpp_get_C2F_index, mpp_update_nest_fine
@@ -261,6 +261,12 @@ subroutine mn_prog_fill_nest_halos_from_parent(Atm, n, child_grid_num, is_fine_p
call fill_nest_halos_from_parent("q", Atm(n)%q, interp_type, Atm(child_grid_num)%neststruct%wt_h, &
Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position, nz)
+ ! Interpolate terrain from coarse grid
+ if (Moving_nest(n)%mn_flag%terrain_smoother .eq. 4) then
+ call fill_nest_halos_from_parent("phis", Atm(n)%phis, interp_type, Atm(child_grid_num)%neststruct%wt_h, &
+ Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position)
+ endif
+
! Move the A-grid winds. TODO consider recomputing them from D grid instead
call fill_nest_halos_from_parent("ua", Atm(n)%ua, interp_type, Atm(child_grid_num)%neststruct%wt_h, &
Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position, nz)
@@ -350,6 +356,8 @@ subroutine mn_prog_fill_intern_nest_halos(Atm, domain_fine, is_fine_pe)
integer :: this_pe
type(fv_moving_nest_prog_type), pointer :: mn_prog
+ integer :: child_grid_num = 2
+
mn_prog => Moving_nest(2)%mn_prog ! TODO allow nest number to vary
this_pe = mpp_pe()
@@ -360,6 +368,10 @@ subroutine mn_prog_fill_intern_nest_halos(Atm, domain_fine, is_fine_pe)
call mn_var_fill_intern_nest_halos(Atm%delp, domain_fine, is_fine_pe)
call mn_var_fill_intern_nest_halos(mn_prog%delz, domain_fine, is_fine_pe)
+ if (Moving_nest(child_grid_num)%mn_flag%terrain_smoother .eq. 4) then
+ call mn_var_fill_intern_nest_halos(Atm%phis, domain_fine, is_fine_pe)
+ endif
+
call mn_var_fill_intern_nest_halos(Atm%ua, domain_fine, is_fine_pe)
call mn_var_fill_intern_nest_halos(Atm%va, domain_fine, is_fine_pe)
@@ -596,14 +608,27 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de
if (first_nest_move) then
if (use_timers) call mpp_clock_begin (id_load1)
+ parent_geo%nxp = Atm(1)%npx
+ parent_geo%nyp = Atm(1)%npy
+
+ parent_geo%nx = parent_geo%nxp - 1
+ parent_geo%ny = parent_geo%nyp - 1
+
call mn_static_filename(surface_dir, parent_tile, 'grid', 1, grid_filename)
- call load_nest_latlons_from_nc(grid_filename, Atm(1)%npx, Atm(1)%npy, 1, pelist, &
+ call load_nest_latlons_from_nc(grid_filename, parent_geo%nxp, parent_geo%nyp, 1, pelist, &
parent_geo, p_istart_fine, p_iend_fine, p_jstart_fine, p_jend_fine)
+
! These are saved between timesteps in fv_moving_nest_main.F90
- allocate(p_grid(1:parent_geo%nxp, 1:parent_geo%nyp,2))
- allocate(p_grid_u(1:parent_geo%nxp, 1:parent_geo%nyp+1,2))
- allocate(p_grid_v(1:parent_geo%nxp+1, 1:parent_geo%nyp,2))
+ allocate(p_grid(parent_geo%nx, parent_geo%ny,2))
+ allocate(p_grid_u(parent_geo%nx, parent_geo%ny+1,2))
+ allocate(p_grid_v(parent_geo%nx+1, parent_geo%ny,2))
+
+ p_grid = 0.0
+ p_grid_u = 0.0
+ p_grid_v = 0.0
+
+
! These are big (parent grid size), and do not change during the model integration.
call assign_p_grids(parent_geo, p_grid, position)
@@ -616,11 +641,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de
if (use_timers) call mpp_clock_begin (id_load2)
- parent_geo%nxp = Atm(1)%npx
- parent_geo%nyp = Atm(1)%npy
-
- parent_geo%nx = Atm(1)%npx - 1
- parent_geo%ny = Atm(1)%npy - 1
!===========================================================
! Begin tile_geo per PE.
@@ -772,6 +792,8 @@ subroutine mn_latlon_read_hires_parent(npx, npy, refine, pelist, fp_super_tile_g
integer :: fp_super_istart_fine, fp_super_jstart_fine,fp_super_iend_fine, fp_super_jend_fine
character(len=256) :: grid_filename
+ integer :: i, j
+
call mn_static_filename(surface_dir, parent_tile, 'grid', refine, grid_filename)
call load_nest_latlons_from_nc(trim(grid_filename), npx, npy, refine, pelist, &
@@ -910,7 +932,7 @@ end subroutine mn_static_read_hires_r8
!>@brief The subroutine 'mn_meta_recalc' recalculates nest halo weights
subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, parent_geo, fp_super_tile_geo, &
- is_fine_pe, nest_domain, position, p_grid, n_grid, wt, istart_coarse, jstart_coarse)
+ is_fine_pe, nest_domain, position, p_grid, n_grid, wt, istart_coarse, jstart_coarse, ind)
integer, intent(in) :: delta_i_c, delta_j_c !< Nest motion in delta i,j
integer, intent(in) :: x_refine, y_refine !< Nest refinement
type(grid_geometry), intent(inout) :: tile_geo, parent_geo, fp_super_tile_geo !< tile geometries
@@ -921,8 +943,10 @@ subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, p
real, allocatable, intent(inout) :: wt(:,:,:) !< Interpolation weights
integer, intent(inout) :: position !< Stagger
integer, intent(in) :: istart_coarse, jstart_coarse !< Initian nest offsets
+ integer, allocatable, intent(in) :: ind(:,:,:)
type(bbox) :: wt_fine, wt_coarse
+ integer :: istag, jstag
integer :: this_pe
this_pe = mpp_pe()
@@ -936,17 +960,41 @@ subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, p
!!
!!===========================================================
+ if (position .eq. CENTER) then
+ istag = 0
+ jstag = 0
+ elseif (position .eq. NORTH) then
+ ! for p_grid_u
+ istag = 1
+ jstag = 0
+
+ !! This aligns with boundary.F90 settings
+ !!istag = 0
+ !!jstag = 1
+
+ elseif (position .eq. EAST) then
+ ! for p_grid_v
+ istag = 0
+ jstag = 1
+
+ !! This aligns with boundary.F90 settings
+ !istag = 1
+ !jstag = 0
+
+ endif
+
+
call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, EAST, position)
- call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine)
+ call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind)
call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, WEST, position)
- call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine)
+ call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind)
call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, NORTH, position)
- call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine)
+ call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind)
call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, SOUTH, position)
- call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine)
+ call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind)
endif
@@ -1034,6 +1082,11 @@ subroutine mn_prog_shift_data(Atm, n, child_grid_num, wt_h, wt_u, wt_v, &
call mn_var_shift_data(mn_prog%delz, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, &
delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz)
+ if (Moving_nest(n)%mn_flag%terrain_smoother .eq. 4) then
+ call mn_var_shift_data(Atm(n)%phis, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, &
+ delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position)
+ endif
+
call mn_var_shift_data(Atm(n)%ua, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, &
delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz)
@@ -2371,36 +2424,123 @@ subroutine assign_p_grids(parent_geo, p_grid, position)
integer :: i,j
+ integer :: num_zeros, num_vals, full_zeros
+
+ num_zeros = 0
+ full_zeros = 0
+ num_vals = 0
+
if (position == CENTER) then
- do j = 1, parent_geo%ny
- do i = 1, parent_geo%nx
+ do j = 1, ubound(p_grid,2)
+ do i = 1, ubound(p_grid,1)
! centered grid version
- p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j)
- p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j)
+ p_grid(i, j, :) = 0.0
+
+ if (2*i .gt. ubound(parent_geo%lons,1)) then
+ print '("[ERROR] WDR PG_CLONi npe=",I0," 2*i=",I0," ubound=",I0)', mpp_pe(), 2*i, ubound(parent_geo%lons,1)
+ elseif (2*j .gt. ubound(parent_geo%lons,2)) then
+ print '("[ERROR] WDR PG_CLONj npe=",I0," 2*j=",I0," ubound=",I0)', mpp_pe(), 2*j, ubound(parent_geo%lons,2)
+ else
+ p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j)
+ p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j)
+ endif
enddo
enddo
+
+
+ do i = 1, ubound(p_grid,1)
+ do j = 1, ubound(p_grid,2)
+
+ if (p_grid(i,j,1) .eq. 0.0) then
+ num_zeros = num_zeros + 1
+ if (p_grid(i,j,2) .eq. 0.0) then
+ full_zeros = full_zeros + 1
+ print '("[INFO] WDR set p_grid FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j
+ endif
+ endif
+ if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1
+ enddo
+ enddo
+
+ if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2)
+
! u(npx, npy+1)
elseif (position == NORTH) then ! u wind on D-stagger
- do j = 1, parent_geo%ny
- do i = 1, parent_geo%nx
- ! centered grid version
- p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j-1)
+ do j = 1, ubound(p_grid,2)
+ do i = 1, ubound(p_grid,1)
+ ! u grid version
+ p_grid(i, j, :) = 0.0
+ if (2*i .gt. ubound(parent_geo%lons,1)) then
+ print '("[ERROR] WDR PG_ULONi npe=",I0," 2*i=",I0," ubound=",I0)', mpp_pe(), 2*i, ubound(parent_geo%lons,1)
+ elseif (2*j-1 .gt. ubound(parent_geo%lons,2)) then
+ print '("[ERROR] WDR PG_ULONj npe=",I0," 2*j-1=",I0," ubound=",I0)', mpp_pe(), 2*j-1, ubound(parent_geo%lons,2)
+ else
+
+ ! This seems correct
+ p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j-1)
p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j-1)
- enddo
+ endif
enddo
+ enddo
+
+
+ do i = 1, ubound(p_grid,1)
+ do j = 1, ubound(p_grid,2)
+
+ if (p_grid(i,j,1) .eq. 0.0) then
+ num_zeros = num_zeros + 1
+ if (p_grid(i,j,2) .eq. 0.0) then
+ full_zeros = full_zeros + 1
+ print '("[INFO] WDR set p_grid_u FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j
+ endif
+ endif
+ if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1
+ enddo
+ enddo
+
+ if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid_u npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2)
+
! v(npx+1, npy)
elseif (position == EAST) then ! v wind on D-stagger
- do j = 1, parent_geo%ny
- do i = 1, parent_geo%nx
- ! centered grid version
- p_grid(i, j, 1) = parent_geo%lons(2*i-1, 2*j)
- p_grid(i, j, 2) = parent_geo%lats(2*i-1, 2*j)
+ do j = 1, ubound(p_grid,2)
+ do i = 1, ubound(p_grid,1)
+ ! v grid version
+ p_grid(i, j, :) = 0.0
+ if (2*i-1 .gt. ubound(parent_geo%lons,1)) then
+ print '("[ERROR] WDR PG_VLONi npe=",I0," 2*i-1=",I0," ubound=",I0)', mpp_pe(), 2*i-1, ubound(parent_geo%lons,1)
+ elseif (2*j .gt. ubound(parent_geo%lons,2)) then
+ print '("[ERROR] WDR PG_VLONj npe=",I0," 2*j=",I0," ubound=",I0)', mpp_pe(), 2*j, ubound(parent_geo%lons,2)
+ else
+ ! This seems correct
+ p_grid(i, j, 1) = parent_geo%lons(2*i-1, 2*j)
+ p_grid(i, j, 2) = parent_geo%lats(2*i-1, 2*j)
+ endif
+ enddo
+ enddo
+
+
+ do i = 1, ubound(p_grid,1)
+ do j = 1, ubound(p_grid,2)
+
+ if (p_grid(i,j,1) .eq. 0.0) then
+ num_zeros = num_zeros + 1
+ if (p_grid(i,j,2) .eq. 0.0) then
+ full_zeros = full_zeros + 1
+ print '("[INFO] WDR set p_grid_v FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j
+ endif
+ endif
+ if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1
enddo
enddo
+
+ if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid_v npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2)
+
+
+
endif
-
+
end subroutine assign_p_grids
@@ -2447,28 +2587,65 @@ subroutine assign_n_grids(tile_geo, n_grid, position)
end subroutine assign_n_grids
+ subroutine calc_inside(p_grid, ic, jc, n_grid1, n_grid2, istag, jstag, is_inside, verbose)
+ real(kind=R_GRID), allocatable, intent(in) :: p_grid(:,:,:)
+ real(kind=R_GRID), intent(in) :: n_grid1, n_grid2
+ integer, intent(in) :: ic, jc, istag, jstag
+ logical, intent(out) :: is_inside
+ logical, intent(in) :: verbose
+
+
+ real(kind=R_GRID) :: max1, max2, min1, min2, eps
+
+ max1 = max(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1))
+ max2 = max(p_grid(ic,jc,2), p_grid(ic,jc+1,2), p_grid(ic,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc,2))
+
+ min1 = min(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1))
+ min2 = min(p_grid(ic,jc,2), p_grid(ic,jc+1,2), p_grid(ic,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc,2))
+
+ is_inside = .False.
+
+ eps = 0.00001
+ !eps = 0.000001
+
+ if (n_grid1 .le. max1+eps .and. n_grid1 .ge. min1-eps) then
+ if (n_grid2 .le. max2+eps .and. n_grid2 .ge. min2-eps) then
+ is_inside = .True.
+ !if (verbose) print '("[INFO] WDR is_inside TRUE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F12.8," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag
+ endif
+ endif
+
+ if (verbose .and. .not. is_inside) then
+ print '("[INFO] WDR is_inside FALSE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F10.4," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag
+ endif
+ end subroutine calc_inside
+
!>@brief The subroutine 'calc_nest_halo_weights' calculates the interpolation weights
!>@details Computationally demanding; target for optimization after nest moves
- subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine)
+ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind)
implicit none
type(bbox), intent(in) :: bbox_coarse, bbox_fine !< Bounding boxes of parent and nest
real(kind=R_GRID), allocatable, intent(in) :: p_grid(:,:,:), n_grid(:,:,:) !< Latlon rids of parent and nest in radians
real, allocatable, intent(inout) :: wt(:,:,:) !< Interpolation weight array
integer, intent(in) :: istart_coarse, jstart_coarse, x_refine, y_refine !< Offsets and nest refinements
+ integer, intent(in) :: istag, jstag !< Staggers
+ integer, allocatable, intent(in) :: ind(:,:,:)
- integer :: i,j, ic, jc
+ integer :: i, j, k, ic, jc
real :: dist1, dist2, dist3, dist4, sum
logical :: verbose = .false.
!logical :: verbose = .true.
-
+ logical :: is_inside, adjusted
integer :: this_pe
real(kind=R_GRID) :: pi = 4 * atan(1.0d0)
real :: pi180
real :: rad2deg, deg2rad
+ real :: old_weight(4), diff_weight(4)
+ integer :: di, dj
pi180 = pi / 180.0
deg2rad = pi / 180.0
@@ -2486,25 +2663,75 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is
!
do j = bbox_fine%js, bbox_fine%je
- ! F90 integer division truncates
- jc = jstart_coarse + (j + y_refine/2 + 1) / y_refine
do i = bbox_fine%is, bbox_fine%ie
- ic = istart_coarse + (i + x_refine/2 + 1) / x_refine
-
+ jc = ind(i,j,2) ! reset this if the UPDATE code altered it
+ ic = ind(i,j,1)
+ if (ic+1 .gt. ubound(p_grid, 1)) print '("[ERROR] WDR CALCWT off end of p_grid i npe=",I0," ic+1=",I0," bound=",I0)', mpp_pe(), ic+1, ubound(p_grid,1)
+ if (jc+1 .gt. ubound(p_grid, 2)) print '("[ERROR] WDR CALCWT off end of p_grid j npe=",I0," jc+1=",I0," bound=",I0)', mpp_pe(), jc+1, ubound(p_grid,2)
+
! dist2side_latlon takes points in longitude-latitude coordinates.
dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:))
dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:))
dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:))
dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:))
+ call calc_inside(p_grid, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .True.)
+
+! if (.not. is_inside) then
+! adjusted = .False.
+!
+! do di = -2,2
+! do dj = -2,1
+! if (.not. adjusted) then
+! call calc_inside(p_grid, ic+di, jc+dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .False.)
+! if (is_inside) then
+! ic = ic + di
+! jc = jc + dj
+!
+! print '("[INFO] WDR is_inside UPDATED npe=",I0," ic=",I0," jc=",I0," istart_coarse=",I0," jstart_coarse=",I0," i=",I0," j=",I0," di=",I0," dj=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, istart_coarse, jstart_coarse, i, j, di, dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag
+
+! dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:))
+! dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:))
+! dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:))
+! dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:))
+!
+! adjusted = .True.
+! endif
+! endif
+! enddo
+! enddo
+! if (.not. adjusted) print '("[ERROR] WDR is_inside UPDATE FAILED npe=",I0," i=",I0," j=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), i, j, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag
+!
+! endif
+
+ old_weight = wt(i,j,:)
+
wt(i,j,1)=dist2*dist3 ! ic, jc weight
wt(i,j,2)=dist3*dist4 ! ic, jc+1 weight
wt(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight
wt(i,j,4)=dist1*dist2 ! ic+1, jc weight
+ ! If sum is zero, this is a problem
+
sum=wt(i,j,1)+wt(i,j,2)+wt(i,j,3)+wt(i,j,4)
- wt(i,j,:)=wt(i,j,:)/sum
+ if (sum .eq. 0.0) then
+
+ call mpp_error(WARNING, "WARNING: calc_nest_halo_weights sum of weights is zero.")
+ wt(i,j,:) = 0.25
+
+ else
+ wt(i,j,:)=wt(i,j,:)/sum
+ endif
+
+ diff_weight = old_weight - wt(i,j,:)
+ do k=1,4
+ if (abs(diff_weight(k)) .ge. 0.01) then
+ print '("[WARN] WDR DIFFWT npe=",I0," old_wt=",F10.6," wt(",I0,",",I0,",",I0,")=",F10.6," diff=",F10.6," istag=",I0," jstag=",I0)', &
+ mpp_pe(), old_weight(k), i, j, k, wt(i,j,k), diff_weight(k), istag, jstag
+
+ endif
+ enddo
enddo
enddo
endif
diff --git a/moving_nest/fv_moving_nest_main.F90 b/moving_nest/fv_moving_nest_main.F90
index c848e313f..4f1ce7aeb 100644
--- a/moving_nest/fv_moving_nest_main.F90
+++ b/moving_nest/fv_moving_nest_main.F90
@@ -682,6 +682,12 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c,
allocate(wt_v(Atm(child_grid_num)%bd%isd:Atm(child_grid_num)%bd%ied+1, Atm(child_grid_num)%bd%jsd:Atm(child_grid_num)%bd%jed, 4))
wt_v = real_snan
+
+ ! Fill in the local weights with the ones from Atm just to be safe
+ call fill_weight_grid(wt_h, Atm(n)%neststruct%wt_h)
+ call fill_weight_grid(wt_u, Atm(n)%neststruct%wt_u)
+ call fill_weight_grid(wt_v, Atm(n)%neststruct%wt_v)
+
else
allocate(wt_h(1,1,4))
wt_h = 0.0
@@ -912,20 +918,28 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c,
call mn_reset_phys_latlon(Atm, n, tile_geo, fp_super_tile_geo, Atm_block, IPD_control, IPD_data)
if (use_timers) call mpp_clock_end (id_movnest5_2)
- if (use_timers) call mpp_clock_begin (id_movnest5_3)
+ endif
!!============================================================================
!! Step 5.2 -- Fill the wt* variables for each stagger
!!============================================================================
+ call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h)
+ call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u)
+ call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v)
+ call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b)
+
+ if (is_fine_pe) then
+ if (use_timers) call mpp_clock_begin (id_movnest5_3)
+
call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, parent_geo, fp_super_tile_geo, &
- is_fine_pe, global_nest_domain, position, p_grid, n_grid, wt_h, istart_coarse, jstart_coarse)
+ is_fine_pe, global_nest_domain, position, p_grid, n_grid, wt_h, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_h)
call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo_u, parent_geo, fp_super_tile_geo, &
- is_fine_pe, global_nest_domain, position_u, p_grid_u, n_grid_u, wt_u, istart_coarse, jstart_coarse)
+ is_fine_pe, global_nest_domain, position_u, p_grid_u, n_grid_u, wt_u, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_u)
call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo_v, parent_geo, fp_super_tile_geo, &
- is_fine_pe, global_nest_domain, position_v, p_grid_v, n_grid_v, wt_v, istart_coarse, jstart_coarse)
+ is_fine_pe, global_nest_domain, position_v, p_grid_v, n_grid_v, wt_v, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_v)
if (use_timers) call mpp_clock_end (id_movnest5_3)
endif
@@ -936,10 +950,10 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c,
!! Step 5.3 -- Adjust the indices by the values of delta_i_c, delta_j_c
!!============================================================================
- call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h)
- call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u)
- call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v)
- call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b)
+ !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h)
+ !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u)
+ !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v)
+ !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b)
if (debug_sync) call mpp_sync(full_pelist) ! Used to make debugging easier. Can be removed.
@@ -996,6 +1010,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c,
case (2)
! Static nest smoothing algorithm - interpolation of coarse terrain in halo zone and 5 point blending zone of coarse and fine data
call set_blended_terrain(Atm(n), mn_static%parent_orog_grid, mn_static%orog_grid, x_refine, Atm(n)%bd%ng, 10, a_step)
+ case (4) ! Use coarse terrain; no-op here.
+ ;
case (5)
! 5 pt smoother. blend zone of 5 to match static nest
call set_smooth_nest_terrain(Atm(n), mn_static%orog_grid, x_refine, 5, Atm(n)%bd%ng, 5)
diff --git a/moving_nest/fv_moving_nest_types.F90 b/moving_nest/fv_moving_nest_types.F90
index 0a2a088fc..45bd12504 100644
--- a/moving_nest/fv_moving_nest_types.F90
+++ b/moving_nest/fv_moving_nest_types.F90
@@ -46,7 +46,7 @@ module fv_moving_nest_types_mod
! Moving Nest Namelist Variables
logical :: is_moving_nest = .false.
character(len=120) :: surface_dir = "INPUT/moving_nest"
- integer :: terrain_smoother = 1
+ integer :: terrain_smoother = 4
integer :: vortex_tracker = 0
integer :: ntrack = 1
integer :: corral_x = 5
@@ -221,7 +221,7 @@ module fv_moving_nest_types_mod
! Moving Nest Namelist Variables
logical, dimension(MAX_NNEST) :: is_moving_nest = .False.
character(len=120) :: surface_dir = "INPUT/moving_nest"
- integer, dimension(MAX_NNEST) :: terrain_smoother = 1 ! 0 -- all high-resolution data, 1 - static nest smoothing algorithm with blending zone of 5 points, 2 - blending zone of 10 points, 5 - 5 point smoother, 9 - 9 point smoother
+ integer, dimension(MAX_NNEST) :: terrain_smoother = 4 ! 0 -- all high-resolution data, 1 - static nest smoothing algorithm with blending zone of 5 points, 2 - blending zone of 10 points, 5 - 5 point smoother, 9 - 9 point smoother
integer, dimension(MAX_NNEST) :: vortex_tracker = 0 ! 0 - not a moving nest, tracker not needed
! 1 - prescribed nest moving
! 2 - following child domain center
diff --git a/moving_nest/fv_moving_nest_utils.F90 b/moving_nest/fv_moving_nest_utils.F90
index 609a327ce..11485aa88 100644
--- a/moving_nest/fv_moving_nest_utils.F90
+++ b/moving_nest/fv_moving_nest_utils.F90
@@ -930,7 +930,7 @@ subroutine load_nest_latlons_from_nc(nc_filename, nxp, nyp, refine, pelist, &
character(*), intent(in) :: nc_filename
integer, intent(in) :: nxp, nyp, refine
integer, allocatable, intent(in) :: pelist(:)
- type(grid_geometry), intent(out) :: fp_tile_geo
+ type(grid_geometry), intent(inout) :: fp_tile_geo
integer, intent(out) :: fp_istart_fine, fp_iend_fine, fp_jstart_fine, fp_jend_fine
!========================================================================================
diff --git a/upp b/upp
index e22724738..b37f8ab7b 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit e22724738fd104327fee7c3c7ffc805ccabd619f
+Subproject commit b37f8ab7b0f298346d79a37e0c5d4a64037fd4d4