Skip to content

Commit

Permalink
FV3: this commits #refs 65359,fix some issues when running FV3 with N…
Browse files Browse the repository at this point in the history 
…oah MP, #refs 65281, Issue about failure of hydrostatic run in FV3 and #refs 65794 IAU updates and adding enhanced divergence damping for EnKF DA

Change-Id: Ieb67ca9c62b54af3b627ceebaedf970a6b4d793a
  • Loading branch information
@junwang-noaa
junwang-noaa committed Jul 15, 2019
1 parent 3dc7e9b commit 64fc1f1
Show file tree
Hide file tree
Showing 14 changed files with 336 additions and 239 deletions.
26 changes: 24 additions & 2 deletions atmos_cubed_sphere/model/dyn_core.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,7 @@ module dyn_core_mod
real, allocatable :: rf(:)
integer:: k_rf = 0
logical:: RFF_initialized = .false.
logical:: first_call = .true.
integer :: kmax=1

contains
Expand Down Expand Up @@ -739,13 +740,14 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
endif

if (flagstruct%regional) call exch_uv(domain, bd, npz, vc, uc)
if (first_call .and. is_master() .and. last_step) write(6,*) 'Sponge layer divergence damping coefficent:'

call timing_on('d_sw')
!$OMP parallel do default(none) shared(npz,flagstruct,nord_v,pfull,damp_vt,hydrostatic,last_step, &
!$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, &
!$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, &
!$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, &
!$OMP heat_source,diss_est) &
!$OMP heat_source,diss_est,ptop,first_call) &
!$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, &
!$OMP d_con_k,kgb, hord_m, hord_v, hord_t, hord_p, wk, heat_s,diss_e, z_rat)
do k=1,npz
Expand Down Expand Up @@ -782,14 +784,18 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
else
! Sponge layers with del-2 damping on divergence, vorticity, w, z, and air mass (delp).
! no special damping of potential temperature in sponge layers
! enhanced del-2 divergence damping has same vertical structure as Rayleigh
! damping if d2_bg_k2<=0.
if (flagstruct%d2_bg_k2 > 0) then ! old version, only applied at top two or three levels

if ( k==1 ) then
! Divergence damping:
nord_k=0; d2_divg = max(0.01, flagstruct%d2_bg, flagstruct%d2_bg_k1)
! Vertical velocity:
nord_w=0; damp_w = d2_divg
if ( flagstruct%do_vort_damp ) then
! damping on delp and vorticity:
nord_v(k)=0;
nord_v(k)=0;
#ifndef HIWPP
damp_vt(k) = 0.5*d2_divg
#endif
Expand All @@ -810,6 +816,21 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
nord_w=0; damp_w = d2_divg
d_con_k = 0.
endif
else ! new version, uses d2_bg_k1 and sponge layer vertical structure
if ( pfull(k) < flagstruct%rf_cutoff ) then
nord_k=0; nord_w=0
d2_divg = max(flagstruct%d2_bg, flagstruct%d2_bg_k1* &
sin(0.5*pi*log(flagstruct%rf_cutoff/pfull(k))/log(flagstruct%rf_cutoff/ptop))**2)
if (first_call .and. is_master() .and. last_step) write(6,*) k, 0.01*pfull(k), d2_divg
damp_w = d2_divg
if ( flagstruct%do_vort_damp ) then
! damping on delp and vorticity
nord_v(k)=0
damp_vt(k) = 0.5*d2_divg
endif
endif
endif

endif

if( hydrostatic .and. (.not.flagstruct%use_old_omega) .and. last_step ) then
Expand Down Expand Up @@ -1300,6 +1321,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,

!-----------------------------------------------------
enddo ! time split loop
first_call=.false.
!-----------------------------------------------------
if ( nq > 0 .and. .not. flagstruct%inline_q ) then
call timing_on('COMM_TOTAL')
Expand Down
4 changes: 3 additions & 1 deletion atmos_cubed_sphere/model/fv_dynamics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,7 +448,9 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
call moist_cp(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, q_con(is:ie,j,k), cvm)
#endif

do i=is,ie
dp1(i,j,k) = zvir*q(i,j,k,sphum)
enddo
enddo
enddo
else
Expand Down
2 changes: 1 addition & 1 deletion atmos_cubed_sphere/tools/fv_iau_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -485,7 +485,7 @@ subroutine interp_inc(field_name,var,jbeg,jend)
if (ierr == 0) then
call get_var3_r4( ncid, field_name, 1,im, jbeg,jend, 1,km, wk3 )
else
print *,'warning: no increment for ',trim(field_name),' found, assuming zero'
if (is_master()) print *,'warning: no increment for ',trim(field_name),' found, assuming zero'
wk3 = 0.
endif
do k=1,km
Expand Down
2 changes: 1 addition & 1 deletion atmos_cubed_sphere/tools/fv_treat_da_inc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -436,7 +436,7 @@ subroutine apply_inc_on_3d_scalar(field_name,var)
if (ierr == 0) then
call get_var3_r4( ncid, field_name, 1,im, jbeg,jend, 1,km, wk3 )
else
print *,'warning: no increment for ',trim(field_name),' found, assuming zero'
if (is_master()) print *,'warning: no increment for ',trim(field_name),' found, assuming zero'
wk3 = 0.
endif

