From cc28aa5aa94817b24aea033cc19c1792e5927459 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Tue, 5 Nov 2019 15:37:56 +0000 Subject: [PATCH 01/13] restart reproducibility fix for regional fv3 --- model/fv_regional_bc.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 8b6d11db8..67d207c59 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -1183,6 +1183,12 @@ subroutine start_regional_restart(Atm & endif ! !----------------------------------------------------------------------- +!*** Get the source of the input data. +!----------------------------------------------------------------------- +! + call get_data_source(data_source,Atm%flagstruct%regional) +! +!----------------------------------------------------------------------- !*** Preliminary setup for the forecast. !----------------------------------------------------------------------- ! From 0e84f88b494b9e0a4097da50abe6b143330e8a2f Mon Sep 17 00:00:00 2001 From: junwang-noaa <37633869+junwang-noaa@users.noreply.github.com> Date: Tue, 24 Dec 2019 12:40:12 -0500 Subject: [PATCH 02/13] make writing last restart file optional (#9) --- driver/fvGFS/atmosphere.F90 | 5 +++-- model/fv_control.F90 | 5 +++-- tools/fv_restart.F90 | 14 +++++++++----- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 48b8d1d54..0d81b6cb1 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -712,7 +712,7 @@ end subroutine atmosphere_dynamics !>@brief The subroutine 'atmosphere_end' is an API for the termination of the !! FV3 dynamical core responsible for writing out a restart and final diagnostic state. - subroutine atmosphere_end (Time, Grid_box) + subroutine atmosphere_end (Time, Grid_box, restart_endfcst) #ifdef CCPP #ifdef STATIC ! For static builds, the ccpp_physics_{init,run,finalize} calls @@ -727,6 +727,7 @@ subroutine atmosphere_end (Time, Grid_box) #endif type (time_type), intent(in) :: Time type(grid_box_type), intent(inout) :: Grid_box + logical, intent(in) :: restart_endfcst #ifdef CCPP integer :: ierr @@ -754,7 +755,7 @@ subroutine atmosphere_end (Time, Grid_box) call timing_off('FV_DIAG') endif - call fv_end(Atm, grids_on_this_pe) + call fv_end(Atm, grids_on_this_pe, restart_endfcst) deallocate (Atm) deallocate( u_dt, v_dt, t_dt, pref, dum1d ) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index e87b01963..6a36bd835 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -592,17 +592,18 @@ end subroutine fv_init !>@brief The subroutine 'fv_end' terminates FV3, deallocates memory, !! saves restart files, and stops I/O. - subroutine fv_end(Atm, grids_on_this_pe) + subroutine fv_end(Atm, grids_on_this_pe, restart_endfcst) type(fv_atmos_type), intent(inout) :: Atm(:) logical, intent(INOUT) :: grids_on_this_pe(:) + logical, intent(in) :: restart_endfcst integer :: n call timing_off('TOTAL') call timing_prt( gid ) - call fv_restart_end(Atm, grids_on_this_pe) + call fv_restart_end(Atm, grids_on_this_pe, restart_endfcst) call fv_io_exit() ! Free temporary memory from sw_core routines diff --git a/tools/fv_restart.F90 b/tools/fv_restart.F90 index 1632212ed..72071cea6 100644 --- a/tools/fv_restart.F90 +++ b/tools/fv_restart.F90 @@ -1434,9 +1434,10 @@ end subroutine fv_write_restart !>@brief The subroutine 'fv_restart_end' writes ending restart files, !! terminates I/O, and prints out diagnostics including global totals !! and checksums. - subroutine fv_restart_end(Atm, grids_on_this_pe) + subroutine fv_restart_end(Atm, grids_on_this_pe, restart_endfcst) type(fv_atmos_type), intent(inout) :: Atm(:) logical, intent(INOUT) :: grids_on_this_pe(:) + logical, intent(in) :: restart_endfcst integer :: isc, iec, jsc, jec integer :: iq, n, ntileMe, ncnst, ntprog, ntdiag @@ -1516,10 +1517,13 @@ subroutine fv_restart_end(Atm, grids_on_this_pe) enddo - call fv_io_write_restart(Atm, grids_on_this_pe) - do n=1,ntileMe - if (Atm(n)%neststruct%nested .and. grids_on_this_pe(n)) call fv_io_write_BCs(Atm(n)) - end do + if ( restart_endfcst ) then + call fv_io_write_restart(Atm, grids_on_this_pe) +! print *,'af call fv_io_write_restart, restart_endfcst=',restart_endfcst + do n=1,ntileMe + if (Atm(n)%neststruct%nested .and. grids_on_this_pe(n)) call fv_io_write_BCs(Atm(n)) + end do + endif module_is_initialized = .FALSE. From a56907a44461c7151e0ba266e160c8f1a1685882 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 6 Jan 2020 08:48:15 -0700 Subject: [PATCH 03/13] Update dev/emc from dtc/develop: add dynamics code for Ferrier-Aligo microphysics (#11) --- model/fv_mapz.F90 | 27 ++++++++++++++++++ model/fv_sg.F90 | 72 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 99 insertions(+) diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index 4994538b9..2c72074a1 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -3465,12 +3465,25 @@ subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai enddo case(4) ! K_warm_rain with fake ice do i=is,ie +#ifndef CCPP qv(i) = q(i,j,k,sphum) qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) #ifdef MULTI_GASES cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + qd(i)*c_liq #else cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + qd(i)*c_liq +#endif +#else + qv(i) = q(i,j,k,sphum) + ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) + qs(i) = q(i,j,k,ice_wat) + qd(i) = ql(i) + qs(i) +#ifdef MULTI_GASES + cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice +#else + cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice +#endif + #endif enddo case(5) @@ -3574,12 +3587,26 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai enddo case(4) ! K_warm_rain scheme with fake ice do i=is,ie +#ifndef CCPP qv(i) = q(i,j,k,sphum) qd(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) #ifdef MULTI_GASES cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + qd(i)*c_liq #else cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + qd(i)*c_liq +#endif +#else + qv(i) = q(i,j,k,sphum) + ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat) + qs(i) = q(i,j,k,ice_wat) + qd(i) = ql(i) + qs(i) +#ifdef MULTI_GASES + cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice +#else + cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice +#endif + + #endif enddo case(5) diff --git a/model/fv_sg.F90 b/model/fv_sg.F90 index 594e54873..a88cfdfe1 100644 --- a/model/fv_sg.F90 +++ b/model/fv_sg.F90 @@ -299,6 +299,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat==4 ) then do i=is,ie +#ifndef CCPP q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) #ifdef MULTI_GASES cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq @@ -306,6 +307,20 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq +#endif + +#else + q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) + q_sol = q0(i,k,ice_wat) +#ifdef MULTI_GASES + cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#else + cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#endif + + #endif enddo elseif ( nwat==5 ) then @@ -382,7 +397,11 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==4 ) then do k=1,kbot do i=is,ie +#ifndef CCPP qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) +#else + qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) + q0(i,k,ice_wat) +#endif enddo enddo elseif ( nwat==5 ) then @@ -451,7 +470,12 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==3 ) then ! AM3/AM4 qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice +#ifndef CCPP qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) +#else + qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & + q0(i,km1,rainwat) +#endif elseif ( nwat==5 ) then ! K_warm_rain scheme with fake ice qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & q0(i,km1,snowwat) + q0(i,km1,rainwat) @@ -572,6 +596,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat == 4 ) then do i=is,ie +#ifndef CCPP q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) #ifdef MULTI_GASES cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq @@ -579,6 +604,19 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq +#endif +#else + q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) + q_sol = q0(i,kk,ice_wat) +#ifdef MULTI_GASES + cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#else + cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#endif + + #endif enddo elseif ( nwat == 5 ) then @@ -850,6 +888,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat==4 ) then do i=is,ie +#ifndef CCPP q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) #ifdef MULTI_GASES cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq @@ -857,6 +896,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,k,sphum)+q_liq))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq cvm(i) = (1.-(q0(i,k,sphum)+q_liq))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq +#endif +#else + q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat) + q_sol = q0(i,k,ice_wat) +#ifdef MULTI_GASES + cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#else + cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#endif + #endif enddo elseif ( nwat==5 ) then @@ -933,7 +984,11 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==4 ) then do k=1,kbot do i=is,ie +#ifndef CCPP qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) +#else + qcon(i,k) = q0(i,k,liq_wat) + q0(i,k,rainwat) + q0(i,k,ice_wat) +#endif enddo enddo elseif ( nwat==5 ) then @@ -998,7 +1053,11 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & elseif ( nwat==3 ) then ! AM3/AM4 qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) elseif ( nwat==4 ) then ! K_warm_rain scheme with fake ice +#ifndef CCPP qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) +#else + qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,rainwat) + q0(i,km1,ice_wat) +#endif elseif ( nwat==5 ) then qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + & q0(i,km1,snowwat) + q0(i,km1,rainwat) @@ -1118,6 +1177,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & enddo elseif ( nwat == 4 ) then do i=is,ie +#ifndef CCPP q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) #ifdef MULTI_GASES cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq @@ -1125,6 +1185,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, & #else cpm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq cvm(i) = (1.-(q0(i,kk,sphum)+q_liq))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq +#endif +#else + q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat) + q_sol = q0(i,kk,ice_wat) +#ifdef MULTI_GASES + cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#else + cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice + cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice +#endif + #endif enddo elseif ( nwat == 5 ) then From c92e0eee7d297f6050465672dcd6f9f88bd1ca58 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 6 Feb 2020 20:58:06 +0000 Subject: [PATCH 04/13] add dry mass fixer for IAU (requires change to GFS_typedefs.F90) --- driver/fvGFS/atmosphere.F90 | 75 +++++++++++++++++++++++++++++++++++-- 1 file changed, 72 insertions(+), 3 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 0d81b6cb1..49b2ac072 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -179,11 +179,11 @@ module atmosphere_mod use fv_fill_mod, only: fill_gfs use fv_dynamics_mod, only: fv_dynamics use fv_nesting_mod, only: twoway_nesting -use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height +use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height, prt_mass 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 @@ -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 @@ -1395,6 +1396,8 @@ 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, psuma, qsuma, betad + real, allocatable, dimension(:,:,:,:) :: qtmp Time_prev = Time Time_next = Time + Time_step_atmos rdt = 1.d0 / dt_atmos @@ -1407,6 +1410,32 @@ 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 + allocate ( qtmp(isd:ied ,jsd:jed ,npz, nwat) ) + ! 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,1) + do nb = 1,Atm_block%nblks + blen = Atm_block%blksz(nb) + do k=1,npz + 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) + qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat) + enddo + enddo + enddo + qsumb = g_sum(Atm(n)%domain,& + sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),& + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + endif + ! IAU increments are in units of 1/sec ! add analysis increment to u,v,t tendencies @@ -1447,6 +1476,35 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc enddo enddo enddo + + if (IAU_data%drymassfixer) then + ! global mean total pressure and water after IAU + psuma = 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,1) + do nb = 1,Atm_block%nblks + blen = Atm_block%blksz(nb) + do k=1,npz + 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) + qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat) + enddo + enddo + enddo + qsuma = g_sum(Atm(n)%domain,& + sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),& + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + ! rescale water consitituents to ensure IAU doesn't change dry mass + betad = (psuma - (psumb - qsumb))/qsuma + IPD_Data(nb)%Stateout%gq0 = betad*IPD_Data(nb)%Stateout%gq0 + deallocate(qtmp) + endif + endif call set_domain ( Atm(mytile)%domain ) @@ -1501,7 +1559,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc #ifdef MULTI_GASES q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:max(nwat,num_gas)))) + sum(qwat(1:max(nwat,num_gas))) #else - qt = sum(qwat(1:nwat)) + qt = sum(qwat(1:nwat)) q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:nwat))) + qt #endif Atm(n)%delp(i,j,k1) = q0 @@ -1575,6 +1633,17 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc call timing_off('FV_UPDATE_PHYS') call mpp_clock_end (id_dynam) + ! 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,1) + Atm(n)%ptop + !! global mean total water + !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,1) + !if( is_master() ) then + ! print *,'dry surface pressure after IAU/physics = ',psum-qsum + !endif + !--- nesting update after updating atmospheric variables with !--- physics tendencies if (ngrids > 1 .and. p_split > 0) then From 6f5f192d9822a9aed5bb745cc23e95966a296780 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Thu, 6 Feb 2020 23:58:14 +0000 Subject: [PATCH 05/13] update - still not quite working --- driver/fvGFS/atmosphere.F90 | 68 ++++++++++++++++++++++++------------- 1 file changed, 45 insertions(+), 23 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 49b2ac072..115c4e02b 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -1398,6 +1398,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc real(kind=kind_phys):: rcp, q0, qwat(nq), qt, rdt real psum, qsum, psumb, qsumb, psuma, qsuma, betad real, allocatable, dimension(:,:,:,:) :: qtmp + logical in_iau_interval,use_iau_massfixer Time_prev = Time Time_next = Time + Time_step_atmos rdt = 1.d0 / dt_atmos @@ -1406,9 +1407,18 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc nwat = Atm(n)%flagstruct%nwat dnats = Atm(mytile)%flagstruct%dnats nq_adv = nq - dnats + in_iau_interval=IAU_Data%in_interval + use_iau_massfixer=IAU_Data%drymassfixer if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers') +!SJL: perform vertical filling to fix the negative humidity if the SAS convection scheme is used +! This call may be commented out if RAS or other positivity-preserving CPS is used. +! do nb = 1,Atm_block%nblks +! blen = Atm_block%blksz(nb) +! call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys) +! enddo + if (IAU_Data%in_interval) then if (IAU_Data%drymassfixer) then @@ -1500,8 +1510,13 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),& isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) ! rescale water consitituents to ensure IAU doesn't change dry mass + ! following + ! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x betad = (psuma - (psumb - qsumb))/qsuma - IPD_Data(nb)%Stateout%gq0 = betad*IPD_Data(nb)%Stateout%gq0 + if (is_master()) print *,'betad = ',betad + do nb = 1,Atm_block%nblks + IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) = betad*IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) + enddo deallocate(qtmp) endif @@ -1517,7 +1532,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc #ifdef MULTI_GASES !$OMP num_gas, & #endif -!$OMP snowwat, graupel, nq_adv, flip_vc) & +!$OMP in_iau_interval, use_iau_massfixer,snowwat, graupel, nq_adv, flip_vc) & !$OMP private (nb, blen, i, j, k, k1, ix, q0, qwat, qt) do nb = 1,Atm_block%nblks @@ -1527,11 +1542,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) @@ -1551,6 +1566,9 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc q0 = IPD_Data(nb)%Statein%prsi(ix,k+1) - IPD_Data(nb)%Statein%prsi(ix,k) endif qwat(1:nq_adv) = q0*IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv) + if (in_iau_interval .and. use_iau_massfixer) then + Atm(n)%q(i,j,k1,1:nq_adv) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv) + else ! ********************************************************************************************************** ! Dry mass: the following way of updating delp is key to mass conservation with hybrid 32-64 bit computation ! ********************************************************************************************************** @@ -1565,6 +1583,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc Atm(n)%delp(i,j,k1) = q0 Atm(n)%q(i,j,k1,1:nq_adv) = qwat(1:nq_adv) / q0 ! if (dnats .gt. 0) Atm(n)%q(i,j,k1,nq_adv+1:nq) = IPD_Data(nb)%Stateout%gq0(ix,k,nq_adv+1:nq) + endif enddo enddo @@ -1573,11 +1592,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) @@ -1588,6 +1607,20 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc enddo ! nb-loop +! print out dry mass inside IAU interval. + if (IAU_Data%in_interval) 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,1) + Atm(n)%ptop + ! global mean total water + 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,1) + if( is_master() ) then + print *,'dry surface pressure in IAU interval = ',psum-qsum + endif + endif + call timing_off('GFS_TENDENCIES') w_diff = get_tracer_index (MODEL_ATMOS, 'w_diff' ) @@ -1633,17 +1666,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc call timing_off('FV_UPDATE_PHYS') call mpp_clock_end (id_dynam) - ! 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,1) + Atm(n)%ptop - !! global mean total water - !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,1) - !if( is_master() ) then - ! print *,'dry surface pressure after IAU/physics = ',psum-qsum - !endif - !--- nesting update after updating atmospheric variables with !--- physics tendencies if (ngrids > 1 .and. p_split > 0) then From fd00e4abfc7658649525bbc957bc9ef1f0f6debe Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 7 Feb 2020 01:47:38 +0000 Subject: [PATCH 06/13] update 2 - now working --- driver/fvGFS/atmosphere.F90 | 95 ++++++++++--------------------------- 1 file changed, 24 insertions(+), 71 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 115c4e02b..48a023125 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -1396,9 +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, psuma, qsuma, betad - real, allocatable, dimension(:,:,:,:) :: qtmp - logical in_iau_interval,use_iau_massfixer + real psum, qsum, psumb, qsumb, betad Time_prev = Time Time_next = Time + Time_step_atmos rdt = 1.d0 / dt_atmos @@ -1407,43 +1405,22 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc nwat = Atm(n)%flagstruct%nwat dnats = Atm(mytile)%flagstruct%dnats nq_adv = nq - dnats - in_iau_interval=IAU_Data%in_interval - use_iau_massfixer=IAU_Data%drymassfixer if( nq<3 ) call mpp_error(FATAL, 'GFS phys must have 3 interactive tracers') -!SJL: perform vertical filling to fix the negative humidity if the SAS convection scheme is used -! This call may be commented out if RAS or other positivity-preserving CPS is used. -! do nb = 1,Atm_block%nblks -! blen = Atm_block%blksz(nb) -! call fill_gfs(blen, npz, IPD_Data(nb)%Statein%prsi, IPD_Data(nb)%Stateout%gq0, 1.e-9_kind_phys) -! enddo if (IAU_Data%in_interval) then if (IAU_Data%drymassfixer) then - allocate ( qtmp(isd:ied ,jsd:jed ,npz, nwat) ) ! 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,1) - do nb = 1,Atm_block%nblks - blen = Atm_block%blksz(nb) - do k=1,npz - 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) - qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat) - enddo - enddo - enddo qsumb = g_sum(Atm(n)%domain,& - sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),& + 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,1) + if( is_master() ) then + print *,'dry surface pressure in IAU interval (before) = ',psumb++Atm%ptop-qsumb + endif endif ! IAU increments are in units of 1/sec @@ -1487,39 +1464,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc enddo enddo - if (IAU_data%drymassfixer) then - ! global mean total pressure and water after IAU - psuma = 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,1) - do nb = 1,Atm_block%nblks - blen = Atm_block%blksz(nb) - do k=1,npz - 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) - qtmp(i,j,k1,1:nwat) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nwat) - enddo - enddo - enddo - qsuma = g_sum(Atm(n)%domain,& - sum(Atm(n)%delp(isc:iec,jsc:jec,1:npz)*sum(qtmp(isc:iec,jsc:jec,1:npz,1:nwat),4),dim=3),& - isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) - ! rescale water consitituents to ensure IAU doesn't change dry mass - ! following - ! https://onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0870.2007.00299.x - betad = (psuma - (psumb - qsumb))/qsuma - if (is_master()) print *,'betad = ',betad - do nb = 1,Atm_block%nblks - IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) = betad*IPD_Data(nb)%Stateout%gq0(:,:,1:nwat) - enddo - deallocate(qtmp) - endif - endif call set_domain ( Atm(mytile)%domain ) @@ -1532,7 +1476,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc #ifdef MULTI_GASES !$OMP num_gas, & #endif -!$OMP in_iau_interval, use_iau_massfixer,snowwat, graupel, nq_adv, flip_vc) & +!$OMP snowwat, graupel, nq_adv, flip_vc) & !$OMP private (nb, blen, i, j, k, k1, ix, q0, qwat, qt) do nb = 1,Atm_block%nblks @@ -1566,9 +1510,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc q0 = IPD_Data(nb)%Statein%prsi(ix,k+1) - IPD_Data(nb)%Statein%prsi(ix,k) endif qwat(1:nq_adv) = q0*IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv) - if (in_iau_interval .and. use_iau_massfixer) then - Atm(n)%q(i,j,k1,1:nq_adv) = IPD_Data(nb)%Stateout%gq0(ix,k,1:nq_adv) - else ! ********************************************************************************************************** ! Dry mass: the following way of updating delp is key to mass conservation with hybrid 32-64 bit computation ! ********************************************************************************************************** @@ -1583,7 +1524,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc Atm(n)%delp(i,j,k1) = q0 Atm(n)%q(i,j,k1,1:nq_adv) = qwat(1:nq_adv) / q0 ! if (dnats .gt. 0) Atm(n)%q(i,j,k1,nq_adv+1:nq) = IPD_Data(nb)%Stateout%gq0(ix,k,nq_adv+1:nq) - endif enddo enddo @@ -1607,18 +1547,31 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc enddo ! nb-loop -! print out dry mass inside IAU interval. - if (IAU_Data%in_interval) then +! 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,1) + Atm(n)%ptop - ! global mean total water + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area,1) + ! 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,1) if( is_master() ) then - print *,'dry surface pressure in IAU interval = ',psum-qsum + print *,'dry surface pressure in IAU interval (after) =',& + psum+Atm%ptop-qsum endif + betad = (psum - (psumb - qsumb))/qsum + !if (is_master()) print *,'betad = ',betad + Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat) + ! global mean total water after 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,1) + !if( is_master() ) then + ! print *,'dry surface pressure in IAU interval (after 2) =',& + ! psum+Atm%ptop-qsum + !endif endif call timing_off('GFS_TENDENCIES') From a66f3a28c9cbaf8804247889bd9cd10c384ca023 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 7 Feb 2020 05:02:32 +0000 Subject: [PATCH 07/13] remove whitespace changes --- driver/fvGFS/atmosphere.F90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 48a023125..74c342d47 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -179,7 +179,7 @@ module atmosphere_mod use fv_fill_mod, only: fill_gfs use fv_dynamics_mod, only: fv_dynamics use fv_nesting_mod, only: twoway_nesting -use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height, prt_mass +use fv_diagnostics_mod, only: fv_diag_init, fv_diag, fv_time, prt_maxmin, prt_height 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 @@ -1408,9 +1408,7 @@ 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),& @@ -1463,7 +1461,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc enddo enddo enddo - endif call set_domain ( Atm(mytile)%domain ) @@ -1518,7 +1515,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc #ifdef MULTI_GASES q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:max(nwat,num_gas)))) + sum(qwat(1:max(nwat,num_gas))) #else - qt = sum(qwat(1:nwat)) + qt = sum(qwat(1:nwat)) q0 = Atm(n)%delp(i,j,k1)*(1.-sum(Atm(n)%q(i,j,k1,1:nwat))) + qt #endif Atm(n)%delp(i,j,k1) = q0 From d4c23401834035bb3aaf64ce8b3c646e0aa819a7 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 7 Feb 2020 13:52:11 +0000 Subject: [PATCH 08/13] remove debug prints --- driver/fvGFS/atmosphere.F90 | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 74c342d47..bfef8f290 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -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, is_master +use fv_mp_mod, only: switch_current_Atm 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 @@ -1416,9 +1416,6 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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,1) - if( is_master() ) then - print *,'dry surface pressure in IAU interval (before) = ',psumb++Atm%ptop-qsumb - endif endif ! IAU increments are in units of 1/sec @@ -1554,21 +1551,8 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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,1) - if( is_master() ) then - print *,'dry surface pressure in IAU interval (after) =',& - psum+Atm%ptop-qsum - endif betad = (psum - (psumb - qsumb))/qsum - !if (is_master()) print *,'betad = ',betad Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat) - ! global mean total water after 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,1) - !if( is_master() ) then - ! print *,'dry surface pressure in IAU interval (after 2) =',& - ! psum+Atm%ptop-qsum - !endif endif call timing_off('GFS_TENDENCIES') From 7a03030c627db5654d4bfa2a7ac5957e41643629 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Fri, 7 Feb 2020 13:56:47 +0000 Subject: [PATCH 09/13] add drymassfixer to IAU_data structure --- tools/fv_iau_mod.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index e9b961cbb..bf2a73ea1 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -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 @@ -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 From 15abb3d9030a66078f2f07c14f30eac1898c8205 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sat, 8 Feb 2020 14:58:57 +0000 Subject: [PATCH 10/13] use area_64 in g_sum call --- driver/fvGFS/atmosphere.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index bfef8f290..685146e07 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -1412,10 +1412,10 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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,1) + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) 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,1) + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) endif ! IAU increments are in units of 1/sec @@ -1546,11 +1546,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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,1) + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) ! 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,1) + isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) betad = (psum - (psumb - qsumb))/qsum Atm(n)%q(:,:,:,1:nwat) = betad*Atm(n)%q(:,:,:,1:nwat) endif From 0c3016190c2851d3db54a7878279ec2869579f1b Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sat, 8 Feb 2020 19:37:31 +0000 Subject: [PATCH 11/13] add debug prints --- driver/fvGFS/atmosphere.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 685146e07..a044895dc 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -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 @@ -1416,6 +1416,9 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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) + if (is_master()) then + print *,'dry ps before IAU',psumb+Atm(n)%ptop-qsumb + endif endif ! IAU increments are in units of 1/sec @@ -1552,7 +1555,16 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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) 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') From 0e30846b5ddacdc8585a288bf8d38dfd30df6a12 Mon Sep 17 00:00:00 2001 From: jswhit2 Date: Sun, 9 Feb 2020 03:09:33 +0000 Subject: [PATCH 12/13] remove one debug print --- driver/fvGFS/atmosphere.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index a044895dc..39ba8c4c8 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -1417,7 +1417,7 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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 before IAU',psumb+Atm(n)%ptop-qsumb + print *,'dry ps before IAU/physics',psumb+Atm(n)%ptop-qsumb endif endif @@ -1556,15 +1556,15 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc isc,iec,jsc,jec,Atm(n)%ng,Atm(n)%gridstruct%area_64,1) betad = (psum - (psumb - qsumb))/qsum if (is_master()) then - print *,'dry ps after IAU+physics',psum+Atm(n)%ptop-qsum + 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 + !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') From da9077d09a6e294d5b383710c6c5fdddecff46e6 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Mon, 17 Feb 2020 21:36:06 +0000 Subject: [PATCH 13/13] add reproduce in g_sum call --- driver/fvGFS/atmosphere.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index 39ba8c4c8..0fc5ba08e 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -1412,10 +1412,10 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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) + 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) + 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 @@ -1549,11 +1549,11 @@ subroutine atmosphere_state_update (Time, IPD_Data, IAU_Data, Atm_block, flip_vc 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) + 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) + 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