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.

Sign up for GitHub

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

FA in fv_regional_bc and external_ic #17

Merged
merged 3 commits into from
Aug 6, 2020
Merged

FA in fv_regional_bc and external_ic #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5dd06bb
- fix the ntrac> ntracers limitation in external_ic.
mzhangw Jul 31, 2020
1951678
add negative hydrometeors adjustment for F-A scheme
mzhangw Aug 1, 2020
1ae0439
- add temperature rebalance in the last step of lagrangian_to_Euleri…
mzhangw Aug 5, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 40 additions & 7 deletions model/fv_dynamics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -275,6 +275,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
integer :: i,j,k, n, iq, n_map, nq, nr, nwat, k_split
integer :: sphum, liq_wat = -999, ice_wat = -999 ! GFDL physics
integer :: rainwat = -999, snowwat = -999, graupel = -999, cld_amt = -999
integer :: q_rimef = -999 ! FA mp
integer :: theta_d = -999
#ifdef CCPP
logical used, do_omega
Expand Down Expand Up @@ -408,6 +409,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
ice_wat = get_tracer_index (MODEL_ATMOS, 'ice_wat')
rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat')
snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat')
q_rimef = get_tracer_index (MODEL_ATMOS, 'q_rimef')
graupel = get_tracer_index (MODEL_ATMOS, 'graupel')
cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt')
endif
Expand All @@ -434,7 +436,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
#else
!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,nwat,q,q_con,sphum,liq_wat, &
#endif
!$OMP rainwat,ice_wat,snowwat,graupel) private(cvm,i,j,k)
!$OMP rainwat,ice_wat,snowwat,q_rimef,graupel) private(cvm,i,j,k)
do k=1,npz
do j=js,je
#ifdef USE_COND
Expand All @@ -452,7 +454,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
#else
!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,q,q_con,sphum,liq_wat, &
#endif
!$OMP rainwat,ice_wat,snowwat,graupel,pkz,flagstruct, &
!$OMP rainwat,ice_wat,snowwat,graupel,q_rimef,pkz,flagstruct, &
#ifdef MULTI_GASES
!$OMP kapad, &
#endif
Expand Down Expand Up @@ -730,6 +732,8 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,snowwat), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)
if ( graupel > 0 ) &
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,graupel), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)
if ( q_rimef > 0 ) &
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,q_rimef), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)
call timing_off('Fill2D')
endif
#endif
Expand Down Expand Up @@ -775,11 +779,18 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
if (is_master()) write(*,'(A, I3, A1, I3)') 'finished k_split ', n_map, '/', k_split
call prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( liq_wat > 0 ) &
call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( rainwat > 0 ) &
call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( ice_wat > 0 ) &
call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( snowwat > 0 ) &
call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( graupel > 0 ) &
call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if ( q_rimef > 0 ) &
call prt_mxm('q_rimef_dyn', q(isd,jsd,1,q_rimef), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
endif
#ifdef AVEC_TIMERS
call avec_timer_stop(6)
Expand Down Expand Up @@ -909,6 +920,28 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
endif
endif


!mzhang:F-A nwat=4
if( nwat == 4 ) then
call neg_adj2(is, ie, js, je, ng, npz, &
flagstruct%hydrostatic, &
peln, delz, &
pt, delp, q(isd,jsd,1,sphum), &
q(isd,jsd,1,liq_wat), &
q(isd,jsd,1,rainwat), &
q(isd,jsd,1,ice_wat), &
check_negative=flagstruct%check_negative)
if ( flagstruct%fv_debug ) then
call prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
endif
endif

!mzhang

if( (flagstruct%consv_am.or.idiag%id_amdt>0.or.idiag%id_aam>0) .and. (.not.do_adiabatic_init) ) then
call compute_aam(npz, is, ie, js, je, isd, ied, jsd, jed, gridstruct, bd, &
ptop, ua, va, u, v, delp, te_2d, ps, m_fac)
Expand Down
33 changes: 17 additions & 16 deletions model/fv_mapz.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,7 +232,8 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
integer:: i,j,k
integer:: kdelz
#ifdef CCPP
integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kp, k_next
integer :: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kp, k_next
integer :: q_rimef ! HWRF
integer :: ierr
#else
integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kmp, kp, k_next
Expand All @@ -255,6 +256,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
graupel = get_tracer_index (MODEL_ATMOS, 'graupel')
cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt')
ccn_cm3 = get_tracer_index (MODEL_ATMOS, 'ccn_cm3')
q_rimef = get_tracer_index (MODEL_ATMOS, 'q_rimef')

if ( do_adiabatic_init .or. do_sat_adj ) then
fast_mp_consv = (.not.do_adiabatic_init) .and. consv>consv_min
Expand Down Expand Up @@ -630,7 +632,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
#if defined(CCPP) && defined(__GFORTRAN__)
!$OMP parallel default(none) shared(is,ie,js,je,km,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, &
!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat,q_rimef, &
!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
Expand All @@ -645,7 +647,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
#elif defined(CCPP)
!$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, &
!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, q_rimef, &
!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
Expand All @@ -661,7 +663,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
#else
!$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, &
!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln,adiabatic, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat,q_rimef, &
!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
Expand Down Expand Up @@ -896,6 +898,17 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / ((1.+r_vir*qv(i))*(1.-gz(i)))
#endif
enddo
!HWRF
elseif ( nwat==4 ) then
do i=is,ie
gz(i) = q(i,j,k,liq_wat)+q(i,j,k,rainwat)+q(i,j,k,ice_wat)
#ifdef MULTI_GASES
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/ virq(q(i,j,k,1:num_gas))
#else
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i)))
#endif
enddo
!HWRF
elseif ( nwat==6 ) then
do i=is,ie
gz(i) = q(i,j,k,liq_wat)+q(i,j,k,rainwat)+q(i,j,k,ice_wat)+q(i,j,k,snowwat)+q(i,j,k,graupel)
Expand Down Expand Up @@ -3595,15 +3608,6 @@ 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)
Expand All @@ -3612,9 +3616,6 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai
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)
Expand Down
83 changes: 83 additions & 0 deletions model/fv_regional_bc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,7 @@ module fv_regional_mod
,o3mr_index & ! in the tracers
,rain_water_index & ! array.
,snow_water_index & !
,q_rimef_index & ! F-A MP
,sphum_index !<--
!
integer,save :: lbnd_x_tracers,lbnd_y_tracers & !<-- Local lower bounds of x,y for tracer arrays
Expand Down Expand Up @@ -741,6 +742,7 @@ subroutine setup_regional_BC(Atm, dt_atmos &
ice_water_index = get_tracer_index(MODEL_ATMOS, 'ice_wat')
rain_water_index = get_tracer_index(MODEL_ATMOS, 'rainwat')
snow_water_index = get_tracer_index(MODEL_ATMOS, 'snowwat')
q_rimef_index = get_tracer_index(MODEL_ATMOS, 'q_rimef')
graupel_index = get_tracer_index(MODEL_ATMOS, 'graupel')
cld_amt_index = get_tracer_index(MODEL_ATMOS, 'cld_amt')
o3mr_index = get_tracer_index(MODEL_ATMOS, 'o3mr')
Expand Down Expand Up @@ -3437,6 +3439,8 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
!!! High-precision
integer i,ie,is,j,je,js,k,l,m, k2,iq
integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
!mzhang
integer q_rimef
!
!---------------------------------------------------------------------------------
!
Expand All @@ -3445,6 +3449,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
ice_wat = ice_water_index
rainwat = rain_water_index
snowwat = snow_water_index
q_rimef = q_rimef_index
graupel = graupel_index
cld_amt = cld_amt_index
o3mr = o3mr_index
Expand Down Expand Up @@ -3745,6 +3750,51 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
enddo
enddo
endif
!zhang nwat=4
if ( Atm%flagstruct%nwat .eq. 4 ) then
do k=1,npz
do i=is,ie
qn1(i,k) = BC_side%q_BC(i,j,k,liq_wat)
BC_side%q_BC(i,j,k,rainwat) = 0.
BC_side%q_BC(i,j,k,q_rimef) = 0.
BC_side%q_BC(i,j,k,snowwat) = 0.
BC_side%q_BC(i,j,k,graupel) = 0.
if ( BC_side%pt_BC(i,j,k) > 273.16 ) then ! > 0C all liq_wat
BC_side%q_BC(i,j,k,liq_wat) = qn1(i,k)
BC_side%q_BC(i,j,k,ice_wat) = 0.
#ifdef ORIG_CLOUDS_PART
else if ( BC_side%pt_BC(i,j,k) < 258.16 ) then ! < -15C all ice_wat
BC_side%q_BC(i,j,k,liq_wat) = 0.
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
else ! between -15~0C: linear interpolation
BC_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((BC_side%pt_BC(i,j,k)-258.16)/15.)
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - BC_side%q_BC(i,j,k,liq_wat)
endif
#else
else if ( BC_side%pt_BC(i,j,k) < 233.16 ) then ! < -40C all ice_wat
BC_side%q_BC(i,j,k,liq_wat) = 0.
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
else
if ( k.eq.1 ) then ! between [-40,0]: linear interpolation
BC_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((BC_side%pt_BC(i,j,k)-233.16)/40.)
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - BC_side%q_BC(i,j,k,liq_wat)
else
if (BC_side%pt_BC(i,j,k)<258.16 .and. BC_side%q_BC(i,j,k-1,ice_wat)>1.e-5 ) then
BC_side%q_BC(i,j,k,liq_wat) = 0.
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
else ! between [-40,0]: linear interpolation
BC_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((BC_side%pt_BC(i,j,k)-233.16)/40.)
BC_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - BC_side%q_BC(i,j,k,liq_wat)
endif
endif
endif
#endif
call mp_auto_conversion_fa(BC_side%q_BC(i,j,k,liq_wat), BC_side%q_BC(i,j,k,rainwat), &
BC_side%q_BC(i,j,k,ice_wat), BC_side%q_BC(i,j,k,snowwat) )
enddo
enddo
endif
!zhang
endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
!
! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
Expand Down Expand Up @@ -5328,6 +5378,26 @@ subroutine mp_auto_conversion(ql, qr, qi, qs)

end subroutine mp_auto_conversion

!zhang the version for F-A
subroutine mp_auto_conversion_fa(ql, qr, qi, qs)
real, intent(inout):: ql, qr, qi, qs
real, parameter:: qi0_max = 2.0e-3
real, parameter:: ql0_max = 2.5e-3

! Convert excess cloud water into rain:
if ( ql > ql0_max ) then
qr = ql - ql0_max
ql = ql0_max
endif
! Convert excess cloud ice into snow:
! if ( qi > qi0_max ) then
! qs = qi - qi0_max
! qi = qi0_max
! endif

end subroutine mp_auto_conversion_fa


!-----------------------------------------------------------------------
!
subroutine nudge_qv_bc(Atm,isd,ied,jsd,jed)
Expand Down Expand Up @@ -6657,6 +6727,8 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
!
integer :: k, j, i, iq, is, ie, js, je
integer :: liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
!mzhang
integer :: q_rimef
real :: qt, wt, m_fac

is=lbound(BC_side%delp_BC,1)
Expand All @@ -6668,6 +6740,7 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
ice_wat = get_tracer_index(MODEL_ATMOS, 'ice_wat')
rainwat = get_tracer_index(MODEL_ATMOS, 'rainwat')
snowwat = get_tracer_index(MODEL_ATMOS, 'snowwat')
q_rimef = get_tracer_index(MODEL_ATMOS, 'q_rimef')
graupel = get_tracer_index(MODEL_ATMOS, 'graupel')
cld_amt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
!
Expand All @@ -6684,6 +6757,11 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
BC_side%q_BC(i,j,k,rainwat) + &
BC_side%q_BC(i,j,k,snowwat) + &
BC_side%q_BC(i,j,k,graupel))
!zhang nwat =4
else if ( nwat == 4 ) then
qt = wt*(1. + BC_side%q_BC(i,j,k,liq_wat) + &
BC_side%q_BC(i,j,k,ice_wat) + &
BC_side%q_BC(i,j,k,rainwat))
else ! all other values of nwat
qt = wt*(1. + sum(BC_side%q_BC(i,j,k,2:nwat)))
endif
Expand All @@ -6706,6 +6784,11 @@ subroutine set_delp_and_tracers(BC_side,npz,nwat)
BC_side%q_BC(i,j,k,rainwat) + &
BC_side%q_BC(i,j,k,snowwat) + &
BC_side%q_BC(i,j,k,graupel))
!zhang nwat=4
else if ( nwat == 4 ) then
qt = wt*(1. + BC_side%q_BC(i,j,k,liq_wat) + &
BC_side%q_BC(i,j,k,ice_wat) + &
BC_side%q_BC(i,j,k,rainwat))
else ! all other values of nwat
qt = wt*(1. + sum(BC_side%q_BC(i,j,k,2:nwat)))
endif
Expand Down
Loading