Expand Down
2 changes: 2 additions & 0 deletions atmos_model.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,7 @@ module atmos_model_mod
logical :: regional ! true if domain is regional
logical :: nested ! true if there is a nest
integer :: mlon, mlat
integer :: iau_offset ! iau running window length
logical :: pe ! current pe.
real(kind=8), pointer, dimension(:) :: ak, bk
real, pointer, dimension(:,:) :: lon_bnd => null() ! local longitude axis grid box corners in radians.
Expand Down Expand Up @@ -553,6 +554,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
Init_parm%cdat(:) = cdat(:)
Init_parm%dt_dycore = dt_phys
Init_parm%dt_phys = dt_phys
Init_parm%iau_offset = Atmos%iau_offset
Init_parm%blksz => Atm_block%blksz
Init_parm%ak => Atmos%ak
Init_parm%bk => Atmos%bk
Expand Down
22 changes: 18 additions & 4 deletions fv3_cap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ module fv3gfs_cap_mod
write_fsyncflag, nsout_io, &
cen_lon, cen_lat, &
lon1, lat1, lon2, lat2, dlon, dlat, &
stdlat1, stdlat2, dx, dy
stdlat1, stdlat2, dx, dy, iau_offset
!
use module_fcst_grid_comp, only: fcstSS => SetServices, &
fcstGrid, numLevels, numSoilLayers, &
Expand Down Expand Up @@ -229,7 +229,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
type(ESMF_Time) :: CurrTime, starttime, StopTime
type(ESMF_Time) :: alarm_output_hf_ring, alarm_output_ring
type(ESMF_Time) :: alarm_output_hf_stop, alarm_output_stop
type(ESMF_TimeInterval) :: RunDuration, timeStep, rsthour
type(ESMF_TimeInterval) :: RunDuration, timeStep, rsthour, IAU_offsetTI
type(ESMF_Config) :: cf
type(ESMF_RegridMethod_Flag) :: regridmethod
type(ESMF_TimeInterval) :: earthStep
Expand Down Expand Up @@ -299,7 +299,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
default=.false., label ='output_1st_tstep_rst:',rc=rc)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return

if(mype == 0) print *,'af nems config,quilting=',quilting,'calendar=', trim(calendar)
call ESMF_ConfigGetAttribute(config=CF,value=iau_offset,default=0,label ='iau_offset:',rc=rc)
if (iau_offset < 0) iau_offset=0

if(mype == 0) print *,'af nems config,quilting=',quilting,'calendar=', trim(calendar),'iau_offset=',iau_offset
!
nfhout = 0 ; nfhmax_hf = 0 ; nfhout_hf = 0 ; nsout = 0
if ( quilting ) then
Expand Down Expand Up @@ -345,6 +348,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
end if
write_nemsioflip =.false.
write_fsyncflag =.false.

if(trim(output_grid) == 'gaussian_grid') then
call ESMF_ConfigGetAttribute(config=CF, value=imo, label ='imo:',rc=rc)
call ESMF_ConfigGetAttribute(config=CF, value=jmo, label ='jmo:',rc=rc)
Expand Down Expand Up @@ -572,6 +576,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
!
allocate(petList(wrttasks_per_group))
if(mype == 0) print *,'af allco wrtComp,write_groups=',write_groups

! set up ESMF time interval at center of iau window
call ESMF_TimeIntervalSet(IAU_offsetTI, h=iau_offset, rc=rc)
!
allocate(originPetList(num_pes_fcst+wrttasks_per_group))
allocate(targetPetList(num_pes_fcst+wrttasks_per_group))
Expand Down Expand Up @@ -735,6 +742,9 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
if (currtime <= starttime+output_hfmax) then
nhf = (currtime-starttime)/output_interval_hf
alarm_output_hf_ring = startTime + (nhf+1_ESMF_KIND_I4)*output_interval_hf
if(iau_offset > 0) then
alarm_output_hf_ring = startTime + IAU_offsetTI
endif
alarm_output_hf = ESMF_AlarmCreate(clock_fv3,name='ALARM_OUTPUT_HF', &
ringTime =alarm_output_hf_ring, &
ringInterval =output_interval_hf, & !<-- Time interval between
Expand All @@ -744,14 +754,17 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
rc =RC)
if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return

