diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 16c0e5753..3296a32e9 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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) diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90 index f4c507fa8..1cf0829fa 100644 --- a/model/fv_mapz.F90 +++ b/model/fv_mapz.F90 @@ -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 @@ -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 @@ -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, & @@ -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, & @@ -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, & @@ -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) @@ -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) @@ -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) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 98285b7a3..9d06033bd 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -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 @@ -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') @@ -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 ! !--------------------------------------------------------------------------------- ! @@ -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 @@ -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 @@ -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) @@ -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) @@ -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') ! @@ -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 @@ -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 diff --git a/model/fv_sg.F90 b/model/fv_sg.F90 index 992ad22f9..8aad05e04 100644 --- a/model/fv_sg.F90 +++ b/model/fv_sg.F90 @@ -1941,12 +1941,14 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & integer, intent(in):: is, ie, js, je, ng, kbot logical, intent(in):: hydrostatic real, intent(in):: dp(is-ng:ie+ng,js-ng:je+ng,kbot) !< total delp-p - real, intent(in):: delz(is-ng:,js-ng:,1:) + !mzhang: bug + !real, intent(in):: delz(is-ng:,js-ng:,1:) + real, intent(in):: delz(is:,js:,1:) real, intent(in):: peln(is:ie,kbot+1,js:je) !< ln(pe) logical, intent(in), OPTIONAL :: check_negative real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: & - pt, qv, ql, qr, qi, qs - real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa + pt, qv, ql, qr, qi + real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa,qs ! Local: logical:: sat_adj = .false. real, parameter :: t48 = tice - 48. @@ -1965,7 +1967,9 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & call prt_negative('liq_wat', ql, is, ie, js, je, ng, kbot, -1.e-7) call prt_negative('rainwat', qr, is, ie, js, je, ng, kbot, -1.e-7) call prt_negative('ice_wat', qi, is, ie, js, je, ng, kbot, -1.e-7) - call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7) + if (present(qs)) then + call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7) + endif endif endif @@ -1986,7 +1990,9 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & qv2(i,j) = qv(i,j,k) ql2(i,j) = ql(i,j,k) qi2(i,j) = qi(i,j,k) - qs2(i,j) = qs(i,j,k) + if (present(qs)) then + qs2(i,j) = qs(i,j,k) + endif qr2(i,j) = qr(i,j,k) dp2(i,j) = dp(i,j,k) pt2(i,j) = pt(i,j,k) @@ -1998,7 +2004,11 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & do i=is, ie p2(i,j) = dp2(i,j)/(peln(i,k+1,j)-peln(i,k,j)) q_liq = max(0., ql2(i,j) + qr2(i,j)) - q_sol = max(0., qi2(i,j) + qs2(i,j)) + if (present(qs)) then + q_sol = max(0., qi2(i,j) + qs2(i,j)) + else + q_sol = max(0., qi2(i,j) ) + endif oneocpm = 1.0 / ((1.-(qv2(i,j)+q_liq+q_sol))*cp_air + qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice) lcpk(i,j) = hlv * oneocpm icpk(i,j) = hlf * oneocpm @@ -2009,7 +2019,11 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & do i=is, ie p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*(1.+zvir*qv2(i,j)) q_liq = max(0., ql2(i,j) + qr2(i,j)) - q_sol = max(0., qi2(i,j) + qs2(i,j)) + if (present(qs)) then + q_sol = max(0., qi2(i,j) + qs2(i,j)) + else + q_sol = max(0., qi2(i,j) ) + endif oneocpm = 1.0 / ((1.-(qv2(i,j)+q_liq+q_sol))*cv_air + qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice) lcpk(i,j) = (lv00+d0_vap*pt2(i,j)) * oneocpm icpk(i,j) = (Li0+dc_ice*pt2(i,j)) * oneocpm @@ -2023,18 +2037,26 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & !----------- do j=js, je do i=is, ie - qsum = qi2(i,j) + qs2(i,j) + if (present(qs)) then + qsum = qi2(i,j) + qs2(i,j) + else + qsum = qi2(i,j) + endif if ( qsum > 0. ) then if ( qi2(i,j) < 0. ) then qi2(i,j) = 0. - qs2(i,j) = qsum + if(present(qs)) then + qs2(i,j) = qsum + endif elseif ( qs2(i,j) < 0. ) then qs2(i,j) = 0. qi2(i,j) = qsum endif else qi2(i,j) = 0. + if(present(qs)) then qs2(i,j) = qsum + endif ! If qsum is negative then borrow from rain water: phase change if ( qs2(i,j) < 0. .and. qr2(i,j) > 0. ) then @@ -2076,10 +2098,18 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & qr2(i,j) = qsum ! rain water is still negative if ( qr(i,j,k) < 0. ) then ! fill negative rain with available qi & qs (cooling) - dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) ) + if (present(qs)) then + dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) ) + else + dq = min( qi2(i,j), -qr2(i,j) ) + endif qr2(i,j) = qr2(i,j) + dq - dq1 = min( dq, qs2(i,j) ) - qs2(i,j) = qs2(i,j) - dq1 + if (present(qs)) then + dq1 = min( dq, qs2(i,j) ) + qs2(i,j) = qs2(i,j) - dq1 + else + dq1 = dq + endif qi2(i,j) = qi2(i,j) + dq1 - dq pt2(i,j) = pt2(i,j) - dq*icpk(i,j) endif @@ -2138,7 +2168,9 @@ subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, & qv(i,j,k) = qv2(i,j) ql(i,j,k) = ql2(i,j) qi(i,j,k) = qi2(i,j) - qs(i,j,k) = qs2(i,j) + if (present(qs)) then + qs(i,j,k) = qs2(i,j) + endif qr(i,j,k) = qr2(i,j) pt(i,j,k) = pt2(i,j) enddo diff --git a/tools/external_ic.F90 b/tools/external_ic.F90 index 432703501..15e686548 100644 --- a/tools/external_ic.F90 +++ b/tools/external_ic.F90 @@ -570,8 +570,9 @@ subroutine get_nggps_ic (Atm, fv_domain, dt_atmos ) !--- read in the number of tracers in the NCEP NGGPS ICs call read_data ('INPUT/'//trim(fn_gfs_ctl), 'ntrac', ntrac, no_domain=.TRUE.) - if (ntrac > ntracers) call mpp_error(FATAL,'==> External_ic::get_nggps_ic: more NGGPS tracers & - &than defined in field_table '//trim(fn_gfs_ctl)//' for NGGPS IC') +!mzhang: not in FA +! if (ntrac > ntracers) call mpp_error(FATAL,'==> External_ic::get_nggps_ic: more NGGPS tracers & +! &than defined in field_table '//trim(fn_gfs_ctl)//' for NGGPS IC') ! call get_data_source(source,Atm%flagstruct%regional)