Skip to content

Commit

Permalink
Merge pull request #17 from mzhangw/hafs_develop_mzhang
Browse files Browse the repository at this point in the history 
FA in fv_regional_bc and external_ic
  • Loading branch information
@climbfuji
climbfuji authored Aug 6, 2020
2 parents b85c7ec + 1ae0439 commit fd80f7d
Show file tree
Hide file tree
Showing 5 changed files with 188 additions and 38 deletions.
47 changes: 40 additions & 7 deletions model/fv_dynamics.F90
View file
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
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
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

0 comments on commit fd80f7d

Please sign in to comment.