alarm_output_ring = starttime + output_hfmax + output_interval
alarm_output_ring = startTime + output_hfmax + output_interval
else
nrg = (currtime-starttime-output_hfmax)/output_interval
alarm_output_ring = startTime + output_hfmax + (nrg+1_ESMF_KIND_I4) * output_interval
endif
else
nrg = (currtime-starttime)/output_interval
alarm_output_ring = startTime + (nrg+1_ESMF_KIND_I4) * output_interval
if(iau_offset > 0) then
alarm_output_ring = startTime + IAU_offsetTI
endif
endif

call ESMF_TimeIntervalSet(output_interval, h=nfhout, m=nfmout, &
Expand Down Expand Up @@ -1037,6 +1050,7 @@ subroutine ModelAdvance(gcomp, rc)
output: IF(lalarm .or. na==first_kdt ) then

timerhi = mpi_wtime()
! if (mype == 0 .or. mype == lead_wrttask(1)) print *,' aft fcst run alarm is on, na=',na,'mype=',mype
do i=1, FBCount
!
! get fcst fieldbundle
Expand Down
19 changes: 16 additions & 3 deletions gfsphysics/GFS_layer/GFS_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,8 +186,8 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, &
Init_parm%levs, Init_parm%cnx, Init_parm%cny, &
Init_parm%gnx, Init_parm%gny, &
Init_parm%dt_dycore, Init_parm%dt_phys, &
Init_parm%bdat, Init_parm%cdat, &
Init_parm%tracer_names, &
Init_parm%iau_offset, Init_parm%bdat, &
Init_parm%cdat, Init_parm%tracer_names, &
Init_parm%input_nml_file, Init_parm%tile_num &
#ifdef CCPP
,Init_parm%ak, Init_parm%bk, Init_parm%blksz, &
Expand Down Expand Up @@ -505,7 +505,8 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, &


!--- local variables
integer :: nb, nblks, k, kdt_rad, blocksize
integer :: nb, nblks, k, kdt_rad, kdt_iau, blocksize
logical :: iauwindow_center
real(kind=kind_phys) :: rinc(5)
real(kind=kind_phys) :: sec, sec_zero
real(kind=kind_phys), parameter :: cn_hr = 3600._kind_phys
Expand Down Expand Up @@ -586,6 +587,18 @@ subroutine GFS_time_vary_step (Model, Statein, Stateout, Sfcprop, Coupling, &
enddo
endif
endif
!
if (Model%iau_offset > 0) then
kdt_iau = nint(Model%iau_offset*con_hr/Model%dtp)
if (Model%kdt == kdt_iau+1) then
iauwindow_center = .true.
do nb = 1,nblks
call Diag(nb)%rad_zero (Model)
call Diag(nb)%phys_zero (Model,iauwindow_center=iauwindow_center)
enddo
if(Model%me == Model%master) print *,'in gfs_driver, at iau_center, zero out rad/phys accumulated diag fields, kdt=',Model%kdt,'kdt_iau=',kdt_iau,'iau_offset=',Model%iau_offset
endif
endif
call run_stochastic_physics(nblks,Model,Grid(:),Coupling(:))

if(Model%do_ca)then
Expand Down
53 changes: 31 additions & 22 deletions gfsphysics/GFS_layer/GFS_typedefs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,7 @@ module GFS_typedefs
!! | logunit | | fortran unit number for writing logfile | none | 0 | integer | | none | F |
!! | bdat | | model begin date in GFS format (same as idat) | none | 0 | integer | | none | F |
!! | cdat | | model current date in GFS format (same as jdat) | none | 0 | integer | | none | F |
!! | iau_offset | | iau running window length in hour | none | 0 | integer | | none | F |
!! | dt_dycore | | dynamics time step in seconds | s | 0 | real | kind_phys | none | F |
!! | dt_phys | | physics time step in seconds | s | 0 | real | kind_phys | none | F |
!! | restart | | flag for restart (warmstart) or coldstart | flag | 0 | logical | | none | F |
Expand Down Expand Up @@ -196,6 +197,7 @@ module GFS_typedefs
integer :: logunit !< fortran unit number for writing logfile
integer :: bdat(8) !< model begin date in GFS format (same as idat)
integer :: cdat(8) !< model current date in GFS format (same as jdat)
integer :: iau_offset !< iau running window length
real(kind=kind_phys) :: dt_dycore !< dynamics time step in seconds
real(kind=kind_phys) :: dt_phys !< physics time step in seconds
#ifdef CCPP
Expand Down Expand Up @@ -1654,6 +1656,7 @@ module GFS_typedefs
#endif

!--- IAU
integer :: iau_offset
real(kind=kind_phys) :: iau_delthrs ! iau time interval (to scale increments) in hours
character(len=240) :: iau_inc_files(7)! list of increment files
real(kind=kind_phys) :: iaufhrs(7) ! forecast hours associated with increment files
Expand Down Expand Up @@ -3676,7 +3679,8 @@ end subroutine coupling_create
subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
logunit, isc, jsc, nx, ny, levs, &
cnx, cny, gnx, gny, dt_dycore, &
dt_phys, idat, jdat, tracer_names, &
dt_phys, iau_offset, idat, jdat, &
tracer_names, &
input_nml_file, tile_num &
#ifdef CCPP
,ak, bk, blksz, &
Expand Down Expand Up @@ -3721,6 +3725,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
integer, intent(in) :: gny
real(kind=kind_phys), intent(in) :: dt_dycore
real(kind=kind_phys), intent(in) :: dt_phys
integer, intent(in) :: iau_offset
integer, intent(in) :: idat(8)
integer, intent(in) :: jdat(8)
character(len=32), intent(in) :: tracer_names(:)
Expand Down Expand Up @@ -4069,7 +4074,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &


!--- IAU options
real(kind=kind_phys) :: iau_delthrs = 6 ! iau time interval (to scale increments)
real(kind=kind_phys) :: iau_delthrs = 0 ! iau time interval (to scale increments)
character(len=240) :: iau_inc_files(7)='' ! list of increment files
real(kind=kind_phys) :: iaufhrs(7)=-1 ! forecast hours associated with increment files
logical :: iau_filter_increments = .false. ! filter IAU increments
Expand Down Expand Up @@ -4292,6 +4297,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%idate(2) = Model%idat(2)
Model%idate(3) = Model%idat(3)
Model%idate(4) = Model%idat(1)
Model%iau_offset = iau_offset

!--- radiation control parameters
Model%fhswr = fhswr
Expand Down Expand Up @@ -4620,6 +4626,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
Model%iau_inc_files = iau_inc_files
Model%iau_delthrs = iau_delthrs
Model%iau_filter_increments = iau_filter_increments
if(Model%me==0) print *,' model init,iaufhrs=',Model%iaufhrs

!--- tracer handling
Model%ntrac = size(tracer_names)
Expand Down Expand Up @@ -6087,9 +6094,10 @@ subroutine diag_create (Diag, IM, Model)
if (Model%cplchm) call Diag%chem_init(IM,Model)

call Diag%rad_zero (Model)
! print *,'in diag_create, call phys_zero'
! if(Model%me==0) print *,'in diag_create, call rad_zero'
linit = .true.
call Diag%phys_zero (Model, linit=linit)
! if(Model%me==0) print *,'in diag_create, call phys_zero'
linit = .false.

end subroutine diag_create
Expand All @@ -6116,10 +6124,12 @@ end subroutine diag_rad_zero
!------------------------
! GFS_diag%phys_zero
!------------------------
subroutine diag_phys_zero (Diag, Model, linit)
subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
class(GFS_diag_type) :: Diag
type(GFS_control_type), intent(in) :: Model
logical,optional, intent(in) :: linit
logical,optional, intent(in) :: linit, iauwindow_center

logical set_totprcp

!--- In/Out
Diag%srunoff = zero
Expand Down Expand Up @@ -6295,23 +6305,22 @@ subroutine diag_phys_zero (Diag, Model, linit)
!-----------------------------

! max hourly diagnostics
Diag%refl_10cm = zero
Diag%refdmax = -35.
Diag%refdmax263k = -35.
Diag%t02max = -999.
Diag%t02min = 999.
Diag%rh02max = -999.
Diag%rh02min = 999.
if (present(linit)) then
if (linit) then
Diag%totprcp = zero
Diag%cnvprcp = zero
Diag%totice = zero
Diag%totsnw = zero
Diag%totgrp = zero
! if(Model%me == Model%master) print *,'in diag_phys_zero, called in init step,set precip diag variable to zero',&
! 'size(Diag%totprcp)=',size(Diag%totprcp),'me=',Model%me,'kdt=',Model%kdt
endif
Diag%refl_10cm = zero
Diag%refdmax = -35.
Diag%refdmax263k = -35.
Diag%t02max = -999.
Diag%t02min = 999.
Diag%rh02max = -999.
Diag%rh02min = 999.
set_totprcp = .false.
if (present(linit) ) set_totprcp = linit
if (present(iauwindow_center) ) set_totprcp = iauwindow_center
if (set_totprcp) then
Diag%totprcp = zero
Diag%cnvprcp = zero
Diag%totice = zero
Diag%totsnw = zero
Diag%totgrp = zero
endif
end subroutine diag_phys_zero

Expand Down
Loading

0 comments on commit 64fc1f1

Please sign in to comment.