From c863adc8c9445a2f3bc46190e4b6d505399f1068 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 21 Oct 2023 16:13:35 -0700 Subject: [PATCH 1/2] test(NetCDF): improve unit test description --- test/netCDF_file_test_m.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/netCDF_file_test_m.f90 b/test/netCDF_file_test_m.f90 index f203befc8..56035d396 100644 --- a/test/netCDF_file_test_m.f90 +++ b/test/netCDF_file_test_m.f90 @@ -45,13 +45,13 @@ function results() result(test_results) type(test_result_t), allocatable :: test_results(:) character(len=*), parameter :: longest_description = & - "written and then read gives input matching the output for a 1D double precision array" + "writing and then reading gives input matching the output for a 1D double precision array" associate( & descriptions => & [ character(len=len(longest_description)) :: & - "written and then read gives input matching the output for a 2D integer array", & - "written and then read gives input matching the output for a 1D double precision array" & + "writing and then reading gives input matching the output for a 2D integer array", & + "writing and then reading gives input matching the output for a 1D double precision array" & ], & outcomes => & [ write_then_read_2D_integer(), & From 2892d1189ca0504c5cd44d4946ce0c488486b3c8 Mon Sep 17 00:00:00 2001 From: Damian Rouson Date: Sat, 21 Oct 2023 23:14:34 -0700 Subject: [PATCH 2/2] example: learn two Thompson microphysics functions --- example/learn-microphysics-procedures.f90 | 231 + example/supporting-modules/mp_thompson.f90 | 3868 +++++++++++++++++ .../supporting-modules/thompson_tensors_m.f90 | 33 + setup.sh | 2 +- 4 files changed, 4133 insertions(+), 1 deletion(-) create mode 100644 example/learn-microphysics-procedures.f90 create mode 100644 example/supporting-modules/mp_thompson.f90 create mode 100644 example/supporting-modules/thompson_tensors_m.f90 diff --git a/example/learn-microphysics-procedures.f90 b/example/learn-microphysics-procedures.f90 new file mode 100644 index 000000000..094eacb99 --- /dev/null +++ b/example/learn-microphysics-procedures.f90 @@ -0,0 +1,231 @@ +! Copyright (c), The Regents of the University of California +! Terms of use are as specified in LICENSE.txt +program learn_microphysics_procedures + !! Train a neural network proxies for procedures in the Thompson microphysics model + !! in of ICAR (https://github.com/BerkeleyLab/icar). + use inference_engine_m, only : & + inference_engine_t, trainable_engine_t, mini_batch_t, tensor_t, input_output_pair_t, shuffle, sigmoid_t + use sourcery_m, only : string_t, file_t, command_line_t, bin_t, csv + use assert_m, only : assert, intrinsic_array_t + use thompson_tensors_m, only : y, T, p + use iso_fortran_env, only : int64, output_unit + implicit none + + type(string_t) network_file + type(command_line_t) command_line + integer(int64) counter_start, counter_end, clock_rate + + network_file = string_t(command_line%flag_value("--output-file")) + + if (len(network_file%string())==0) then + error stop new_line('a') // new_line('a') // & + 'Usage: ./build/run-fpm.sh run learn_microphysics_procedures -- --output-file ""' + end if + + call system_clock(counter_start, clock_rate) + + block + integer, parameter :: max_num_epochs = 10000000, num_mini_batches = 10 + integer num_pairs ! number of input/output pairs + + type(mini_batch_t), allocatable :: mini_batches(:) + type(input_output_pair_t), allocatable :: input_output_pairs(:) + type(tensor_t), allocatable :: inputs(:), desired_outputs(:) + type(trainable_engine_t) trainable_engine + type(bin_t), allocatable :: bins(:) + real, allocatable :: cost(:), random_numbers(:) + integer io_status, network_unit, plot_unit + integer, parameter :: io_success=0, diagnostics_print_interval = 1000, network_save_interval = 10000 + integer, parameter :: nodes_per_layer(*) = [2, 72, 2] + real, parameter :: cost_tolerance = 1.E-08 + + call random_init(image_distinct=.true., repeatable=.true.) + + open(newunit=network_unit, file=network_file%string(), form='formatted', status='old', iostat=io_status, action='read') + + if (io_status == io_success) then + print *,"Reading network from file " // network_file%string() + trainable_engine = trainable_engine_t(inference_engine_t(file_t(network_file))) + close(network_unit) + else + close(network_unit) + print *,"Initializing a new network" + trainable_engine = perturbed_identity_network(perturbation_magnitude=0.05, n = nodes_per_layer) + end if + call output(trainable_engine%to_inference_engine(), string_t("initial-network.json")) + + associate(num_inputs => trainable_engine%num_inputs(), num_outputs => trainable_engine%num_outputs()) + + block + integer i, j + integer, allocatable :: output_sizes(:) + inputs = [( [(tensor_t([T(i), p(j)]), j=1,size(p))], i = 1,size(T))] + num_pairs = size(inputs) + call assert(num_pairs == size(T)*size(p), "train_cloud_microphysics: inputs tensor array complete") + desired_outputs = y(inputs) + output_sizes = [(size(desired_outputs(i)%values()),i=1,size(desired_outputs))] + call assert(all([num_outputs==output_sizes]), "fit-polynomials: # outputs", intrinsic_array_t([num_outputs,output_sizes])) + end block + + input_output_pairs = input_output_pair_t(inputs, desired_outputs) + + block + integer b + bins = [(bin_t(num_items=num_pairs, num_bins=num_mini_batches, bin_number=b), b = 1, num_mini_batches)] + end block + + block + integer e, b, stop_unit, previous_epoch + real previous_clock_time + + call open_plot_file_for_appending("cost.plt", plot_unit, previous_epoch, previous_clock_time) + print *, " Epoch | Cost Function| System_Clock | Nodes per Layer" + allocate(random_numbers(2:size(input_output_pairs))) + + do e = previous_epoch + 1, previous_epoch + max_num_epochs + call random_number(random_numbers) + call shuffle(input_output_pairs, random_numbers) + mini_batches = [(mini_batch_t(input_output_pairs(bins(b)%first():bins(b)%last())), b = 1, size(bins))] + call trainable_engine%train(mini_batches, cost, adam=.true.) + call system_clock(counter_end, clock_rate) + + associate( & + cost_avg => sum(cost)/size(cost), & + cumulative_clock_time => previous_clock_time + real(counter_end - counter_start) / real(clock_rate), & + loop_ending => e == previous_epoch + max_num_epochs & + ) + write_and_exit_if_converged: & + if (cost_avg < cost_tolerance) then + call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer) + call output(trainable_engine%to_inference_engine(), network_file) + exit + end if write_and_exit_if_converged + + open(newunit=stop_unit, file="stop", form='formatted', status='old', iostat=io_status) + + write_and_exit_if_stop_file_exists: & + if (io_status==0) then + call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer) + call output(trainable_engine%to_inference_engine(), network_file) + exit + end if write_and_exit_if_stop_file_exists + + if (mod(e,diagnostics_print_interval)==0 .or. loop_ending) & + call print_diagnostics(plot_unit, e, cost_avg, cumulative_clock_time, nodes_per_layer) + if (mod(e,network_save_interval)==0 .or. loop_ending) call output(trainable_engine%to_inference_engine(), network_file) + end associate + end do + + close(plot_unit) + + report_network_performance: & + block + integer p + + associate(network_outputs => trainable_engine%infer(inputs)) + print *," Inputs (normalized) | Outputs | Desired outputs" + do p = 1, num_pairs + print "(6(G13.5,2x))", inputs(p)%values(), network_outputs(p)%values(), desired_outputs(p)%values() + end do + end associate + end block report_network_performance + + end block + + end associate + + call output(trainable_engine%to_inference_engine(), network_file) + + end block + +contains + + subroutine print_diagnostics(plot_file_unit, epoch, cost, clock, nodes) + integer, intent(in) :: plot_file_unit, epoch, nodes(:) + real, intent(in) :: cost, clock + + write(unit=output_unit, fmt='(3(g13.5,2x))', advance='no') epoch, cost, clock + write(unit=output_unit, fmt=csv) nodes + write(unit=plot_file_unit, fmt='(3(g13.5,2x))', advance='no') epoch, cost, clock + write(unit=plot_file_unit, fmt=csv) nodes + end subroutine + + subroutine output(inference_engine, file_name) + type(inference_engine_t), intent(in) :: inference_engine + type(string_t), intent(in) :: file_name + type(file_t) json_file + json_file = inference_engine%to_json() + call json_file%write_lines(file_name) + end subroutine + + pure function e(j,n) result(unit_vector) + integer, intent(in) :: j, n + integer k + real, allocatable :: unit_vector(:) + unit_vector = real([(merge(1,0,j==k),k=1,n)]) + end function + + function perturbed_identity_network(perturbation_magnitude, n) result(trainable_engine) + type(trainable_engine_t) trainable_engine + real, intent(in) :: perturbation_magnitude + integer, intent(in) :: n(:) + integer j, k, l + real, allocatable :: identity(:,:,:), w_harvest(:,:,:), b_harvest(:,:) + + associate(n_max => maxval(n), layers => size(n)) + identity = reshape( [( [(e(k,n_max), k=1,n_max)], l = 1, layers-1 )], [n_max, n_max, layers-1]) + + allocate(w_harvest, mold = identity) + allocate(b_harvest(size(identity,1), size(identity,3))) + + call random_number(w_harvest) + call random_number(b_harvest) + + associate(w => identity + perturbation_magnitude*(w_harvest-0.5)/0.5, b => perturbation_magnitude*(b_harvest-0.5)/0.5) + + trainable_engine = trainable_engine_t( & + nodes = n, weights = w, biases = b, differentiable_activation_strategy = sigmoid_t(), & + metadata = & + [string_t("Thompson microphysics procedures"), string_t("Damian Rouson"), string_t("2023-09-23"), string_t("sigmoid"), & + string_t("false")] & + ) + end associate + end associate + end function + + subroutine open_plot_file_for_appending(plot_file_name, plot_unit, previous_epoch, previous_clock) + character(len=*), intent(in) :: plot_file_name + integer, intent(out) :: plot_unit, previous_epoch + real, intent(out) :: previous_clock + + type(file_t) plot_file + type(string_t), allocatable :: lines(:) + character(len=:), allocatable :: last_line + integer io_status + integer, parameter :: io_success = 0 + logical preexisting_plot_file + real cost + + inquire(file=plot_file_name, exist=preexisting_plot_file) + open(newunit=plot_unit,file="cost.plt",status="unknown",position="append") + + associate(header => " Epoch | Cost Function| System_Clock | Nodes per Layer") + if (.not. preexisting_plot_file) then + write(plot_unit,*) header + previous_epoch = 0 + previous_clock = 0 + else + plot_file = file_t(string_t(plot_file_name)) + lines = plot_file%lines() + last_line = lines(size(lines))%string() + read(last_line,*, iostat=io_status) previous_epoch, cost, previous_clock + if ((io_status /= io_success .and. last_line == header) .or. len(trim(last_line))==0) then + previous_epoch = 0 + previous_clock = 0 + end if + end if + end associate + + end subroutine + +end program diff --git a/example/supporting-modules/mp_thompson.f90 b/example/supporting-modules/mp_thompson.f90 new file mode 100644 index 000000000..4de3637ed --- /dev/null +++ b/example/supporting-modules/mp_thompson.f90 @@ -0,0 +1,3868 @@ +! Copyright (c), The Regents of the University of California +! Terms of use are as specified in LICENSE.txt + +! MIT License +! +! Copyright (c) 2017 National Center for Atmospheric Research +! +! Permission is hereby granted, free of charge, to any person obtaining a copy +! of this software and associated documentation files (the "Software"), to deal +! in the Software without restriction, including without limitation the rights +! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +! copies of the Software, and to permit persons to whom the Software is +! furnished to do so, subject to the following conditions: +! +! The above copyright notice and this permission notice shall be included in all +! copies or substantial portions of the Software. +! +! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +! SOFTWARE. + +!>---------------------------------------------------------------------- +!! This subroutine computes the moisture tendencies of water vapor, +!! cloud droplets, rain, cloud ice (pristine), snow, and graupel. +!! +!! Prior to WRFv2.2 this code was based on Reisner et al (1998), but +!! few of those pieces remain. A complete description is now found in +!! Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall, 2008: +!! Explicit Forecasts of winter precipitation using an improved bulk +!! microphysics scheme. Part II: Implementation of a new snow +!! parameterization. Mon. Wea. Rev., 136, 5095-5115. +!! Prior to WRFv3.1, this code was single-moment rain prediction as +!! described in the reference above, but in v3.1 and higher, the +!! scheme is two-moment rain (predicted rain number concentration). +!! +!! Most importantly, users may wish to modify the prescribed number of +!! cloud droplets (Nt_c; see guidelines mentioned below). Otherwise, +!! users may alter the rain and graupel size distribution parameters +!! to use exponential (Marshal-Palmer) or generalized gamma shape. +!! The snow field assumes a combination of two gamma functions (from +!! Field et al. 2005) and would require significant modifications +!! throughout the entire code to alter its shape as well as accretion +!! rates. Users may also alter the constants used for density of rain, +!! graupel, ice, and snow, but the latter is not constant when using +!! Paul Field's snow distribution and moments methods. Other values +!! users can modify include the constants for mass and/or velocity +!! power law relations and assumed capacitances used in deposition/ +!! sublimation/evaporation/melting. +!! Remaining values should probably be left alone. +!! +!! Modifications for ICAR include OPENMP paralellization and the ability +!! to input many important parameters so they are no longer hard coded +!! +!! @author Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805 +!! Last modified: 27 Jul 2012 +!! +!!-------------------------------------------------------------------- +!wrft:model_layer:physics +!+---+-----------------------------------------------------------------+ +! + MODULE module_mp_thompson +!!! use options_types, only: mp_options_type +!!! commented for compilation outside ICAR + +! USE module_wrf_error +! USE module_mp_radar +! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm +! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep + + IMPLICIT NONE + + LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. + INTEGER, PARAMETER, PRIVATE:: IFDRY = 0 + REAL, PARAMETER, PRIVATE:: T_0 = 273.15 +! REAL, PARAMETER, PRIVATE:: PI = 3.1415926536 !trude added. Note, pi is defined in data_structures, and conflict with definition here. Need to determine what to to about it. + REAL, PARAMETER, PRIVATE:: PI2 = 3.1415926536 !trude added. Note, pi is defined in data_structures, and conflict with definition here. Need to determine what to to about it. + +!..Densities of rain, snow, graupel, and cloud ice. + REAL, PARAMETER, PRIVATE:: rho_w = 1000.0 + REAL, PARAMETER, PRIVATE:: rho_s = 100.0 + ! REAL, PARAMETER, PRIVATE:: rho_g = 500.0 ! trude commented out for changing from parameter to input variable + REAL, PARAMETER, PRIVATE:: rho_i = 890.0 + +!..Prescribed number of cloud droplets. Set according to known data or +!.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and +!.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter, +!.. mu_c, calculated based on Nt_c is important in autoconversion +!.. scheme. +! ++ trude comment out +! REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6 +! -- trude comment out +!..Generalized gamma distributions for rain, graupel and cloud ice. +!.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential. +! REAL, PARAMETER, PRIVATE:: mu_r = 0.0 + REAL, PARAMETER, PRIVATE:: mu_g = 0.0 + REAL, PARAMETER, PRIVATE:: mu_i = 0.0 + REAL, PRIVATE:: mu_c + +!..Sum of two gamma distrib for snow (Field et al. 2005). +!.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3) +!.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)] +!.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively +!.. calculated as function of ice water content and temperature. + REAL, PARAMETER, PRIVATE:: mu_s = 0.6357 + REAL, PARAMETER, PRIVATE:: Kap0 = 490.6 + REAL, PARAMETER, PRIVATE:: Kap1 = 17.46 + REAL, PARAMETER, PRIVATE:: Lam0 = 20.78 + REAL, PARAMETER, PRIVATE:: Lam1 = 3.29 + +!..Y-intercept parameter for graupel is not constant and depends on +!.. mixing ratio. Also, when mu_g is non-zero, these become equiv +!.. y-intercept for an exponential distrib and proper values are +!.. computed based on same mixing ratio and total number concentration. + REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4 + REAL, PARAMETER, PRIVATE:: gonv_max = 3.E6 + +!..Mass power law relations: mass = am*D**bm +!.. Snow from Field et al. (2005), others assume spherical form. + REAL, PARAMETER, PRIVATE:: am_r = PI2*rho_w/6.0 + REAL, PARAMETER, PRIVATE:: bm_r = 3.0 +! REAL, PARAMETER, PRIVATE:: am_s = 0.069 ! trude commented out for changing from parameter to input variable + REAL, PARAMETER, PRIVATE:: bm_s = 2.0 +! REAL, PARAMETER, PRIVATE:: am_g = PI2*rho_g/6.0 ! trude commented out. am_g need to be calculated later since rho_g is an imput variable + REAL, PARAMETER, PRIVATE:: bm_g = 3.0 + REAL, PARAMETER, PRIVATE:: am_i = PI2*rho_i/6.0 + REAL, PARAMETER, PRIVATE:: bm_i = 3.0 + +!..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D) +!.. Rain from Ferrier (1994), ice, snow, and graupel from +!.. Thompson et al (2008). Coefficient fv is zero for graupel/ice. + REAL, PARAMETER, PRIVATE:: av_r = 4854.0 + REAL, PARAMETER, PRIVATE:: bv_r = 1.0 + REAL, PARAMETER, PRIVATE:: fv_r = 195.0 +! REAL, PARAMETER, PRIVATE:: av_s = 40.0 ! trude commented out. Will be used as input vaiable. +! REAL, PARAMETER, PRIVATE:: bv_s = 0.55 ! trude commented out. Will be used as input vaiable. +! REAL, PARAMETER, PRIVATE:: fv_s = 100.0 ! trude commented out. Will be used as input vaiable. +! REAL, PARAMETER, PRIVATE:: av_g = 442.0 ! trude commented out. Will be used as input vaiable. +! REAL, PARAMETER, PRIVATE:: bv_g = 0.89 ! trude commented out. Will be used as input vaiable. +! REAL, PARAMETER, PRIVATE:: av_i = 1847.5 + REAL, PARAMETER, PRIVATE:: bv_i = 1.0 + +!..Capacitance of sphere and plates/aggregates: D**3, D**2 + REAL, PARAMETER, PRIVATE:: C_cube = 0.5 +! REAL, PARAMETER, PRIVATE:: C_sqrd = 0.3 + +!..Collection efficiencies. Rain/snow/graupel collection of cloud +!.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and +!.. get computed elsewhere because they are dependent on stokes +!.. number. +! REAL, PARAMETER, PRIVATE:: Ef_si = 0.05 +! REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95 +! REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75 +! REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95 + +!..Minimum microphys values +!.. R1 value, 1.E-12, cannot be set lower because of numerical +!.. problems with Paul Field's moments and should not be set larger +!.. because of truncation problems in snow/ice growth. + REAL, PARAMETER, PRIVATE:: R1 = 1.E-12 + REAL, PARAMETER, PRIVATE:: R2 = 1.E-6 + REAL, PARAMETER, PRIVATE:: eps = 1.E-15 + +!..Constants in Cooper curve relation for cloud ice number. +! REAL, PARAMETER, PRIVATE:: TNO = 5.0 ! trude comment this out for use as input variable instead + REAL, PARAMETER, PRIVATE:: ATO = 0.304 + +!..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment. + REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0) + +!..Schmidt number + REAL, PARAMETER, PRIVATE:: Sc = 0.632 + REAL, PRIVATE:: Sc3 + +!..Homogeneous freezing temperature + REAL, PARAMETER, PRIVATE:: HGFR = 235.16 + +!..Water vapor and air gas constants at constant pressure + REAL, PARAMETER, PRIVATE:: Rv = 461.5 + REAL, PARAMETER, PRIVATE:: oRv = 1./Rv +! REAL, PARAMETER, PRIVATE:: R = 287.04 !trude added. Note, R = 287.058 is defined in data_structures, and conflict with definition here. Need to determine what to to about it. + REAL, PARAMETER, PRIVATE:: RR2 = 287.04 !trude added. Note, R = 287.058 is defined in data_structures, and conflict with definition here. Need to determine what to to about it. +! REAL, PARAMETER, PRIVATE:: Cp = 1004.0 !trude added. Note, Cp = 1012.0 is defined in data_structures, and conflict with definition here. Need to determine what to to about it. + REAL, PARAMETER, PRIVATE:: Cp2 = 1004.0 !trude added. Note, Cp = 1012.0 is defined in data_structures, and conflict with definition here. Need to determine what to to about it. + +!..Enthalpy of sublimation, vaporization, and fusion at 0C. + REAL, PARAMETER, PRIVATE:: lsub = 2.834E6 + REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6 + REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0 + REAL, PARAMETER, PRIVATE:: olfus = 1./lfus + +!..Ice initiates with this mass (kg), corresponding diameter calc. +!..Min diameters and mass of cloud, rain, snow, and graupel (m, kg). + REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12 + REAL, PARAMETER, PRIVATE:: D0c = 1.E-6 + REAL, PARAMETER, PRIVATE:: D0r = 50.E-6 + REAL, PARAMETER, PRIVATE:: D0s = 200.E-6 + REAL, PARAMETER, PRIVATE:: D0g = 250.E-6 + REAL, PRIVATE:: D0i, xm0s, xm0g + +!..Lookup table dimensions + INTEGER, PARAMETER, PRIVATE:: nbins = 100 + INTEGER, PARAMETER, PRIVATE:: nbc = nbins + INTEGER, PARAMETER, PRIVATE:: nbi = nbins + INTEGER, PARAMETER, PRIVATE:: nbr = nbins + INTEGER, PARAMETER, PRIVATE:: nbs = nbins + INTEGER, PARAMETER, PRIVATE:: nbg = nbins + INTEGER, PARAMETER, PRIVATE:: ntb_c = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_i = 64 + INTEGER, PARAMETER, PRIVATE:: ntb_r = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_s = 28 + INTEGER, PARAMETER, PRIVATE:: ntb_g = 28 + INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28 + INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37 + INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55 + INTEGER, PARAMETER, PRIVATE:: ntb_t = 9 + INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3 + + DOUBLE PRECISION, DIMENSION(nbins+1):: xDx + DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc + DOUBLE PRECISION, DIMENSION(nbi):: Di, dti + DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr + DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts + DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg + +!..Lookup tables for cloud water content (kg/m**3). + REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: & + r_c = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for cloud ice content (kg/m**3). + REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: & + r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, & + 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, & + 1.e-9,2.e-9,3.e-9,4.e-9,5.e-9,6.e-9,7.e-9,8.e-9,9.e-9, & + 1.e-8,2.e-8,3.e-8,4.e-8,5.e-8,6.e-8,7.e-8,8.e-8,9.e-8, & + 1.e-7,2.e-7,3.e-7,4.e-7,5.e-7,6.e-7,7.e-7,8.e-7,9.e-7, & + 1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3/) + +!..Lookup tables for rain content (kg/m**3). + REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: & + r_r = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & + 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for graupel content (kg/m**3). + REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: & + r_g = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for snow content (kg/m**3). + REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: & + r_s = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & + 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & + 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & + 1.e-2/) + +!..Lookup tables for rain y-intercept parameter (/m**4). + REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: & + N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & + 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, & + 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, & + 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, & + 1.e10/) + +!..Lookup tables for graupel y-intercept parameter (/m**4). + REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: & + N0g_exp = (/1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & + 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & + 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & + 1.e7/) + +!..Lookup tables for ice number concentration (/m**3). + REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: & + Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, & + 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, & + 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & + 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & + 1.e6/) + +!..For snow moments conversions (from Field et al. 2005) + REAL, DIMENSION(10), PARAMETER, PRIVATE:: & + sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, & + 0.31255, 0.000204, 0.003199, 0.0, -0.015952/) + REAL, DIMENSION(10), PARAMETER, PRIVATE:: & + sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, & + 0.060366, 0.000079, 0.000594, 0.0, -0.003577/) + +!..Temperatures (5 C interval 0 to -40) used in lookup tables. + REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: & + Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./) + +!..Lookup tables for various accretion/collection terms. +!.. ntb_x refers to the number of elements for rain, snow, graupel, +!.. and temperature array indices. Variables beginning with t-p/c/m/n +!.. represent lookup tables. Save compile-time memory by making +!.. allocatable (2009Jun12, J. Michalakes). + INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8 + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, & + tnr_racg, tnr_gacr + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & + tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, & + tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, & + tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2 + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: & + tpi_qcfz, tni_qcfz + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: & + tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: & + tps_iaus, tni_iaus, tpi_ide + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw + REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev + +!..Variables holding a bunch of exponents and gamma values (cloud water, +!.. cloud ice, rain, snow, then graupel). + REAL, DIMENSION(3), PRIVATE:: cce, ccg + REAL, PRIVATE:: ocg1, ocg2 + REAL, DIMENSION(7), PRIVATE:: cie, cig + REAL, PRIVATE:: oig1, oig2, obmi + REAL, DIMENSION(13), PRIVATE:: cre, crg + REAL, PRIVATE:: ore1, org1, org2, org3, obmr + REAL, DIMENSION(18), PRIVATE:: cse, csg + REAL, PRIVATE:: oams, obms, ocms + REAL, DIMENSION(12), PRIVATE:: cge, cgg + REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg + +!..Declaration of precomputed constants in various rate eqns. + REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi + REAL:: t1_qr_ev, t2_qr_ev + REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd + REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me + + CHARACTER*256:: mp_debug + +! TRUDE + REAL, PRIVATE :: Nt_c, TNO, am_s,rho_g,av_s,bv_s,fv_s,av_g,bv_g,av_i,Ef_si,Ef_rs,Ef_rg,Ef_ri + REAL, PRIVATE :: C_cubes,C_sqrd, mu_r, t_adjust + LOGICAL, PRIVATE :: Ef_rw_l, Ef_sw_l + REAL, PRIVATE :: am_g + + +!+---+ +!+---+-----------------------------------------------------------------+ +!..END DECLARATIONS +!+---+-----------------------------------------------------------------+ +!+---+ +!ctrlL + + CONTAINS + +!!! The subroutine below is commented to support compiling this module outside ICAR + +!!! SUBROUTINE thompson_init(mp_options) +!!! +!!! IMPLICIT NONE +!!! +!!! ! ++ trude +!!! type(mp_options_type), intent(in) :: mp_options +!!! ! -- trude +!!! INTEGER:: i, j, k, m, n +!!! LOGICAL:: micro_init +!!! !..Allocate space for lookup tables (J. Michalakes 2009Jun08). +!!! micro_init = .FALSE. +!!! if (.NOT. ALLOCATED(tcg_racg) ) then +!!! ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! micro_init = .TRUE. +!!! endif +!!! +!!! if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) +!!! +!!! if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) +!!! +!!! if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,45)) +!!! if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,45)) +!!! +!!! if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45)) +!!! if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45)) +!!! if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45)) +!!! if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45)) +!!! +!!! if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1)) +!!! if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1)) +!!! if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1)) +!!! +!!! if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc)) +!!! if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc)) +!!! +!!! if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r)) +!!! +!!! if (micro_init) then +!!! +!!! !++ trude +!!! Nt_c = mp_options%Nt_c +!!! TNO = mp_options%TNO +!!! am_s = mp_options%am_s +!!! rho_g = mp_options%rho_g +!!! av_s = mp_options%av_s +!!! bv_s = mp_options%bv_s +!!! fv_s = mp_options%fv_s +!!! av_g = mp_options%av_g +!!! bv_g = mp_options%bv_g +!!! av_i = mp_options%av_i +!!! Ef_rs = mp_options%Ef_rs +!!! Ef_rg = mp_options%Ef_rg +!!! Ef_si = mp_options%Ef_si +!!! Ef_ri = mp_options%Ef_ri +!!! C_cubes = mp_options%C_cubes +!!! C_sqrd = mp_options%C_sqrd +!!! mu_r = mp_options%mu_r +!!! t_adjust = mp_options%t_adjust +!!! Ef_rw_l = mp_options%Ef_rw_l +!!! Ef_sw_l = mp_options%Ef_sw_l +!!! am_g = PI2*rho_g/6.0 ! trude, nb, this should can be defined in the modular section +!!! !--trude +!!! !..From Martin et al. (1994), assign gamma shape parameter mu for cloud +!!! !.. drops according to general dispersion characteristics (disp=~0.25 +!!! !.. for Maritime and 0.45 for Continental). +!!! !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime +!!! !.. to 2 for really dirty air. +!!! mu_c = MIN(15., (1000.E6/Nt_c + 2.)) +!!! +!!! !..Schmidt number to one-third used numerous times. +!!! Sc3 = Sc**(1./3.) +!!! +!!! !..Compute min ice diam from mass, min snow/graupel mass from diam. +!!! D0i = (xm0i/am_i)**(1./bm_i) +!!! xm0s = am_s * D0s**bm_s +!!! xm0g = am_g * D0g**bm_g +!!! +!!! !..These constants various exponents and gamma() assoc with cloud, +!!! !.. rain, snow, and graupel. +!!! cce(1) = mu_c + 1. +!!! cce(2) = bm_r + mu_c + 1. +!!! cce(3) = bm_r + mu_c + 4. +!!! ccg(1) = WGAMMA(cce(1)) +!!! ccg(2) = WGAMMA(cce(2)) +!!! ccg(3) = WGAMMA(cce(3)) +!!! ocg1 = 1./ccg(1) +!!! ocg2 = 1./ccg(2) +!!! +!!! cie(1) = mu_i + 1. +!!! cie(2) = bm_i + mu_i + 1. +!!! cie(3) = bm_i + mu_i + bv_i + 1. +!!! cie(4) = mu_i + bv_i + 1. +!!! cie(5) = mu_i + 2. +!!! cie(6) = bm_i*0.5 + mu_i + bv_i + 1. +!!! cie(7) = bm_i*0.5 + mu_i + 1. +!!! cig(1) = WGAMMA(cie(1)) +!!! cig(2) = WGAMMA(cie(2)) +!!! cig(3) = WGAMMA(cie(3)) +!!! cig(4) = WGAMMA(cie(4)) +!!! cig(5) = WGAMMA(cie(5)) +!!! cig(6) = WGAMMA(cie(6)) +!!! cig(7) = WGAMMA(cie(7)) +!!! oig1 = 1./cig(1) +!!! oig2 = 1./cig(2) +!!! obmi = 1./bm_i +!!! +!!! cre(1) = bm_r + 1. +!!! cre(2) = mu_r + 1. +!!! cre(3) = bm_r + mu_r + 1. +!!! cre(4) = bm_r*2. + mu_r + 1. +!!! cre(5) = mu_r + bv_r + 1. +!!! cre(6) = bm_r + mu_r + bv_r + 1. +!!! cre(7) = bm_r*0.5 + mu_r + bv_r + 1. +!!! cre(8) = bm_r + mu_r + bv_r + 3. +!!! cre(9) = mu_r + bv_r + 3. +!!! cre(10) = mu_r + 2. +!!! cre(11) = 0.5*(bv_r + 5. + 2.*mu_r) +!!! cre(12) = bm_r*0.5 + mu_r + 1. +!!! cre(13) = bm_r*2. + mu_r + bv_r + 1. +!!! do n = 1, 13 +!!! crg(n) = WGAMMA(cre(n)) +!!! enddo +!!! obmr = 1./bm_r +!!! ore1 = 1./cre(1) +!!! org1 = 1./crg(1) +!!! org2 = 1./crg(2) +!!! org3 = 1./crg(3) +!!! +!!! cse(1) = bm_s + 1. +!!! cse(2) = bm_s + 2. +!!! cse(3) = bm_s*2. +!!! cse(4) = bm_s + bv_s + 1. +!!! cse(5) = bm_s*2. + bv_s + 1. +!!! cse(6) = bm_s*2. + 1. +!!! cse(7) = bm_s + mu_s + 1. +!!! cse(8) = bm_s + mu_s + 2. +!!! cse(9) = bm_s + mu_s + 3. +!!! cse(10) = bm_s + mu_s + bv_s + 1. +!!! cse(11) = bm_s*2. + mu_s + bv_s + 1. +!!! cse(12) = bm_s*2. + mu_s + 1. +!!! cse(13) = bv_s + 2. +!!! cse(14) = bm_s + bv_s +!!! cse(15) = mu_s + 1. +!!! cse(16) = 1.0 + (1.0 + bv_s)/2. +!!! cse(17) = cse(16) + mu_s + 1. +!!! cse(18) = bv_s + mu_s + 3. +!!! do n = 1, 18 +!!! csg(n) = WGAMMA(cse(n)) +!!! enddo +!!! oams = 1./am_s +!!! obms = 1./bm_s +!!! ocms = oams**obms +!!! +!!! cge(1) = bm_g + 1. +!!! cge(2) = mu_g + 1. +!!! cge(3) = bm_g + mu_g + 1. +!!! cge(4) = bm_g*2. + mu_g + 1. +!!! cge(5) = bm_g*2. + mu_g + bv_g + 1. +!!! cge(6) = bm_g + mu_g + bv_g + 1. +!!! cge(7) = bm_g + mu_g + bv_g + 2. +!!! cge(8) = bm_g + mu_g + bv_g + 3. +!!! cge(9) = mu_g + bv_g + 3. +!!! cge(10) = mu_g + 2. +!!! cge(11) = 0.5*(bv_g + 5. + 2.*mu_g) +!!! cge(12) = 0.5*(bv_g + 5.) + mu_g +!!! do n = 1, 12 +!!! cgg(n) = WGAMMA(cge(n)) +!!! enddo +!!! oamg = 1./am_g +!!! obmg = 1./bm_g +!!! ocmg = oamg**obmg +!!! oge1 = 1./cge(1) +!!! ogg1 = 1./cgg(1) +!!! ogg2 = 1./cgg(2) +!!! ogg3 = 1./cgg(3) +!!! +!!! !+---+-----------------------------------------------------------------+ +!!! !..Simplify various rate eqns the best we can now. +!!! !+---+-----------------------------------------------------------------+ +!!! +!!! !..Rain collecting cloud water and cloud ice +!!! t1_qr_qc = PI2*.25*av_r * crg(9) +!!! t1_qr_qi = PI2*.25*av_r * crg(9) +!!! t2_qr_qi = PI2*.25*am_r*av_r * crg(8) +!!! +!!! !..Graupel collecting cloud water +!!! t1_qg_qc = PI2*.25*av_g * cgg(9) +!!! +!!! !..Snow collecting cloud water +!!! t1_qs_qc = PI2*.25*av_s +!!! +!!! !..Snow collecting cloud ice +!!! t1_qs_qi = PI2*.25*av_s +!!! +!!! !..Evaporation of rain; ignore depositional growth of rain. +!!! t1_qr_ev = 0.78 * crg(10) +!!! t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11) +!!! !..Sublimation/depositional growth of snow +!!! t1_qs_sd = 0.86 +!!! t2_qs_sd = 0.28*Sc3*SQRT(av_s) +!!! +!!! !..Melting of snow +!!! t1_qs_me = PI2*4.*C_sqrd*olfus * 0.86 +!!! t2_qs_me = PI2*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s) +!!! +!!! !..Sublimation/depositional growth of graupel +!!! t1_qg_sd = 0.86 * cgg(10) +!!! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11) +!!! +!!! !..Melting of graupel +!!! t1_qg_me = PI2*4.*C_cube*olfus * 0.86 * cgg(10) +!!! t2_qg_me = PI2*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11) +!!! +!!! !..Constants for helping find lookup table indexes. +!!! nic2 = NINT(ALOG10(r_c(1))) +!!! nii2 = NINT(ALOG10(r_i(1))) +!!! nii3 = NINT(ALOG10(Nt_i(1))) +!!! nir2 = NINT(ALOG10(r_r(1))) +!!! nir3 = NINT(ALOG10(N0r_exp(1))) +!!! nis2 = NINT(ALOG10(r_s(1))) +!!! nig2 = NINT(ALOG10(r_g(1))) +!!! nig3 = NINT(ALOG10(N0g_exp(1))) +!!! +!!! !..Create bins of cloud water (from min diameter up to 100 microns). +!!! Dc(1) = D0c*1.0d0 +!!! dtc(1) = D0c*1.0d0 +!!! do n = 2, nbc +!!! Dc(n) = Dc(n-1) + 1.0D-6 +!!! dtc(n) = (Dc(n) - Dc(n-1)) +!!! enddo +!!! +!!! !..Create bins of cloud ice (from min diameter up to 5x min snow size). +!!! xDx(1) = D0i*1.0d0 +!!! xDx(nbi+1) = 5.0d0*D0s +!!! do n = 2, nbi +!!! xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) & +!!! *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1))) +!!! enddo +!!! do n = 1, nbi +!!! Di(n) = DSQRT(xDx(n)*xDx(n+1)) +!!! dti(n) = xDx(n+1) - xDx(n) +!!! enddo +!!! +!!! !..Create bins of rain (from min diameter up to 5 mm). +!!! xDx(1) = D0r*1.0d0 +!!! xDx(nbr+1) = 0.005d0 +!!! do n = 2, nbr +!!! xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) & +!!! *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1))) +!!! enddo +!!! do n = 1, nbr +!!! Dr(n) = DSQRT(xDx(n)*xDx(n+1)) +!!! dtr(n) = xDx(n+1) - xDx(n) +!!! enddo +!!! +!!! !..Create bins of snow (from min diameter up to 2 cm). +!!! xDx(1) = D0s*1.0d0 +!!! xDx(nbs+1) = 0.02d0 +!!! do n = 2, nbs +!!! xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) & +!!! *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1))) +!!! enddo +!!! do n = 1, nbs +!!! Ds(n) = DSQRT(xDx(n)*xDx(n+1)) +!!! dts(n) = xDx(n+1) - xDx(n) +!!! enddo +!!! +!!! !..Create bins of graupel (from min diameter up to 5 cm). +!!! xDx(1) = D0g*1.0d0 +!!! xDx(nbg+1) = 0.05d0 +!!! do n = 2, nbg +!!! xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) & +!!! *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1))) +!!! enddo +!!! do n = 1, nbg +!!! Dg(n) = DSQRT(xDx(n)*xDx(n+1)) +!!! dtg(n) = xDx(n+1) - xDx(n) +!!! enddo +!!! +!!! !+---+-----------------------------------------------------------------+ +!!! !..Create lookup tables for most costly calculations. +!!! !+---+-----------------------------------------------------------------+ +!!! +!!! do m = 1, ntb_r +!!! do k = 1, ntb_r1 +!!! do j = 1, ntb_g +!!! do i = 1, ntb_g1 +!!! tcg_racg(i,j,k,m) = 0.0d0 +!!! tmr_racg(i,j,k,m) = 0.0d0 +!!! tcr_gacr(i,j,k,m) = 0.0d0 +!!! tmg_gacr(i,j,k,m) = 0.0d0 +!!! tnr_racg(i,j,k,m) = 0.0d0 +!!! tnr_gacr(i,j,k,m) = 0.0d0 +!!! enddo +!!! enddo +!!! enddo +!!! enddo +!!! +!!! do m = 1, ntb_r +!!! do k = 1, ntb_r1 +!!! do j = 1, ntb_t +!!! do i = 1, ntb_s +!!! tcs_racs1(i,j,k,m) = 0.0d0 +!!! tmr_racs1(i,j,k,m) = 0.0d0 +!!! tcs_racs2(i,j,k,m) = 0.0d0 +!!! tmr_racs2(i,j,k,m) = 0.0d0 +!!! tcr_sacr1(i,j,k,m) = 0.0d0 +!!! tms_sacr1(i,j,k,m) = 0.0d0 +!!! tcr_sacr2(i,j,k,m) = 0.0d0 +!!! tms_sacr2(i,j,k,m) = 0.0d0 +!!! tnr_racs1(i,j,k,m) = 0.0d0 +!!! tnr_racs2(i,j,k,m) = 0.0d0 +!!! tnr_sacr1(i,j,k,m) = 0.0d0 +!!! tnr_sacr2(i,j,k,m) = 0.0d0 +!!! enddo +!!! enddo +!!! enddo +!!! enddo +!!! +!!! do k = 1, 45 +!!! do j = 1, ntb_r1 +!!! do i = 1, ntb_r +!!! tpi_qrfz(i,j,k) = 0.0d0 +!!! tni_qrfz(i,j,k) = 0.0d0 +!!! tpg_qrfz(i,j,k) = 0.0d0 +!!! tnr_qrfz(i,j,k) = 0.0d0 +!!! enddo +!!! enddo +!!! do i = 1, ntb_c +!!! tpi_qcfz(i,k) = 0.0d0 +!!! tni_qcfz(i,k) = 0.0d0 +!!! enddo +!!! enddo +!!! +!!! do j = 1, ntb_i1 +!!! do i = 1, ntb_i +!!! tps_iaus(i,j) = 0.0d0 +!!! tni_iaus(i,j) = 0.0d0 +!!! tpi_ide(i,j) = 0.0d0 +!!! enddo +!!! enddo +!!! +!!! do j = 1, nbc +!!! do i = 1, nbr +!!! t_Efrw(i,j) = 0.0 +!!! enddo +!!! do i = 1, nbs +!!! t_Efsw(i,j) = 0.0 +!!! enddo +!!! enddo +!!! +!!! do k = 1, ntb_r +!!! do j = 1, ntb_r1 +!!! do i = 1, nbr +!!! tnr_rev(i,j,k) = 0.0d0 +!!! enddo +!!! enddo +!!! enddo +!!! +!!! ! CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ') +!!! ! WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & +!!! ! ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g +!!! ! CALL wrf_debug(150, wrf_err_message) +!!! +!!! !..Collision efficiency between rain/snow and cloud water. +!!! ! CALL wrf_debug(200, ' creating qc collision eff tables') +!!! call table_Efrw +!!! call table_Efsw +!!! +!!! !..Drop evaporation. +!!! ! CALL wrf_debug(200, ' creating rain evap table') +!!! ! call table_dropEvap +!!! +!!! !..Initialize various constants for computing radar reflectivity. +!!! ! xam_r = am_r +!!! ! xbm_r = bm_r +!!! ! xmu_r = mu_r +!!! ! xam_s = am_s +!!! ! xbm_s = bm_s +!!! ! xmu_s = mu_s +!!! ! xam_g = am_g +!!! ! xbm_g = bm_g +!!! ! xmu_g = mu_g +!!! ! call radar_init +!!! +!!! if (.not. iiwarm) then +!!! +!!! !..Rain collecting graupel & graupel collecting rain. +!!! ! CALL wrf_debug(200, ' creating rain collecting graupel table') +!!! call qr_acr_qg +!!! +!!! !..Rain collecting snow & snow collecting rain. +!!! ! CALL wrf_debug(200, ' creating rain collecting snow table') +!!! call qr_acr_qs +!!! +!!! !..Cloud water and rain freezing (Bigg, 1953). +!!! ! CALL wrf_debug(200, ' creating freezing of water drops table') +!!! call freezeH2O +!!! +!!! !..Conversion of some ice mass into snow category. +!!! ! CALL wrf_debug(200, ' creating ice converting to snow table') +!!! call qi_aut_qs +!!! +!!! endif +!!! +!!! ! CALL wrf_debug(150, ' ... DONE microphysical lookup tables') +!!! +!!! endif +!!! +!!! END SUBROUTINE thompson_init +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..This is a wrapper routine designed to transfer values from 3D to 1D. +!+---+-----------------------------------------------------------------+ + SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, & + th, pii, p, dz, dt_in, itimestep, & + RAINNC, RAINNCV, & + SNOWNC, SNOWNCV, & + GRAUPELNC, GRAUPELNCV, SR, & + ids,ide, jds,jde, kds,kde, & ! domain dims + ims,ime, jms,jme, kms,kme, & ! memory dims + its,ite, jts,jte, kts,kte) ! tile dims + + implicit none + +!..Subroutine arguments + INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & + qv, qc, qr, qi, qs, qg, ni, nr, th + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & + pii, p, dz + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & + RAINNC, RAINNCV, SR + REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: & + SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV +! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & +! refl_10cm + REAL, INTENT(IN):: dt_in +!..Local variables + INTEGER, INTENT(IN):: itimestep + REAL, DIMENSION(kts:kte):: & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, t1d, p1d, dz1d, dBZ +! REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic + REAL:: dt, pptrain, pptsnow, pptgraul, pptice + REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max + INTEGER:: i, j, k + INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr + INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr + INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr + INTEGER:: i_start, j_start, i_end, j_end + + +!$OMP PARALLEL DEFAULT(PRIVATE) FIRSTPRIVATE(ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,& +!$OMP kms,kme,its,ite,jts,jte,kts,kte,itimestep) & +!$OMP SHARED(RAINNCV,RAINNC,SNOWNCV,SNOWNC,GRAUPELNCV,GRAUPELNC,SR,th,pii,p,dz,qv,qc,& +!$OMP qi,qr,qs,qg,ni,nr,dt_in,Nt_c,TNO,rho_g,av_s,bv_s,fv_s,av_g,bv_g,EF_si,Ef_ri) +! old w/pcp_xx vars !!$OMP qi,qr,qs,qg,ni,nr,pcp_ra,pcp_sn,pcp_gr,pcp_ic,dt_in) + + i_start = its + j_start = jts + i_end = MIN(ite, ide-1) + j_end = MIN(jte, jde-1) + + +!..For idealized testing by developer. +! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. & +! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then +! i_start = its + 2 +! i_end = ite +! j_start = jts +! j_end = jte +! endif + + dt = dt_in + + qc_max = 0. + qr_max = 0. + qs_max = 0. + qi_max = 0. + qg_max = 0 + ni_max = 0. + nr_max = 0. + imax_qc = 0 + imax_qr = 0 + imax_qi = 0 + imax_qs = 0 + imax_qg = 0 + imax_ni = 0 + imax_nr = 0 + jmax_qc = 0 + jmax_qr = 0 + jmax_qi = 0 + jmax_qs = 0 + jmax_qg = 0 + jmax_ni = 0 + jmax_nr = 0 + kmax_qc = 0 + kmax_qr = 0 + kmax_qi = 0 + kmax_qs = 0 + kmax_qg = 0 + kmax_ni = 0 + kmax_nr = 0 +! do i = 1, 256 +! mp_debug(i:i) = char(0) +! enddo +!$omp do schedule(guided) + j_loop: do j = j_start, j_end + i_loop: do i = i_start, i_end + + pptrain = 0. + pptsnow = 0. + pptgraul = 0. + pptice = 0. +! RAINNCV(i,j) = 0. +! IF ( PRESENT (snowncv) ) THEN +! SNOWNCV(i,j) = 0. +! ENDIF +! IF ( PRESENT (graupelncv) ) THEN +! GRAUPELNCV(i,j) = 0. +! ENDIF +! SR(i,j) = 0. + + do k = kts, kte + t1d(k) = th(i,k,j)*pii(i,k,j) + p1d(k) = p(i,k,j) + dz1d(k) = dz(i,k,j) + qv1d(k) = qv(i,k,j) + qc1d(k) = qc(i,k,j) + qi1d(k) = qi(i,k,j) + qr1d(k) = qr(i,k,j) + qs1d(k) = qs(i,k,j) + qg1d(k) = qg(i,k,j) + ni1d(k) = ni(i,k,j) + nr1d(k) = nr(i,k,j) + enddo + + call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, t1d, p1d, dz1d, & + pptrain, pptsnow, pptgraul, pptice, & + kts, kte, dt, i, j) + +! pcp_ra(i,j) = pptrain +! pcp_sn(i,j) = pptsnow +! pcp_gr(i,j) = pptgraul +! pcp_ic(i,j) = pptice +! RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice + RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + IF ( PRESENT(SNOWNC) ) SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice + IF ( PRESENT(SNOWNCV) ) SNOWNCV(i,j) = pptsnow + pptice + + IF ( PRESENT(GRAUPELNC) ) GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul + IF ( PRESENT(GRAUPELNCV) ) GRAUPELNCV(i,j) = pptgraul + + SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) + + do k = kts, kte + qv(i,k,j) = qv1d(k) + qc(i,k,j) = qc1d(k) + qi(i,k,j) = qi1d(k) + qr(i,k,j) = qr1d(k) + qs(i,k,j) = qs1d(k) + qg(i,k,j) = qg1d(k) + ni(i,k,j) = ni1d(k) + nr(i,k,j) = nr1d(k) + th(i,k,j) = t1d(k)/pii(i,k,j) +! if (qc1d(k) .gt. qc_max) then +! imax_qc = i +! jmax_qc = j +! kmax_qc = k +! qc_max = qc1d(k) +! elseif (qc1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (qr1d(k) .gt. qr_max) then +! imax_qr = i +! jmax_qr = j +! kmax_qr = k +! qr_max = qr1d(k) +! elseif (qr1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (nr1d(k) .gt. nr_max) then +! imax_nr = i +! jmax_nr = j +! kmax_nr = k +! nr_max = nr1d(k) +! elseif (nr1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (qs1d(k) .gt. qs_max) then +! imax_qs = i +! jmax_qs = j +! kmax_qs = k +! qs_max = qs1d(k) +! elseif (qs1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (qi1d(k) .gt. qi_max) then +! imax_qi = i +! jmax_qi = j +! kmax_qi = k +! qi_max = qi1d(k) +! elseif (qi1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (qg1d(k) .gt. qg_max) then +! imax_qg = i +! jmax_qg = j +! kmax_qg = k +! qg_max = qg1d(k) +! elseif (qg1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif +! if (ni1d(k) .gt. ni_max) then +! imax_ni = i +! jmax_ni = j +! kmax_ni = k +! ni_max = ni1d(k) +! elseif (ni1d(k) .lt. 0.0) then +! write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) +! endif + if (qv1d(k) .lt. 1.E-7) then + if (k.lt.kte-2 .and. k.gt.kts+1) then + qv(i,k,j) = 0.5*(qv(i,k-1,j) + qv(i,k+1,j)) + ! note, if qv(i,k+1,j) < 0 then qv(i,k,j) could still be < 0 + if (qv1d(k) .lt. 1.E-7) then + qv(i,k,j) = 1.E-7 + endif + else + qv(i,k,j) = 1.E-7 + endif +! write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), & +! ' at i,j,k=', i,j,k + ! CALL wrf_debug(150, mp_debug) + endif + enddo + +! IF ( PRESENT (diagflag) ) THEN +! if (diagflag .and. do_radar_ref == 1) then +! call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & +! t1d, p1d, dBZ, kts, kte, i, j) +! do k = kts, kte +! refl_10cm(i,k,j) = MAX(-35., dBZ(k)) +! enddo +! endif +! ENDIF + + enddo i_loop + enddo j_loop + !$omp end do + !$omp end parallel + +! DEBUG - GT +! write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & +! 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & +! 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & +! 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & +! 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & +! 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & +! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & +! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' + ! CALL wrf_debug(150, mp_debug) +! END DEBUG - GT + +! do i = 1, 256 +! mp_debug(i:i) = char(0) +! enddo + + END SUBROUTINE mp_gt_driver + +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +!.. This subroutine computes the moisture tendencies of water vapor, +!.. cloud droplets, rain, cloud ice (pristine), snow, and graupel. +!.. Previously this code was based on Reisner et al (1998), but few of +!.. those pieces remain. A complete description is now found in +!.. Thompson et al. (2004, 2008). +!+---+-----------------------------------------------------------------+ +! + subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, t1d, p1d, dzq, & + pptrain, pptsnow, pptgraul, pptice, & + kts, kte, dt, ii, jj) + + implicit none + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte, ii, jj + REAL, DIMENSION(kts:kte), INTENT(INOUT):: & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, t1d, p1d + REAL, DIMENSION(kts:kte), INTENT(IN):: dzq + REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice + REAL, INTENT(IN):: dt + +!..Local variables + REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & + qrten, qsten, qgten, niten, nrten + + DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd + + DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, & + prr_rcg, prr_sml, prr_gml, & + prr_rci, prv_rev, & + pnr_wau, pnr_rcs, pnr_rcg, & + pnr_rci, pnr_sml, pnr_gml, & + pnr_rev, pnr_rcr, pnr_rfz + + DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, & + pni_ihm, pri_wfz, pni_wfz, & + pri_rfz, pni_rfz, pri_ide, & + pni_ide, pri_rci, pni_rci, & + pni_sci, pni_iau + + DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, & + prs_scw, prs_sde, prs_ihm, & + prs_ide + + DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, & + prg_gcw, prg_rci, prg_rcs, & + prg_rcg, prg_ihm + + DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 + + REAL, DIMENSION(kts:kte):: temp, pres, qv + REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr + REAL, DIMENSION(kts:kte):: rho, rhof, rhof2 + REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs + REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati + REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, & + tcond, lvap, ocp, lvt2 + + DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g + REAL, DIMENSION(kts:kte):: mvd_r, mvd_c + REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, & + smoc, smod, smoe, smof + + REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n + + REAL:: rgvm, delta_tp, orho, lfus2 + REAL, DIMENSION(4):: onstep + DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg + DOUBLE PRECISION:: lami, ilami + REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m + DOUBLE PRECISION:: Dr_star + REAL:: zeta1, zeta, taud, tau + REAL:: stoke_r, stoke_s, stoke_g, stoke_i + REAL:: vti, vtr, vts, vtg + REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk + REAL, DIMENSION(kts:kte):: vts_boost + REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow + REAL:: a_, b_, loga_, A1, A2, tf + REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat + REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr + REAL:: xsat, rate_max, sump, ratio + REAL:: clap, fcd, dfcd + REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl + REAL:: r_frac, g_frac + REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr + REAL:: dtsave, odts, odt, odzq + REAL:: xslw1, ygra1, zans1 + INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq + INTEGER, DIMENSION(4):: ksed1 + INTEGER:: nir, nis, nig, nii, nic + INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, & + idx_i1, idx_i, idx_c, idx, idx_d + LOGICAL:: melti, no_micro + LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg + LOGICAL:: debug_flag + +!+---+ + debug_flag = .false. +! if (ii.eq.315 .and. jj.eq.2) debug_flag = .true. + + no_micro = .true. + dtsave = dt + odt = 1./dt + odts = 1./dtsave + iexfrq = 1 + +!+---+-----------------------------------------------------------------+ +!.. Source/sink terms. First 2 chars: "pr" represents source/sink of +!.. mass while "pn" represents source/sink of number. Next char is one +!.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for +!.. cloud water, "s" for snow, and "g" for graupel. Next chars +!.. represent processes: "de" for sublimation/deposition, "ev" for +!.. evaporation, "fz" for freezing, "ml" for melting, "au" for +!.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop +!.. secondary ice production, and "c" for collection followed by the +!.. character for the species being collected. ALL of these terms are +!.. positive (except for deposition/sublimation terms which can switch +!.. signs based on super/subsaturation) and are treated as negatives +!.. where necessary in the tendency equations. +!+---+-----------------------------------------------------------------+ + + do k = kts, kte + tten(k) = 0. + qvten(k) = 0. + qcten(k) = 0. + qiten(k) = 0. + qrten(k) = 0. + qsten(k) = 0. + qgten(k) = 0. + niten(k) = 0. + nrten(k) = 0. + + prw_vcd(k) = 0. + + prv_rev(k) = 0. + prr_wau(k) = 0. + prr_rcw(k) = 0. + prr_rcs(k) = 0. + prr_rcg(k) = 0. + prr_sml(k) = 0. + prr_gml(k) = 0. + prr_rci(k) = 0. + pnr_wau(k) = 0. + pnr_rcs(k) = 0. + pnr_rcg(k) = 0. + pnr_rci(k) = 0. + pnr_sml(k) = 0. + pnr_gml(k) = 0. + pnr_rev(k) = 0. + pnr_rcr(k) = 0. + pnr_rfz(k) = 0. + + pri_inu(k) = 0. + pni_inu(k) = 0. + pri_ihm(k) = 0. + pni_ihm(k) = 0. + pri_wfz(k) = 0. + pni_wfz(k) = 0. + pri_rfz(k) = 0. + pni_rfz(k) = 0. + pri_ide(k) = 0. + pni_ide(k) = 0. + pri_rci(k) = 0. + pni_rci(k) = 0. + pni_sci(k) = 0. + pni_iau(k) = 0. + + prs_iau(k) = 0. + prs_sci(k) = 0. + prs_rcs(k) = 0. + prs_scw(k) = 0. + prs_sde(k) = 0. + prs_ihm(k) = 0. + prs_ide(k) = 0. + + prg_scw(k) = 0. + prg_rfz(k) = 0. + prg_gde(k) = 0. + prg_gcw(k) = 0. + prg_rci(k) = 0. + prg_rcs(k) = 0. + prg_rcg(k) = 0. + prg_ihm(k) = 0. + enddo + +!+---+-----------------------------------------------------------------+ +!..Put column of data into local arrays. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + temp(k) = t1d(k) + qv(k) = MAX(1.E-10, qv1d(k)) + pres(k) = p1d(k) + rho(k) = 0.622*pres(k)/(RR2*temp(k)*(qv(k)+0.622)) + if (qc1d(k) .gt. R1) then + no_micro = .false. + rc(k) = qc1d(k)*rho(k) + L_qc(k) = .true. + else + qc1d(k) = 0.0 + rc(k) = R1 + L_qc(k) = .false. + endif + if (qi1d(k) .gt. R1) then + no_micro = .false. + ri(k) = qi1d(k)*rho(k) + ni(k) = MAX(R2, ni1d(k)*rho(k)) + L_qi(k) = .true. + lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi + ilami = 1./lami + xDi = (bm_i + mu_i + 1.) * ilami + if (xDi.lt. 20.E-6) then + lami = cie(2)/20.E-6 + ni(k) = MIN(250.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + elseif (xDi.gt. 300.E-6) then + lami = cie(2)/300.E-6 + ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i + endif + else + qi1d(k) = 0.0 + ni1d(k) = 0.0 + ri(k) = R1 + ni(k) = R2 + L_qi(k) = .false. + endif + + mvd_r(k) = 0.0 ! must be initialized or a later test can crash where qr1d(k)<=R1 + if (qr1d(k) .gt. R1) then + no_micro = .false. + rr(k) = qr1d(k)*rho(k) + nr(k) = MAX(R2, nr1d(k)*rho(k)) + L_qr(k) = .true. + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + if (mvd_r(k) .gt. 2.5E-3) then + mvd_r(k) = 2.5E-3 + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r + elseif (mvd_r(k) .lt. D0r*0.75) then + mvd_r(k) = D0r*0.75 + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r + endif + else + qr1d(k) = 0.0 + nr1d(k) = 0.0 + rr(k) = R1 + nr(k) = R2 + L_qr(k) = .false. + endif + if (qs1d(k) .gt. R1) then + no_micro = .false. + rs(k) = qs1d(k)*rho(k) + L_qs(k) = .true. + else + qs1d(k) = 0.0 + rs(k) = R1 + L_qs(k) = .false. + endif + if (qg1d(k) .gt. R1) then + no_micro = .false. + rg(k) = qg1d(k)*rho(k) + L_qg(k) = .true. + else + qg1d(k) = 0.0 + rg(k) = R1 + L_qg(k) = .false. + endif + enddo + + +!+---+-----------------------------------------------------------------+ +!..Derive various thermodynamic variables frequently used. +!.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from +!.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from +!.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + tempc = temp(k) - 273.15 + rhof(k) = SQRT(RHO_NOT/rho(k)) + rhof2(k) = SQRT(rhof(k)) + qvs(k) = rslf(pres(k), temp(k)) + delQvs(k) = MAX(0.0, rslf(pres(k), 273.15)-qv(k)) + if (tempc .le. 0.0) then + qvsi(k) = rsif(pres(k), temp(k)) + else + qvsi(k) = qvs(k) + endif + satw(k) = qv(k)/qvs(k) + sati(k) = qv(k)/qvsi(k) + ssatw(k) = satw(k) - 1. + ssati(k) = sati(k) - 1. + if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0 + if (abs(ssati(k)).lt. eps) ssati(k) = 0.0 + if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false. + diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k)) + if (tempc .ge. 0.0) then + visco(k) = (1.718+0.0049*tempc)*1.0E-5 + else + visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5 + endif + ocp(k) = 1./(Cp2*(1.+0.887*qv(k))) + vsc2(k) = SQRT(rho(k)/visco(k)) + lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc + tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936 + enddo + +!+---+-----------------------------------------------------------------+ +!..If no existing hydrometeor species and no chance to initiate ice or +!.. condense cloud water, just exit quickly! +!+---+-----------------------------------------------------------------+ + + if (no_micro) return + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope, and useful moments for snow. +!+---+-----------------------------------------------------------------+ + if (.not. iiwarm) then + do k = kts, kte + if (.not. L_qs(k)) CYCLE + tc0 = MIN(-0.1, temp(k)-273.15) + smob(k) = rs(k)*oams + +!..All other moments based on reference, 2nd moment. If bm_s.ne.2, +!.. then we must compute actual 2nd moment and use as reference. + if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then + smo2(k) = smob(k) + else + loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & + + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & + + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & + + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & + + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & + + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & + + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & + + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & + + sb(10)*bm_s*bm_s*bm_s + smo2(k) = (smob(k)/a_)**(1./b_) + endif + +!..Calculate 0th moment. Represents snow number concentration. + loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0 + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0 + smo0(k) = a_ * smo2(k)**b_ + +!..Calculate 1st moment. Useful for depositional growth and melting. + loga_ = sa(1) + sa(2)*tc0 + sa(3) & + + sa(4)*tc0 + sa(5)*tc0*tc0 & + + sa(6) + sa(7)*tc0*tc0 & + + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 & + + sa(10) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 & + + sb(5)*tc0*tc0 + sb(6) & + + sb(7)*tc0*tc0 + sb(8)*tc0 & + + sb(9)*tc0*tc0*tc0 + sb(10) + smo1(k) = a_ * smo2(k)**b_ + +!..Calculate bm_s+1 (th) moment. Useful for diameter calcs. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) + smoc(k) = a_ * smo2(k)**b_ + +!..Calculate bv_s+2 (th) moment. Useful for riming. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) & + + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 & + + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) & + + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(13)*cse(13)*cse(13) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) & + + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) & + + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13) + smoe(k) = a_ * smo2(k)**b_ + +!..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth. + loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) & + + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 & + + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) & + + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 & + + sa(10)*cse(16)*cse(16)*cse(16) + a_ = 10.0**loga_ + b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) & + + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) & + + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) & + + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16) + smof(k) = a_ * smo2(k)**b_ + + enddo + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for graupel. +!+---+-----------------------------------------------------------------+ + N0_min = gonv_max + do k = kte, kts, -1 + if (temp(k).lt.270.65 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then + xslw1 = 4.01 + alog10(mvd_r(k)) + else + xslw1 = 0.01 + endif + ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) + zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + N0_exp = 10.**(zans1) + N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + N0_min = MIN(N0_exp, N0_min) + N0_exp = N0_min + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ilamg(k) = 1./lamg + N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) +!+---+-----------------------------------------------------------------+ +! if( debug_flag .and. k.lt.42) then +! if (k.eq.41) write(mp_debug,*) 'DEBUG-GT: K, zans1, rc, rr, rg, N0_g' +! if (k.eq.41) CALL wrf_debug(0, mp_debug) +! write(mp_debug, 'a, i2, 1x, f6.3, 1x, 4(1x,e13.6,1x)') & +! ' GT ', k, zans1, rc(k), rr(k), rg(k), N0_g(k) +! CALL wrf_debug(0, mp_debug) +! endif +!+---+-----------------------------------------------------------------+ + enddo + + endif + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for rain. +!+---+-----------------------------------------------------------------+ + do k = kte, kts, -1 + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + ilamr(k) = 1./lamr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + N0_r(k) = nr(k)*org2*lamr**cre(2) + enddo + +!+---+-----------------------------------------------------------------+ +!..Compute warm-rain process terms (except evap done later). +!+---+-----------------------------------------------------------------+ + + do k = kts, kte + +!..Rain self-collection follows Seifert, 1994 and drop break-up +!.. follows Verlinde and Cotton, 1993. RAIN2M + if (L_qr(k) .and. mvd_r(k).gt. D0r) then +!-GT Ef_rr = 1.0 +!-GT if (mvd_r(k) .gt. 1500.0E-6) then + Ef_rr = 2.0 - EXP(2300.0*(mvd_r(k)-1600.0E-6)) +!-GT endif + pnr_rcr(k) = Ef_rr * 4.*nr(k)*rr(k) + endif + + mvd_c(k) = D0c + if (.not. L_qc(k)) CYCLE + xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*Nt_c))**obmr) * 1.E6) + lamc = (Nt_c*am_r* ccg(2) * ocg1 / rc(k))**obmr + mvd_c(k) = (3.0+mu_c+0.672) / lamc + +!..Autoconversion follows Berry & Reinhardt (1974) with characteristic +!.. diameters correctly computed from gamma distrib of cloud droplets. + if (rc(k).gt. 0.01e-3) then + Dc_g = ((ccg(3)*ocg2)**obmr / lamc) * 1.E6 + Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) & + **(1./6.) + zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) & + + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4)) + zeta = 0.027*rc(k)*zeta1 + taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1 + tau = 3.72/(rc(k)*taud) + prr_wau(k) = zeta/tau + prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) + pnr_wau(k) = prr_wau(k) / (am_r*mu_c*D0r*D0r*D0r) ! RAIN2M + endif + +!..Rain collecting cloud water. In CE, assume Dc<1). Either way, only bother to do sedimentation below +!.. 1st level that contains any sedimenting particles (k=ksed1 on down). +!.. New in v3.0+ is computing separate for rain, ice, snow, and +!.. graupel species thus making code faster with credit to J. Schmidt. +!+---+-----------------------------------------------------------------+ + nstep = 0 + onstep(:) = 1.0 + ksed1(:) = 1 + do k = kte+1, kts, -1 + vtrk(k) = 0. + vtnrk(k) = 0. + vtik(k) = 0. + vtnik(k) = 0. + vtsk(k) = 0. + vtgk(k) = 0. + enddo + do k = kte, kts, -1 + vtr = 0. + rhof(k) = SQRT(RHO_NOT/rho(k)) + + if (rr(k).gt. R1) then + lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr + vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) & + *((lamr+fv_r)**(-cre(6))) + vtrk(k) = vtr +! First below is technically correct: +! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) & +! *((lamr+fv_r)**(-cre(5))) +! Test: make number fall faster (but still slower than mass) +! Goal: less prominent size sorting + vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) & + *((lamr+fv_r)**(-cre(7))) + vtnrk(k) = vtr + else + vtrk(k) = vtrk(k+1) + vtnrk(k) = vtnrk(k+1) + endif + + if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then + ksed1(1) = MAX(ksed1(1), k) + delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k))) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(1) .eq. kte) ksed1(1) = kte-1 + if (nstep .gt. 0) onstep(1) = 1./REAL(nstep) + +!+---+-----------------------------------------------------------------+ + + if (.not. iiwarm) then + + nstep = 0 + do k = kte, kts, -1 + vti = 0. + + if (ri(k).gt. R1) then + lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi + ilami = 1./lami + vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i + vtik(k) = vti +! First below is technically correct: +! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i +! Goal: less prominent size sorting + vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i + vtnik(k) = vti + else + vtik(k) = vtik(k+1) + vtnik(k) = vtnik(k+1) + endif + + if (vtik(k) .gt. 1.E-3) then + ksed1(2) = MAX(ksed1(2), k) + delta_tp = dzq(k)/vtik(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(2) .eq. kte) ksed1(2) = kte-1 + if (nstep .gt. 0) onstep(2) = 1./REAL(nstep) + +!+---+-----------------------------------------------------------------+ + + nstep = 0 + do k = kte, kts, -1 + vts = 0. + + if (rs(k).gt. R1) then + xDs = smoc(k) / smob(k) + Mrat = 1./xDs + ils1 = 1./(Mrat*Lam0 + fv_s) + ils2 = 1./(Mrat*Lam1 + fv_s) + t1_vts = Kap0*csg(4)*ils1**cse(4) + t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10) + ils1 = 1./(Mrat*Lam0) + ils2 = 1./(Mrat*Lam1) + t3_vts = Kap0*csg(1)*ils1**cse(1) + t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7) + vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) + if (temp(k).gt. T_0) then + vtsk(k) = MAX(vts*vts_boost(k), vtrk(k)) + else + vtsk(k) = vts*vts_boost(k) + endif + else + vtsk(k) = vtsk(k+1) + endif + + if (vtsk(k) .gt. 1.E-3) then + ksed1(3) = MAX(ksed1(3), k) + delta_tp = dzq(k)/vtsk(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(3) .eq. kte) ksed1(3) = kte-1 + if (nstep .gt. 0) onstep(3) = 1./REAL(nstep) + +!+---+-----------------------------------------------------------------+ + + nstep = 0 + do k = kte, kts, -1 + vtg = 0. + + if (rg(k).gt. R1) then + vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g + if (temp(k).gt. T_0) then + vtgk(k) = MAX(vtg, vtrk(k)) + else + vtgk(k) = vtg + endif + else + vtgk(k) = vtgk(k+1) + endif + + if (vtgk(k) .gt. 1.E-3) then + ksed1(4) = MAX(ksed1(4), k) + delta_tp = dzq(k)/vtgk(k) + nstep = MAX(nstep, INT(DT/delta_tp + 1.)) + endif + enddo + if (ksed1(4) .eq. kte) ksed1(4) = kte-1 + if (nstep .gt. 0) onstep(4) = 1./REAL(nstep) + endif + +!+---+-----------------------------------------------------------------+ +!..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD, +!.. whereas neglect m(D) term for number concentration. Therefore, +!.. cloud ice has proper differential sedimentation. +!.. New in v3.0+ is computing separate for rain, ice, snow, and +!.. graupel species thus making code faster with credit to J. Schmidt. +!+---+-----------------------------------------------------------------+ + + nstep = NINT(1./onstep(1)) + do n = 1, nstep + do k = kte, kts, -1 + sed_r(k) = vtrk(k)*rr(k) + sed_n(k) = vtnrk(k)*nr(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho + nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho + rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1)) + nr(k) = MAX(R2, nr(k) - sed_n(k)*odzq*DT*onstep(1)) + do k = ksed1(1), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) & + *odzq*onstep(1)*orho + nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*onstep(1)*orho + rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) & + *odzq*DT*onstep(1)) + nr(k) = MAX(R2, nr(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*DT*onstep(1)) + enddo + + if (rr(kts).gt.R1*10.) & + pptrain = pptrain + sed_r(kts)*DT*onstep(1) + enddo + +!+---+-----------------------------------------------------------------+ + + nstep = NINT(1./onstep(2)) + do n = 1, nstep + do k = kte, kts, -1 + sed_i(k) = vtik(k)*ri(k) + sed_n(k) = vtnik(k)*ni(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho + niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho + ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2)) + ni(k) = MAX(R2, ni(k) - sed_n(k)*odzq*DT*onstep(2)) + do k = ksed1(2), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) & + *odzq*onstep(2)*orho + niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*onstep(2)*orho + ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) & + *odzq*DT*onstep(2)) + ni(k) = MAX(R2, ni(k) + (sed_n(k+1)-sed_n(k)) & + *odzq*DT*onstep(2)) + enddo + + if (ri(kts).gt.R1*10.) & + pptice = pptice + sed_i(kts)*DT*onstep(2) + enddo + +!+---+-----------------------------------------------------------------+ + + nstep = NINT(1./onstep(3)) + do n = 1, nstep + do k = kte, kts, -1 + sed_s(k) = vtsk(k)*rs(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho + rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3)) + do k = ksed1(3), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) & + *odzq*onstep(3)*orho + rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) & + *odzq*DT*onstep(3)) + enddo + + if (rs(kts).gt.R1*10.) & + pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) + enddo + +!+---+-----------------------------------------------------------------+ + + nstep = NINT(1./onstep(4)) + do n = 1, nstep + do k = kte, kts, -1 + sed_g(k) = vtgk(k)*rg(k) + enddo + k = kte + odzq = 1./dzq(k) + orho = 1./rho(k) + qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho + rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4)) + do k = ksed1(4), kts, -1 + odzq = 1./dzq(k) + orho = 1./rho(k) + qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) & + *odzq*onstep(4)*orho + rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) & + *odzq*DT*onstep(4)) + enddo + + if (rg(kts).gt.R1*10.) & + pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) + enddo + +!+---+-----------------------------------------------------------------+ +!.. Instantly melt any cloud ice into cloud water if above 0C and +!.. instantly freeze any cloud water found below HGFR. +!+---+-----------------------------------------------------------------+ + if (.not. iiwarm) then + do k = kts, kte + xri = MAX(0.0, qi1d(k) + qiten(k)*DT) + if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then + qcten(k) = qcten(k) + xri*odt + qiten(k) = qiten(k) - xri*odt + niten(k) = -ni1d(k)*odt + tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY) + endif + + xrc = MAX(0.0, qc1d(k) + qcten(k)*DT) + if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then + lfus2 = lsub - lvap(k) + qiten(k) = qiten(k) + xrc*odt + niten(k) = niten(k) + xrc/xm0i * odt + qcten(k) = qcten(k) - xrc*odt + tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY) + endif + enddo + endif + +!+---+-----------------------------------------------------------------+ +!.. All tendencies computed, apply and pass back final values to parent. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + t1d(k) = t1d(k) + tten(k)*DT + qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) + qc1d(k) = qc1d(k) + qcten(k)*DT + if (qc1d(k) .le. R1) qc1d(k) = 0.0 + qi1d(k) = qi1d(k) + qiten(k)*DT + ni1d(k) = MAX(R2/rho(k), ni1d(k) + niten(k)*DT) + if (qi1d(k) .le. R1) then + qi1d(k) = 0.0 + ni1d(k) = 0.0 + else + lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi + ilami = 1./lami + xDi = (bm_i + mu_i + 1.) * ilami + if (xDi.lt. 20.E-6) then + lami = cie(2)/20.E-6 + elseif (xDi.gt. 300.E-6) then + lami = cie(2)/300.E-6 + endif + ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & + 250.D3/rho(k)) + endif + qr1d(k) = qr1d(k) + qrten(k)*DT + nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) + if (qr1d(k) .le. R1) then + qr1d(k) = 0.0 + nr1d(k) = 0.0 + else + lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr + mvd_r(k) = (3.0 + mu_r + 0.672) / lamr + if (mvd_r(k) .gt. 2.5E-3) then + mvd_r(k) = 2.5E-3 + elseif (mvd_r(k) .lt. D0r*0.75) then + mvd_r(k) = D0r*0.75 + endif + lamr = (3.0 + mu_r + 0.672) / mvd_r(k) + nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r + endif + qs1d(k) = qs1d(k) + qsten(k)*DT + if (qs1d(k) .le. R1) qs1d(k) = 0.0 + qg1d(k) = qg1d(k) + qgten(k)*DT + if (qg1d(k) .le. R1) qg1d(k) = 0.0 + enddo + + end subroutine mp_thompson +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Creation of the lookup tables and support functions found below here. +!+---+-----------------------------------------------------------------+ +!..Rain collecting graupel (and inverse). Explicit CE integration. +!+---+-----------------------------------------------------------------+ + + subroutine qr_acr_qg + + implicit none + +!..Local variables + INTEGER:: i, j, k, m, n, n2 + INTEGER:: km, km_s, km_e + DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g + DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r + DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr + DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2 + logical :: lexist, lopen + integer :: good + + good = 0 + + + INQUIRE(FILE="qr_acr_qg_mpt.dat",EXIST=lexist) + IF ( lexist ) THEN + if (this_image()==1) print *, "ThompMP: read qr_acr_qg_mpt.dat instead of computing" + OPEN(63,file="qr_acr_qg_mpt.dat",form="unformatted",err=1234) + READ(63,err=1234) tcg_racg + READ(63,err=1234) tmr_racg + READ(63,err=1234) tcr_gacr + READ(63,err=1234) tmg_gacr + READ(63,err=1234) tnr_racg + READ(63,err=1234) tnr_gacr + good = 1 + 1234 CONTINUE + INQUIRE(63,opened=lopen) + IF (lopen) THEN + CLOSE(63) + ENDIF + + ENDIF + + if (good .NE. 1) then + + do n2 = 1, nbr + ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) + vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & + + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) & + - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2) + enddo + do n = 1, nbg + vg(n) = av_g*Dg(n)**bv_g + enddo + + !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for + !.. fortran indices. J. Michalakes, 2009Oct30. + + ! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + ! CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) + ! #else + km_s = 0 + km_e = ntb_r*ntb_r1 - 1 + ! #endif + + !$omp parallel default(shared) & + !$omp private(km,i,j,k,m,n,n2,lam_exp,lamr,N0_r,N_r,lamg,N_g,N0_g,t1,t2,z1,z2,y1,y2,massr,massg,dvg,dvr) & + !$omp shared(tcg_racg,tmr_racg,tcr_gacr,tmg_gacr,tnr_racg,tnr_gacr,mu_r,ef_rs,am_s) & + !$omp firstprivate(crg,cre,cge,cgg,ore1,oge1,org1,org2,ogg2,ogg1,obmg,obmr) & + !$omp firstprivate(km_s,km_e,Dr,Dg,dtg,dtr,vr,vg) + !$omp do + do km = km_s, km_e + m = km / ntb_r1 + 1 + k = mod( km , ntb_r1 ) + 1 + + lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2) + do n2 = 1, nbr + N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2) + enddo + + do j = 1, ntb_g + do i = 1, ntb_g1 + lam_exp = (N0g_exp(i)*am_g*cgg(1)/r_g(j))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + N0_g = N0g_exp(i)/(cgg(2)*lam_exp) * lamg**cge(2) + do n = 1, nbg + N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n) + enddo + + t1 = 0.0d0 + t2 = 0.0d0 + z1 = 0.0d0 + z2 = 0.0d0 + y1 = 0.0d0 + y2 = 0.0d0 + do n2 = 1, nbr + massr = am_r * Dr(n2)**bm_r + do n = 1, nbg + massg = am_g * Dg(n)**bm_g + + dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n))) + dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2))) + + t1 = t1+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg*massg * N_g(n)* N_r(n2) + z1 = z1+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg*massr * N_g(n)* N_r(n2) + y1 = y1+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg * N_g(n)* N_r(n2) + + t2 = t2+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr*massr * N_g(n)* N_r(n2) + y2 = y2+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr * N_g(n)* N_r(n2) + z2 = z2+ PI2*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr*massg * N_g(n)* N_r(n2) + enddo + 97 continue + enddo + tcg_racg(i,j,k,m) = t1 + tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) + tcr_gacr(i,j,k,m) = t2 + tmg_gacr(i,j,k,m) = z2 + tnr_racg(i,j,k,m) = y1 + tnr_gacr(i,j,k,m) = y2 + enddo + enddo + enddo + !$omp end do + !$omp end parallel + IF ( this_image()==1 ) THEN + print *, "Writing qr_acr_qg_mpt.dat in Thompson MP init" + OPEN(63,file="qr_acr_qg_mpt.dat",form="unformatted",err=9234) + WRITE(63,err=9234) tcg_racg + WRITE(63,err=9234) tmr_racg + WRITE(63,err=9234) tcr_gacr + WRITE(63,err=9234) tmg_gacr + WRITE(63,err=9234) tnr_racg + WRITE(63,err=9234) tnr_gacr + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + print *, ("Error writing qr_acr_qg_mpt.dat") + ENDIF + + endif + +!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). + +! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + ! CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) +! #endif + + + end subroutine qr_acr_qg +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Rain collecting snow (and inverse). Explicit CE integration. +!+---+-----------------------------------------------------------------+ + + subroutine qr_acr_qs + + implicit none +!..Local variables + INTEGER:: i, j, k, m, n, n2 + INTEGER:: km, km_s, km_e + DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r + DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s + DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3 + DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2 + DOUBLE PRECISION:: dvs, dvr, masss, massr + DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4 + DOUBLE PRECISION:: y1, y2, y3, y4 + LOGICAL lexist,lopen + INTEGER :: good + ! LOGICAL, EXTERNAL :: wrf_dm_on_monitor + +!+---+ + + ! CALL nl_get_force_read_thompson(1,force_read_thompson) + ! CALL nl_get_write_thompson_tables(1,write_thompson_tables) + + good = 0 + INQUIRE(FILE="qr_acr_qs_mpt.dat",EXIST=lexist) + IF ( lexist ) THEN + IF ( this_image() == 1 ) print *, "ThompMP: read qr_acr_qs_mpt.dat instead of computing" + OPEN(63,file="qr_acr_qs_mpt.dat",form="unformatted",err=1234) + READ(63,err=1234)tcs_racs1 + READ(63,err=1234)tmr_racs1 + READ(63,err=1234)tcs_racs2 + READ(63,err=1234)tmr_racs2 + READ(63,err=1234)tcr_sacr1 + READ(63,err=1234)tms_sacr1 + READ(63,err=1234)tcr_sacr2 + READ(63,err=1234)tms_sacr2 + READ(63,err=1234)tnr_racs1 + READ(63,err=1234)tnr_racs2 + READ(63,err=1234)tnr_sacr1 + READ(63,err=1234)tnr_sacr2 + good = 1 + 1234 CONTINUE + INQUIRE(63,opened=lopen) + IF (lopen) THEN + CLOSE(63) + ENDIF + ENDIF + + if (good .NE. 1) then + if (this_image()==1) print *, "ThompMP: computing qr_acr_qs" + !+---+ + do n2 = 1, nbr + ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) + vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & + + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) & + - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2) + D1(n2) = (vr(n2)/av_s)**(1./bv_s) + enddo + do n = 1, nbs + vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n)) + enddo + + !..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for + !.. fortran indices. J. Michalakes, 2009Oct30. + + ! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + ! CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) + ! #else + km_s = 0 + km_e = ntb_r*ntb_r1 - 1 + ! #endif + + !$omp parallel default(private) & + !$omp shared(tcs_racs1,tmr_racs1,tcs_racs2,tmr_racs2,tcr_sacr1,tms_sacr1)& + !$omp shared(tcr_sacr2,tms_sacr2,tnr_racs1,tnr_racs2,tnr_sacr1,tnr_sacr2,mu_r,ef_rs,am_s) & + !$omp private(second,a_,b_,n2,km,M2,M3,oM3,Mrat,M0,slam1,slam2,i,j,k,m,n,lam_exp,lamr,N0_r,N_r,loga_) & + !$omp private(t1,t2,t3,t4,z1,z2,z3,z4,y1,y2,y3,y4,N_s,massr,masss,dvs,dvr) & + !$omp firstprivate(km_s,km_e,vr,vs) & + !$omp firstprivate(Dr,dtr,dts,Ds,crg,cre,cse,ore1,org2,org1,obmr,oams) + !$omp do + do km = km_s, km_e + m = km / ntb_r1 + 1 + k = mod( km , ntb_r1 ) + 1 + + lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2) + do n2 = 1, nbr + N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2) + enddo + + do j = 1, ntb_t + do i = 1, ntb_s + + !..From the bm_s moment, compute plus one moment. If we are not + !.. using bm_s=2, then we must transform to the pure 2nd moment + !.. (variable called "second") and then to the bm_s+1 moment. + + M2 = r_s(i)*oams *1.0d0 + if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then + loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s & + + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) & + + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s & + + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) & + + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s & + + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) & + + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s & + + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) & + + sb(10)*bm_s*bm_s*bm_s + second = (M2/a_)**(1./b_) + else + second = M2 + endif + + loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) & + + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) & + + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) & + + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) & + + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) & + + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) & + + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1) + M3 = a_ * second**b_ + + oM3 = 1./M3 + Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3) + M0 = (M2*oM3)**mu_s + slam1 = M2 * oM3 * Lam0 + slam2 = M2 * oM3 * Lam1 + + do n = 1, nbs + N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) & + + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n) + enddo + + t1 = 0.0d0 + t2 = 0.0d0 + t3 = 0.0d0 + t4 = 0.0d0 + z1 = 0.0d0 + z2 = 0.0d0 + z3 = 0.0d0 + z4 = 0.0d0 + y1 = 0.0d0 + y2 = 0.0d0 + y3 = 0.0d0 + y4 = 0.0d0 + do n2 = 1, nbr + massr = am_r * Dr(n2)**bm_r + do n = 1, nbs + masss = am_s * Ds(n)**bm_s + + dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n))) + dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2))) + + if (massr .gt. 1.5*masss) then + t1 = t1+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*masss * N_s(n)* N_r(n2) + z1 = z1+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*massr * N_s(n)* N_r(n2) + y1 = y1+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs * N_s(n)* N_r(n2) + else + t3 = t3+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*masss * N_s(n)* N_r(n2) + z3 = z3+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*massr * N_s(n)* N_r(n2) + y3 = y3+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs * N_s(n)* N_r(n2) + endif + + if (massr .gt. 1.5*masss) then + t2 = t2+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*massr * N_s(n)* N_r(n2) + y2 = y2+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr * N_s(n)* N_r(n2) + z2 = z2+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*masss * N_s(n)* N_r(n2) + else + t4 = t4+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*massr * N_s(n)* N_r(n2) + y4 = y4+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr * N_s(n)* N_r(n2) + z4 = z4+ PI2*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*masss * N_s(n)* N_r(n2) + endif + + enddo + enddo + tcs_racs1(i,j,k,m) = t1 + tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) + tcs_racs2(i,j,k,m) = t3 + tmr_racs2(i,j,k,m) = z3 + tcr_sacr1(i,j,k,m) = t2 + tms_sacr1(i,j,k,m) = z2 + tcr_sacr2(i,j,k,m) = t4 + tms_sacr2(i,j,k,m) = z4 + tnr_racs1(i,j,k,m) = y1 + tnr_racs2(i,j,k,m) = y3 + tnr_sacr1(i,j,k,m) = y2 + tnr_sacr2(i,j,k,m) = y4 + enddo + enddo + enddo + !$omp end do + !$omp end parallel + + IF ( this_image()==1 ) THEN + print *, "Writing qr_acr_qs_mpt.dat in Thompson MP init" + OPEN(63,file="qr_acr_qs_mpt.dat",form="unformatted",err=9234) + WRITE(63,err=9234)tcs_racs1 + WRITE(63,err=9234)tmr_racs1 + WRITE(63,err=9234)tcs_racs2 + WRITE(63,err=9234)tmr_racs2 + WRITE(63,err=9234)tcr_sacr1 + WRITE(63,err=9234)tms_sacr1 + WRITE(63,err=9234)tcr_sacr2 + WRITE(63,err=9234)tms_sacr2 + WRITE(63,err=9234)tnr_racs1 + WRITE(63,err=9234)tnr_racs2 + WRITE(63,err=9234)tnr_sacr1 + WRITE(63,err=9234)tnr_sacr2 + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + print*, "Error writing qr_acr_qs_mpt.dat" + ENDIF + ENDIF + +!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). + +! #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + ! CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + ! CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) +! #endif + + + end subroutine qr_acr_qs +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..This is a literal adaptation of Bigg (1954) probability of drops of +!..a particular volume freezing. Given this probability, simply freeze +!..the proportion of drops summing their masses. +!+---+-----------------------------------------------------------------+ + + subroutine freezeH2O + + implicit none +!..Local variables + INTEGER:: i, j, k, n, n2 + DOUBLE PRECISION, DIMENSION(nbr):: N_r, massr + DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc + DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, & + prob, vol, Texp, orho_w, & + lam_exp, lamr, N0_r, lamc, N0_c, y + LOGICAL lexist,lopen + INTEGER :: good + ! LOGICAL, EXTERNAL :: wrf_dm_on_monitor + +!+---+ + + good = 0 + ! IF ( this_image() == 1 ) THEN + INQUIRE(FILE="freezeH2O_mpt.dat",EXIST=lexist) + IF ( lexist ) THEN + IF ( this_image() == 1 ) print *, "ThompMP: read freezeH2O_mpt.dat instead of computing" + OPEN(63,file="freezeH2O_mpt.dat",form="unformatted",err=1234) + READ(63,err=1234)tpi_qrfz + READ(63,err=1234)tni_qrfz + READ(63,err=1234)tpg_qrfz + READ(63,err=1234)tnr_qrfz + READ(63,err=1234)tpi_qcfz + READ(63,err=1234)tni_qcfz + good = 1 + 1234 CONTINUE + INQUIRE(63,opened=lopen) + IF (lopen) THEN + CLOSE(63) + endif + ENDIF + + + IF ( good .NE. 1 ) THEN + if (this_image()==1) print *, "ThompMP: computing freezeH2O" + !+---+ + orho_w = 1./rho_w + + do n2 = 1, nbr + massr(n2) = am_r*Dr(n2)**bm_r + enddo + do n = 1, nbc + massc(n) = am_r*Dc(n)**bm_r + enddo + + !..Freeze water (smallest drops become cloud ice, otherwise graupel). + do k = 1, 45 + ! print*, ' Freezing water for temp = ', -k + ! ++ trude, add tadjust, so to chane temperature for where Bigg freezing starts. Follow approach in WRFV3.6 with IN + Texp = DEXP( DFLOAT(k) -t_adjust*1.0D0) - 1.0D0 ! NB Trude. Check for when texp is negative..... + ! Texp = DEXP( DFLOAT(k) ) - 1.0D0 + ! -- trude + do j = 1, ntb_r1 + do i = 1, ntb_r + lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2) + sum1 = 0.0d0 + sum2 = 0.0d0 + sumn1 = 0.0d0 + sumn2 = 0.0d0 + do n2 = nbr, 1, -1 + N_r(n2) = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2) + vol = massr(n2)*orho_w + prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp) + !++ trude + prob = MAX(prob, 0.0d0 ) + ! -- trude + if (massr(n2) .lt. xm0g) then + sumn1 = sumn1 + prob*N_r(n2) + sum1 = sum1 + prob*N_r(n2)*massr(n2) + else + sumn2 = sumn2 + prob*N_r(n2) + sum2 = sum2 + prob*N_r(n2)*massr(n2) + endif + if ((sum1+sum2) .ge. r_r(i)) EXIT + enddo + tpi_qrfz(i,j,k) = sum1 + tni_qrfz(i,j,k) = sumn1 + tpg_qrfz(i,j,k) = sum2 + tnr_qrfz(i,j,k) = sumn2 + enddo + enddo + do i = 1, ntb_c + lamc = 1.0D-6 * (Nt_c*am_r* ccg(2) * ocg1 / r_c(i))**obmr + N0_c = 1.0D-18 * Nt_c*ocg1 * lamc**cce(1) + sum1 = 0.0d0 + sumn2 = 0.0d0 + do n = nbc, 1, -1 + y = Dc(n)*1.0D6 + vol = massc(n)*orho_w + prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp) + !++ trude + prob = MAX(prob, 0.0d0 ) + ! -- trude + N_c(n) = N0_c* y**mu_c * EXP(-lamc*y)*dtc(n) + N_c(n) = 1.0D24 * N_c(n) + sumn2 = sumn2 + prob*N_c(n) + sum1 = sum1 + prob*N_c(n)*massc(n) + if (sum1 .ge. r_c(i)) EXIT + enddo + tpi_qcfz(i,k) = sum1 + tni_qcfz(i,k) = sumn2 + enddo + enddo + + IF ( this_image() == 1 ) THEN + print *, "Writing freezeH2O_mpt.dat in Thompson MP init" + OPEN(63,file="freezeH2O_mpt.dat",form="unformatted",err=9234) + WRITE(63,err=9234)tpi_qrfz + WRITE(63,err=9234)tni_qrfz + WRITE(63,err=9234)tpg_qrfz + WRITE(63,err=9234)tnr_qrfz + WRITE(63,err=9234)tpi_qcfz + WRITE(63,err=9234)tni_qcfz + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + print*, "Error writing freezeH2O_mpt.dat" + ENDIF + + endif + end subroutine freezeH2O +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Cloud ice converting to snow since portion greater than min snow +!.. size. Given cloud ice content (kg/m**3), number concentration +!.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into +!.. bins and figure out the mass/number of ice with sizes larger than +!.. D0s. Also, compute incomplete gamma function for the integration +!.. of ice depositional growth from diameter=0 to D0s. Amount of +!.. ice depositional growth is this portion of distrib while larger +!.. diameters contribute to snow growth (as in Harrington et al. 1995). +!+---+-----------------------------------------------------------------+ + + subroutine qi_aut_qs + + implicit none + +!..Local variables + INTEGER:: i, j, n2 + DOUBLE PRECISION, DIMENSION(nbi):: N_i + DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2 + REAL:: xlimit_intg + +!+---+ + + do j = 1, ntb_i1 + do i = 1, ntb_i + lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi + Di_mean = (bm_i + mu_i + 1.) / lami + N0_i = Nt_i(j)*oig1 * lami**cie(1) + t1 = 0.0d0 + t2 = 0.0d0 + if (SNGL(Di_mean) .gt. 5.*D0s) then + t1 = r_i(i) + t2 = Nt_i(j) + tpi_ide(i,j) = 0.0D0 + elseif (SNGL(Di_mean) .lt. D0i) then + t1 = 0.0D0 + t2 = 0.0D0 + tpi_ide(i,j) = 1.0D0 + else + xlimit_intg = lami*D0s + tpi_ide(i,j) = GAMMP(mu_i+2.0, xlimit_intg) * 1.0D0 + do n2 = 1, nbi + N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2) + if (Di(n2).ge.D0s) then + t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i + t2 = t2 + N_i(n2) + endif + enddo + endif + tps_iaus(i,j) = t1 + tni_iaus(i,j) = t2 + enddo + enddo + + end subroutine qi_aut_qs +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Variable collision efficiency for rain collecting cloud water using +!.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise +!.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9. +!+---+-----------------------------------------------------------------+ + + subroutine table_Efrw + + implicit none + +!..Local variables + DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw + DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X + INTEGER:: i, j + + do j = 1, nbc + do i = 1, nbr + Ef_rw = 0.0 + p = Dc(j)/Dr(i) + if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then + t_Efrw(i,j) = 0.0 + elseif (p.gt.0.25) then + X = Dc(j)*1.D6 + if (Dr(i) .lt. 75.e-6) then + Ef_rw = 0.026794*X - 0.20604 + elseif (Dr(i) .lt. 125.e-6) then + Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089 + elseif (Dr(i) .lt. 175.e-6) then + Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X & + + 0.0066237*X*X - 0.0013687*X - 0.073022 + elseif (Dr(i) .lt. 250.e-6) then + Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X & + - 0.65988 + elseif (Dr(i) .lt. 350.e-6) then + Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X & + - 0.56125 + else + Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X & + - 0.46929 + endif + else + vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) & + + 0.07934E9*Dr(i)*Dr(i)*Dr(i) & + - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i) + stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i)) + reynolds = 9.*stokes/(p*p*rho_w) + + F = DLOG(reynolds) + G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F + K0 = DEXP(G) + z = DLOG(stokes/(K0+1.D-15)) + H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z + yc0 = 2.0D0/PI2 * ATAN(H) + Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p)) + endif + + t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95)) +! ++ trude + if (Ef_rw_l) then + if (Ef_rw.ne.0.0) then + t_Efrw(i,j) = 1.0 + endif + endif +! -- trude + enddo + enddo + + end subroutine table_Efrw +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Variable collision efficiency for snow collecting cloud water using +!.. method of Wang and Ji, 2000 except equate melted snow diameter to +!.. their "effective collision cross-section." +!+---+-----------------------------------------------------------------+ + + subroutine table_Efsw + + implicit none +!..Local variables + DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw + DOUBLE PRECISION:: p, yc0, F, G, H, z, K0 + INTEGER:: i, j + + + do j = 1, nbc + vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0) + do i = 1, nbs + vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc + Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr + p = Dc(j)/Ds_m + if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 & + .or. vts.lt.1.E-3) then + t_Efsw(i,j) = 0.0 + else + stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m) + reynolds = 9.*stokes/(p*p*rho_w) + + F = DLOG(reynolds) + G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F + K0 = DEXP(G) + z = DLOG(stokes/(K0+1.D-15)) + H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z + yc0 = 2.0D0/PI2 * ATAN(H) + Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p)) + + t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95)) + +! ++ trude + if (Ef_sw_l) then + if (Ef_sw.ne.0.0) then + t_Efsw(i,j) = 1.0 + endif + endif +! -- trude + + endif + + enddo + enddo + + end subroutine table_Efsw +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Integrate rain size distribution from zero to D-star to compute the +!.. number of drops smaller than D-star that evaporate in a single +!.. timestep. Drops larger than D-star dont evaporate entirely so do +!.. not affect number concentration. +!+---+-----------------------------------------------------------------+ + + subroutine table_dropEvap + + implicit none + +!..Local variables + DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam + REAL:: xlimit_intg + INTEGER:: i, j, k + + do k = 1, ntb_r + do j = 1, ntb_r1 + lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1 + lam = lam_exp * (crg(3)*org2*org1)**obmr + N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2) + Nt_r = N0 * crg(2) / lam**cre(2) + + do i = 1, nbr + xlimit_intg = lam*Dr(i) + tnr_rev(i,j,k) = GAMMP(mu_r+1.0, xlimit_intg) * Nt_r + enddo + + enddo + enddo + + end subroutine table_dropEvap + +! TO APPLY TABLE ABOVE +!..Rain lookup table indexes. +! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) & +! * 0.78*4.*diffu(k)*xsat*rvs/rho_w) +! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) & +! / DLOG(Dr(nbr)/D0r)) +! idx_d = MAX(1, MIN(idx_d, nbr)) +! +! nir = NINT(ALOG10(rr(k))) +! do nn = nir-1, nir+1 +! n = nn +! if ( (rr(k)/10.**nn).ge.1.0 .and. & +! (rr(k)/10.**nn).lt.10.0) goto 154 +! enddo +!154 continue +! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2) +! idx_r = MAX(1, MIN(idx_r, ntb_r)) +! +! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr +! lam_exp = lamr * (crg(3)*org2*org1)**bm_r +! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1) +! nir = NINT(DLOG10(N0_exp)) +! do nn = nir-1, nir+1 +! n = nn +! if ( (N0_exp/10.**nn).ge.1.0 .and. & +! (N0_exp/10.**nn).lt.10.0) goto 155 +! enddo +!155 continue +! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3) +! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1)) +! +! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M +! * odts)) +! +!ctrlL +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ + SUBROUTINE GCF(GAMMCF,A,X,GLN) +! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS +! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS +! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY +! --- A MODIFIED LENTZ METHOD. +! --- USES GAMMLN + IMPLICIT NONE + INTEGER, PARAMETER:: ITMAX=100 + REAL, PARAMETER:: gEPS=3.E-7 + REAL, PARAMETER:: FPMIN=1.E-30 + REAL, INTENT(IN):: A, X + REAL:: GAMMCF,GLN + INTEGER:: I + REAL:: AN,B,C,D,DEL,H + GLN=GAMMLN(A) + B=X+1.-A + C=1./FPMIN + D=1./B + H=D + DO 11 I=1,ITMAX + AN=-I*(I-A) + B=B+2. + D=AN*D+B + IF(ABS(D).LT.FPMIN)D=FPMIN + C=B+AN/C + IF(ABS(C).LT.FPMIN)C=FPMIN + D=1./D + DEL=D*C + H=H*DEL + IF(ABS(DEL-1.).LT.gEPS)GOTO 1 + 11 CONTINUE + PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF' + 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H + END SUBROUTINE GCF +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + SUBROUTINE GSER(GAMSER,A,X,GLN) +! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS +! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A)) +! --- AS GLN. +! --- USES GAMMLN + IMPLICIT NONE + INTEGER, PARAMETER:: ITMAX=100 + REAL, PARAMETER:: gEPS=3.E-7 + REAL, INTENT(IN):: A, X + REAL:: GAMSER,GLN + INTEGER:: N + REAL:: AP,DEL,SUM + GLN=GAMMLN(A) + IF(X.LE.0.)THEN + IF(X.LT.0.) PRINT *, 'X < 0 IN GSER' + GAMSER=0. + RETURN + ENDIF + AP=A + SUM=1./A + DEL=SUM + DO 11 N=1,ITMAX + AP=AP+1. + DEL=DEL*X/AP + SUM=SUM+DEL + IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1 + 11 CONTINUE + PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER' + 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) + END SUBROUTINE GSER +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMLN(XX) +! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0. + IMPLICIT NONE + REAL, INTENT(IN):: XX + DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 + DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & + COF = (/76.18009172947146D0, -86.50532032941677D0, & + 24.01409824083091D0, -1.231739572450155D0, & + .1208650973866179D-2, -.5395239384953D-5/) + DOUBLE PRECISION:: SER,TMP,X,Y + INTEGER:: J + + X=XX + Y=X + TMP=X+5.5D0 + TMP=(X+0.5D0)*LOG(TMP)-TMP + SER=1.000000000190015D0 + DO 11 J=1,6 + Y=Y+1.D0 + SER=SER+COF(J)/Y +11 CONTINUE + GAMMLN=TMP+LOG(STP*SER/X) + END FUNCTION GAMMLN +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMP(A,X) +! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X) +! --- SEE ABRAMOWITZ AND STEGUN 6.5.1 +! --- USES GCF,GSER + IMPLICIT NONE + REAL, INTENT(IN):: A,X + REAL:: GAMMCF,GAMSER,GLN + GAMMP = 0. + IF((X.LT.0.) .OR. (A.LE.0.)) THEN + PRINT *, 'BAD ARGUMENTS IN GAMMP' + RETURN + ELSEIF(X.LT.A+1.)THEN + CALL GSER(GAMSER,A,X,GLN) + GAMMP=GAMSER + ELSE + CALL GCF(GAMMCF,A,X,GLN) + GAMMP=1.-GAMMCF + ENDIF + END FUNCTION GAMMP +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION WGAMMA(y) + + IMPLICIT NONE + REAL, INTENT(IN):: y + + WGAMMA = EXP(GAMMLN(y)) + + END FUNCTION WGAMMA +!+---+-----------------------------------------------------------------+ +! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS +! A FUNCTION OF TEMPERATURE AND PRESSURE +! + REAL FUNCTION RSLF(P,T) + + IMPLICIT NONE + REAL, INTENT(IN):: P, T + REAL:: ESL,X + REAL, PARAMETER:: C0= .611583699E03 + REAL, PARAMETER:: C1= .444606896E02 + REAL, PARAMETER:: C2= .143177157E01 + REAL, PARAMETER:: C3= .264224321E-1 + REAL, PARAMETER:: C4= .299291081E-3 + REAL, PARAMETER:: C5= .203154182E-5 + REAL, PARAMETER:: C6= .702620698E-8 + REAL, PARAMETER:: C7= .379534310E-11 + REAL, PARAMETER:: C8=-.321582393E-13 + + X=MAX(-80.,T-273.16) + +! ESL=612.2*EXP(17.67*X/(T-29.65)) + ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) + RSLF=.622*ESL/(P-ESL) + +! ALTERNATIVE +! ; Source: Murphy and Koop, Review of the vapour pressure of ice and +! supercooled water for atmospheric applications, Q. J. R. +! Meteorol. Soc (2005), 131, pp. 1539-1565. +! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T +! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22 +! / T - 9.44523 * ALOG(T) + 0.014025 * T)) + + END FUNCTION RSLF +!+---+-----------------------------------------------------------------+ +! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A +! FUNCTION OF TEMPERATURE AND PRESSURE +! + REAL FUNCTION RSIF(P,T) + + IMPLICIT NONE + REAL, INTENT(IN):: P, T + REAL:: ESI,X + REAL, PARAMETER:: C0= .609868993E03 + REAL, PARAMETER:: C1= .499320233E02 + REAL, PARAMETER:: C2= .184672631E01 + REAL, PARAMETER:: C3= .402737184E-1 + REAL, PARAMETER:: C4= .565392987E-3 + REAL, PARAMETER:: C5= .521693933E-5 + REAL, PARAMETER:: C6= .307839583E-7 + REAL, PARAMETER:: C7= .105785160E-9 + REAL, PARAMETER:: C8= .161444444E-12 + + X=MAX(-80.,T-273.16) + ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) + RSIF=.622*ESI/(P-ESI) + +! ALTERNATIVE +! ; Source: Murphy and Koop, Review of the vapour pressure of ice and +! supercooled water for atmospheric applications, Q. J. R. +! Meteorol. Soc (2005), 131, pp. 1539-1565. +! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T) + + END FUNCTION RSIF +!+---+-----------------------------------------------------------------+ + +!+---+-----------------------------------------------------------------+ +END MODULE module_mp_thompson +!+---+-----------------------------------------------------------------+ diff --git a/example/supporting-modules/thompson_tensors_m.f90 b/example/supporting-modules/thompson_tensors_m.f90 new file mode 100644 index 000000000..548ab5d52 --- /dev/null +++ b/example/supporting-modules/thompson_tensors_m.f90 @@ -0,0 +1,33 @@ +! Copyright (c), The Regents of the University of California +! Terms of use are as specified in LICENSE.txt +module thompson_tensors_m + !! This module supports the program in the file example/learn-microphysics-procedures.f90. + use module_mp_thompson, only : rslf, rsif + use inference_engine_m, only : tensor_t + use assert_m, only : assert + implicit none + + private + public :: T, p, y + + real, parameter :: T_min = 236.352524, T_max = 307.610779 + real, parameter :: p_min = 29671.1348, p_max = 98596.7578 + integer, parameter :: resolution = 10 + integer i + real, parameter :: T(*) = [(real(i)/real(resolution), i=0,resolution)] + real, parameter :: p(*) = [(real(i)/real(resolution), i=0,resolution)] + +contains + + elemental impure function y(x_in) result(a) + type(tensor_t), intent(in) :: x_in + type(tensor_t) a + associate(x => x_in%values()) + call assert(lbound(x,1)==1 .and. ubound(x,1)==2,"y(x) :: sufficient input") + associate(temperature => T_min + (T_max - T_min)*x(1), pressure => p_min + (p_max - p_min)*x(2) ) + a = tensor_t([rslf(pressure, temperature), rsif(pressure, temperature)]) + end associate + end associate + end function + +end module diff --git a/setup.sh b/setup.sh index 85924f0ee..7dc1503b4 100755 --- a/setup.sh +++ b/setup.sh @@ -65,7 +65,7 @@ HDF5_LIB_PATH="`brew --prefix hdf5`/lib" NETCDFF_LIB_PATH="`brew --prefix netcdf-fortran`/lib" FPM_LD_FLAG=" -L$NETCDF_LIB_PATH -L$HDF5_LIB_PATH -L$NETCDFF_LIB_PATH" -FPM_FLAG="-O3 -fallow-argument-mismatch -ffree-line-length-none -L$NETCDF_LIB_PATH -L$HDF5_LIB_PATH" +FPM_FLAG="-fcoarray=single -O3 -fallow-argument-mismatch -ffree-line-length-none -L$NETCDF_LIB_PATH -L$HDF5_LIB_PATH" FPM_FC=${FC:-"gfortran-13"} FPM_CC=${CC:-"gcc-13"}