Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update dtc/develop from dev/emc 2020/03/04 #9

Merged
merged 12 commits into from
Mar 9, 2020
59 changes: 48 additions & 11 deletions driver/fvGFS/atmosphere.F90
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ module atmosphere_mod
use fv_nggps_diags_mod, only: fv_nggps_diag_init, fv_nggps_diag, fv_nggps_tavg
use fv_restart_mod, only: fv_restart, fv_write_restart
use fv_timing_mod, only: timing_on, timing_off
use fv_mp_mod, only: switch_current_Atm
use fv_mp_mod, only: switch_current_Atm, is_master
use fv_sg_mod, only: fv_subgrid_z
use fv_update_phys_mod, only: fv_update_phys
use fv_nwp_nudge_mod, only: fv_nwp_nudge_init, fv_nwp_nudge_end, do_adiabatic_init
Expand All @@ -194,6 +194,7 @@ module atmosphere_mod
a_step, p_step, current_time_in_seconds

use mpp_domains_mod, only: mpp_get_data_domain, mpp_get_compute_domain
use fv_grid_utils_mod, only: g_sum

implicit none
private
Expand Down Expand Up @@ -1395,6 +1396,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
integer :: i, j, ix, k, k1, n, w_diff, nt_dyn, iq
integer :: nb, blen, nwat, dnats, nq_adv
real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt
real psum, qsum, psumb, qsumb, betad
Time_prev = Time
Time_next = Time + Time_step_atmos
rdt = 1.d0 / dt_atmos
Expand All @@ -1407,6 +1409,18 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers')

if (IAU_Data%in_interval) then
if (IAU_Data%drymassfixer) then
! global mean total pressure and water before IAU
psumb = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.)
qsumb = g_sum(Atm(n)%domain,&
sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.)
if (is_master()) then
print *,'dry ps before IAU/physics',psumb+Atm(n)%ptop-qsumb
endif
endif

! IAU increments are in units of 1/sec

! add analysis increment to u,v,t tendencies
Expand Down Expand Up @@ -1469,11 +1483,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys)

do k = 1, npz
if(flip_vc) then
k1 = npz+1-k !reverse the k direction
else
k1 = k
endif
if(flip_vc) then
k1 = npz+1-k !reverse the k direction
else
k1 = k
endif
do ix = 1, blen
i = Atm_block%index(nb)%ii(ix)
j = Atm_block%index(nb)%jj(ix)
Expand Down Expand Up @@ -1515,11 +1529,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc
!--- See Note in statein...
do iq = nq+1, ncnst
do k = 1, npz
if(flip_vc) then
k1 = npz+1-k !reverse the k direction
else
k1 = k
endif
if(flip_vc) then
k1 = npz+1-k !reverse the k direction
else
k1 = k
endif
do ix = 1, blen
i = Atm_block%index(nb)%ii(ix)
j = Atm_block%index(nb)%jj(ix)
Expand All @@ -1530,6 +1544,29 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc

enddo ! nb-loop

! dry mass fixer in IAU interval following
! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x
if (IAU_Data%in_interval .and. IAU_data%drymassfixer) then
! global mean total pressure
psum = g_sum(Atm(n)%domain,sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz),dim=3),&
isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.)
! global mean total water (before adjustment)
qsum = g_sum(Atm(n)%domain,&
sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1,reproduce=.true.)
betad = (psum - (psumb - qsumb))/qsum
if (is_master()) then
print *,'dry ps after IAU/physics',psum+Atm(n)%ptop-qsum
endif
Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat)
!qsum = g_sum(Atm(n)%domain,&
! sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(Atm(n)%q(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),&
! isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1)
!if (is_master()) then
! print *,'dry ps after iau_drymassfixer',psum+Atm(n)%ptop-qsum
!endif
endif

call timing_off('GFS_TENDENCIES')

w_diff = get_tracer_index (MODEL_ATMOS, 'w_diff' )
Expand Down
2 changes: 2 additions & 0 deletions tools/fv_iau_mod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,7 @@ module fv_iau_mod
real,allocatable :: delz_inc(:,:,:)
real,allocatable :: tracer_inc(:,:,:,:)
logical :: in_interval = .false.
logical :: drymassfixer = .false.
end type iau_external_data_type
type iau_state_type
type(iau_internal_data_type):: inc1
Expand Down Expand Up @@ -280,6 +281,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm)
call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2)))
endif
! print*,'in IAU init',dt,rdt
IAU_data%drymassfixer = IPD_control%iau_drymassfixer

end subroutine IAU_initialize

Expand Down