From 999da97595bc60b47fd69793a502e15fb6419b87 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2020 17:56:00 +0000 Subject: [PATCH] feature/chgres_cube_grib2_release: This commit references #7. Remove more code for GSD physics, which will not be supported for the next public release. --- sorc/chgres_cube.fd/atmosphere.F90 | 303 +------------------------- sorc/chgres_cube.fd/input_data.F90 | 4 - sorc/chgres_cube.fd/program_setup.f90 | 2 - sorc/chgres_cube.fd/surface.F90 | 8 - 4 files changed, 1 insertion(+), 316 deletions(-) diff --git a/sorc/chgres_cube.fd/atmosphere.F90 b/sorc/chgres_cube.fd/atmosphere.F90 index e10af0b15..e8624be9e 100644 --- a/sorc/chgres_cube.fd/atmosphere.F90 +++ b/sorc/chgres_cube.fd/atmosphere.F90 @@ -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, & @@ -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. !----------------------------------------------------------------------------------- @@ -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, & @@ -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 diff --git a/sorc/chgres_cube.fd/input_data.F90 b/sorc/chgres_cube.fd/input_data.F90 index df0f99d29..7c2f0dda4 100644 --- a/sorc/chgres_cube.fd/input_data.F90 +++ b/sorc/chgres_cube.fd/input_data.F90 @@ -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(:) diff --git a/sorc/chgres_cube.fd/program_setup.f90 b/sorc/chgres_cube.fd/program_setup.f90 index ae961418e..835f0db9f 100644 --- a/sorc/chgres_cube.fd/program_setup.f90 +++ b/sorc/chgres_cube.fd/program_setup.f90 @@ -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) diff --git a/sorc/chgres_cube.fd/surface.F90 b/sorc/chgres_cube.fd/surface.F90 index bab7b508e..01f9f1033 100644 --- a/sorc/chgres_cube.fd/surface.F90 +++ b/sorc/chgres_cube.fd/surface.F90 @@ -603,10 +603,8 @@ subroutine interp(localpet) call error_handler("IN FieldGather", rc) if (localpet == 0) then - print*,'bf search, ice target(14,86) = ',data_one_tile(14,86) call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 91, & latitude=latitude_one_tile) - print*,'af search, ice target(14,86) = ',data_one_tile(14,86) endif print*,"- CALL FieldGather FOR TARGET LANDMASK TILE: ", tile @@ -616,7 +614,6 @@ subroutine interp(localpet) if (localpet == 0) then - print*,'bf ice repl, ice, slmsk_target(14,86) = ', data_one_tile(14,86), mask_target_one_tile(14,86) do j = 1, j_target do i = 1, i_target if (data_one_tile(i,j) > 1.0_esmf_kind_r8) then @@ -626,7 +623,6 @@ subroutine interp(localpet) if (data_one_tile(i,j) >= 0.15_esmf_kind_r8) mask_target_one_tile(i,j) = 2 enddo enddo - print*,'af ice repl, ice, slmsk_target(14,86) = ', data_one_tile(14,86),mask_target_one_tile(14,86) endif print*,"- CALL FieldScatter FOR TARGET GRID SEAICE FRACTION TILE: ", tile @@ -807,9 +803,7 @@ subroutine interp(localpet) call error_handler("IN FieldGather", rc) if (localpet == 0) then - print*, "snod target bf search ice ", data_one_tile(61,80) call search(data_one_tile, mask_target_one_tile, i_target, j_target, tile, 66) - print*, "snod target af search ice ", data_one_tile(61,80) endif print*,"- CALL FieldScatter FOR TARGET GRID SNOW DEPTH TILE: ", tile @@ -1672,12 +1666,10 @@ subroutine interp(localpet) call error_handler("IN FieldGather", rc) if (localpet == 0) then - print*, "snod input bf search land = ",data_one_tile(61,80) allocate(land_target_one_tile(i_target,j_target)) land_target_one_tile = 0 where(mask_target_one_tile == 1) land_target_one_tile = 1 call search(data_one_tile, land_target_one_tile, i_target, j_target, tile, 66) - print*, "snod af search land =", data_one_tile(61,80) endif print*,"- CALL FieldScatter FOR TARGET GRID SNOW DEPTH: ", tile