Skip to content

Commit

Permalink
feature/chgres_cube_grib2_release: This commit references #7.
Browse files Browse the repository at this point in the history 
Remove more code for GSD physics, which will not be supported
for the next public release.
  • Loading branch information
@GeorgeGayno-NOAA
GeorgeGayno-NOAA committed Jan 14, 2020
1 parent ec372ba commit 999da97
Show file tree
Hide file tree
Showing 4 changed files with 1 addition and 316 deletions.
303 changes: 1 addition & 302 deletions sorc/chgres_cube.fd/atmosphere.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,7 @@ module atmosphere
pres_input_grid, &
terrain_input_grid, &
read_input_atm_data, &
cleanup_input_atm_data, &
P_QC, P_QNC, P_QI, P_QV, &
P_QNI, P_QR, P_QNR, P_QNWFA
cleanup_input_atm_data

use model_grid, only : target_grid, &
latitude_s_target_grid, &
Expand Down Expand Up @@ -392,21 +390,6 @@ subroutine atmosphere_driver(localpet)
call convert_winds

!-----------------------------------------------------------------------------------
!.. Let us ensure that double-moment microphysics variables have numbers
!.. where there is mass. Currently doing this for Thompson-MP only, but
!.. can consider doing it for every MP scheme that has 2-moment variables.
!.. This is important because pressure-level RAP/HRRR files have mass but
!.. not number values for example (whereas native model level files have
!.. both).
!-----------------------------------------------------------------------------------

if ((trim(phys_suite)=="RAP" .or. trim(phys_suite)=="GSD") .and. &
trim(input_type) == "grib2") then
call create_number_concentrations
endif

!call ESMF_FieldDestroy(landmask_target_grid, rc=rc)
!-----------------------------------------------------------------------------------
! Write target data to file.
!-----------------------------------------------------------------------------------

Expand Down Expand Up @@ -512,14 +495,6 @@ subroutine create_atm_esmf_fields

do n = 1, num_tracers
print*,"- CALL FieldCreate FOR TARGET GRID TRACERS ", trim(tracers(n))
if (trim(tracers(n)) == "sphum") P_QV = n
if (trim(tracers(n)) == "liq_wat") P_QC = n
if (trim(tracers(n)) == "ice_wat") P_QI = n
if (trim(tracers(n)) == "rainwat") P_QR = n
if (trim(tracers(n)) == "ice_nc") P_QNI = n
if (trim(tracers(n)) == "rain_nc") P_QNR = n
if (trim(tracers(n)) == "water_nc") P_QNC = n
if (trim(tracers(n)) == "liq_aero") P_QNWFA = n
tracers_target_grid(n) = ESMF_FieldCreate(target_grid, &
typekind=ESMF_TYPEKIND_R8, &
staggerloc=ESMF_STAGGERLOC_CENTER, &
Expand Down Expand Up @@ -1759,282 +1734,6 @@ subroutine compute_zh

end subroutine compute_zh

subroutine create_number_concentrations

use esmf
implicit none

real(esmf_kind_r8), pointer :: tempptr(:,:,:), presptr(:,:,:), &
cloudptr(:,:,:), vaporptr(:,:,:), &
iceptr(:,:,:), rainptr(:,:,:), &
qncptr(:,:,:), qniptr(:,:,:), &
qnrptr(:,:,:), qnwfaptr(:,:,:)
integer(esmf_kind_i8), pointer :: landptr(:,:)
real(esmf_kind_r8) :: alt, temp_rho
integer :: i,j,k, clb(3), cub(3), rc
real(esmf_kind_r8), parameter :: cvpm = -0.714285731, &
Rd = 287.05, &
Rv =461.51

call ESMF_FieldGet(temp_target_grid, &
computationalLBound=clb, &
computationalUBound=cub, &
farrayPtr=tempptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet temp", rc)

call ESMF_FieldGet(pres_target_grid, &
farrayPtr=presptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet pres", rc)

call ESMF_FieldGet(landmask_target_grid, &
farrayPtr=landptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet landmask", rc)

call ESMF_FieldGet(tracers_target_grid(P_QV), farrayPtr=vaporptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qc", rc)

call ESMF_FieldGet(tracers_target_grid(P_QC), farrayPtr=cloudptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qc", rc)

if (P_QI > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QI), farrayPtr=iceptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qi", rc)
endif

if (P_QR > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QR), farrayPtr=rainptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qr", rc)
endif

if (P_QNC > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QNC), farrayPtr=qncptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qnc", rc)
endif

if (P_QNI > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QNI), farrayPtr=qniptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qni", rc)
endif

if (P_QNR > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QNR), farrayPtr=qnrptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qnr", rc)
endif

if (P_QNWFA > 0) then
call ESMF_FieldGet(tracers_target_grid(P_QNWFA), farrayPtr=qnwfaptr, rc=rc)
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
call error_handler("IN FieldGet tracers_qnr", rc)
endif



do i = clb(1),cub(1)
do j = clb(2),cub(2)
do k = clb(3),cub(3)

alt = Rd/presptr(i,j,k)*tempptr(i,j,k)*(1+Rv/Rd*vaporptr(i,j,k))
temp_rho = 1./alt

!..Produce a sensible cloud droplet number concentration

if (P_QNC.gt.1 .AND. cloudptr(i,j,k).gt.0.0 .AND. qncptr(i,j,k).le.0.0) then
if (P_QNWFA .gt. 1) then
qncptr(i,j,k) = make_DropletNumber (vaporptr(i,j,k)*temp_rho, &
& qnwfaptr(i,j,k)*temp_rho, real(landptr(i,j)))
else
qncptr(i,j,k) = make_DropletNumber (vaporptr(i,j,k)*temp_rho, &
& 0.0, real(landptr(i,j)))
endif
qncptr(i,j,k) = qncptr(i,j,k) / temp_rho
endif

!..Produce a sensible cloud ice number concentration

if (P_QNI.gt.1 .AND. iceptr(i,j,k).gt.0.0 .AND. qniptr(i,j,k).le.0.0) then
qniptr(i,j,k) = make_IceNumber (iceptr(i,j,k)*temp_rho, tempptr(i,j,k))
qniptr(i,j,k) = qniptr(i,j,k) / temp_rho
endif

!..Produce a sensible rain number concentration

if (P_QNR.gt.1 .AND. rainptr(i,j,k).gt.0.0 .AND. qnrptr(i,j,k).le.0.0) then
qnrptr(i,j,k) = make_RainNumber (rainptr(i,j,k)*temp_rho, tempptr(i,j,k))
qnrptr(i,j,k) = qnrptr(i,j,k) / temp_rho
endif

enddo

enddo
enddo

end subroutine create_number_concentrations

!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+

real function make_IceNumber (Q_ice, temp)

IMPLICIT NONE
REAL, PARAMETER:: Ice_density = 890.0
REAL, PARAMETER:: PI = 3.1415926536
integer idx_rei
real corr, reice, deice, Q_ice, temp
double precision lambda

!+---+-----------------------------------------------------------------+
!..Table of lookup values of radiative effective radius of ice crystals
!.. as a function of Temperature from -94C to 0C. Taken from WRF RRTMG
!.. radiation code where it is attributed to Jon Egill Kristjansson
!.. and coauthors.
!+---+-----------------------------------------------------------------+

real retab(95)
data retab / &
5.92779, 6.26422, 6.61973, 6.99539, 7.39234, &
7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, &
10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, &
15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, &
20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, &
27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, &
31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, &
34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, &
38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, &
42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, &
50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, &
65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, &
93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, &
124.954, 130.630, 136.457, 142.446, 148.608, 154.956, &
161.503, 168.262, 175.248, 182.473, 189.952, 197.699, &
205.728, 214.055, 222.694, 231.661, 240.971, 250.639/

!+---+-----------------------------------------------------------------+
!..From the model 3D temperature field, subtract 179K for which
!.. index value of retab as a start. Value of corr is for
!.. interpolating between neighboring values in the table.
!+---+-----------------------------------------------------------------+

idx_rei = int(temp-179.)
idx_rei = min(max(idx_rei,1),94)
corr = temp - int(temp)
reice = retab(idx_rei)*(1.-corr) + retab(idx_rei+1)*corr
deice = 2.*reice * 1.E-6

!+---+-----------------------------------------------------------------+
!..Now we have the final radiative effective size of ice (as function
!.. of temperature only). This size represents 3rd moment divided by
!.. second moment of the ice size distribution, so we can compute a
!.. number concentration from the mean size and mass mixing ratio.
!.. The mean (radiative effective) diameter is 3./Slope for an inverse
!.. exponential size distribution. So, starting with slope, work
!.. backwords to get number concentration.
!+---+-----------------------------------------------------------------+

lambda = 3.0 / deice
make_IceNumber = Q_ice * lambda*lambda*lambda / (PI*Ice_density)

!+---+-----------------------------------------------------------------+
!..Example1: Common ice size coming from Thompson scheme is about 30 microns.
!.. An example ice mixing ratio could be 0.001 g/kg for a temperature of -50C.
!.. Remember to convert both into MKS units. This gives N_ice=357652 per kg.
!..Example2: Lower in atmosphere at T=-10C matching ~162 microns in retab,
!.. and assuming we have 0.1 g/kg mixing ratio, then N_ice=28122 per kg,
!.. which is 28 crystals per liter of air if the air density is 1.0.
!+---+-----------------------------------------------------------------+

return
end function make_IceNumber

!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+

real function make_DropletNumber (Q_cloud, qnwfa, xland)

IMPLICIT NONE

real:: Q_cloud, qnwfa, xland

real, parameter:: PI = 3.1415926536
real, parameter:: am_r = PI*1000./6.
real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, &
& 504,720,990,1320,1716,2184,2730,3360,4080,4896/)
double precision:: lambda, qnc
real:: q_nwfa, x1, xDc
integer:: nu_c

!+---+

if (qnwfa .le. 0.0) then

if ((xland-1.5).gt.0.) then !--- Ocean
xDc = 17.E-6
nu_c = 12
else !--- Land
xDc = 11.E-6
nu_c = 4
endif

else
q_nwfa = MAX(99.E6, MIN(qnwfa,5.E10))
nu_c = MAX(2, MIN(NINT(2.5E10/q_nwfa), 15))

x1 = MAX(1., MIN(q_nwfa*1.E-9, 10.)) - 1.
xDc = (30. - x1*20./9.) * 1.E-6
endif

lambda = (4.0D0 + nu_c) / xDc
qnc = Q_cloud / g_ratio(nu_c) * lambda*lambda*lambda / am_r
make_DropletNumber = SNGL(qnc)

return
end function make_DropletNumber

!+---+-----------------------------------------------------------------+
!+---+-----------------------------------------------------------------+

real function make_RainNumber (Q_rain, temp)

IMPLICIT NONE

real, intent(in):: Q_rain, temp
double precision:: lambda, N0, qnr
real, parameter:: PI = 3.1415926536
real, parameter:: am_r = PI*1000./6.

!+---+-----------------------------------------------------------------+
!.. Not thrilled with it, but set Y-intercept parameter to Marshal-Palmer value
!.. that basically assumes melting snow becomes typical rain. However, for
!.. -2C < T < 0C, make linear increase in exponent to attempt to keep
!.. supercooled collision-coalescence (warm-rain) similar to drizzle rather
!.. than bigger rain drops. While this could also exist at T>0C, it is
!.. more difficult to assume it directly from having mass and not number.
!+---+-----------------------------------------------------------------+

N0 = 8.E6

if (temp .le. 271.15) then
N0 = 8.E8
elseif (temp .gt. 271.15 .and. temp.lt.273.15) then
N0 = 8. * 10**(279.15-temp)
endif

lambda = SQRT(SQRT(N0*am_r*6.0/Q_rain))
qnr = Q_rain / 6.0 * lambda*lambda*lambda / am_r
make_RainNumber = SNGL(qnr)

return
end function make_RainNumber

subroutine cleanup_target_atm_b4adj_data

implicit none
Expand Down
4 changes: 0 additions & 4 deletions sorc/chgres_cube.fd/input_data.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,10 +102,6 @@ module input_data

integer, public :: lsoil_input=4 ! # of soil layers,
! # hardwire for now

integer, public :: P_QC, P_QNC,P_QI, & ! Use for keeping track of tracer
P_QNI, P_QR, P_QNR, & ! locations in tracer fields for
P_QNWFA, P_QV ! use in creating concentrations

character(len=50), allocatable :: slevs(:)

Expand Down
2 changes: 0 additions & 2 deletions sorc/chgres_cube.fd/program_setup.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,8 +294,6 @@ subroutine read_varmap
character(len=20),allocatable :: var_type(:)

if (trim(input_type) == "grib2") then
!varmap_table_file = trim(base_install_dir) // "/" // trim(varmap_tables_dir) // "/" &
! // trim(phys_suite) // "phys_var_map.txt"

print*,"OPEN VARIABLE MAPPING FILE: ", trim(varmap_file)
open(14, file=trim(varmap_file), form='formatted', iostat=istat)
Expand Down
Loading

0 comments on commit 999da97

Please sign in to comment.