diff --git a/README.md b/README.md index b5e292d3fc..8cd81d8f99 100644 --- a/README.md +++ b/README.md @@ -7,5 +7,3 @@ For an overview of the WRF modeling system, along with information regarding dow Information regarding WRF Model citations (including a DOI) can be found here: [https://www2.mmm.ucar.edu/wrf/users/citing_wrf.html](https://www2.mmm.ucar.edu/wrf/users/citing_wrf.html). The WRF Model is open-source code in the public domain, and its use is unrestricted. The name "WRF", however, is a registered trademark of the University Corporation for Atmospheric Research. The WRF public domain notice and related information may be found here: [https://www2.mmm.ucar.edu/wrf/users/public.html](https://www2.mmm.ucar.edu/wrf/users/public.html). - - diff --git a/Registry/registry.diags b/Registry/registry.diags index 83c63e4b4a..2d2ba5025f 100644 --- a/Registry/registry.diags +++ b/Registry/registry.diags @@ -11,7 +11,7 @@ dimspec np 2 namelist=num_press_levels z num_press_level rconfig integer p_lev_diags namelist,diags 1 0 - "flag to process vertical interp diagnostics: 0=nope, 1=yep" "flag" rconfig integer p_lev_diags_dfi namelist,diags 1 0 - "when doing p_level diags and dfi, turn off diags during 'non forecast'" rconfig integer num_press_levels namelist,diags 1 0 - "number of pressure levels to interpolate diagnostics to" "index" -rconfig real press_levels namelist,diags max_plevs 0 - "array of pressure levels to interpolate diagnostics to" "Pa" +rconfig real press_levels namelist,diags max_plevs 0 - "array of pressure levels to interpolate diagnostics to" "Pa" rconfig integer use_tot_or_hyd_p namelist,diags 1 2 - "1=use total pressure, 2=use hydrostatic pressure" "flag" rconfig integer extrap_below_grnd namelist,diags 1 1 - "1=no extrapolation, 2=extrapolate adiabatically" "flag" rconfig real p_lev_missing namelist,diags 1 -999 - "missing values below ground, no extrapolation" "constant" diff --git a/Registry/registry.em_shared_collection b/Registry/registry.em_shared_collection index 28a75f175d..e51b2a8f56 100644 --- a/Registry/registry.em_shared_collection +++ b/Registry/registry.em_shared_collection @@ -30,3 +30,6 @@ include registry.solar_fields include registry.diags include registry.iau include registry.CMAQ +#ifdef SIMRADAR +include registry.simradar +#endif diff --git a/Registry/registry.simradar b/Registry/registry.simradar new file mode 100644 index 0000000000..09ff5451d3 --- /dev/null +++ b/Registry/registry.simradar @@ -0,0 +1,255 @@ +# registry for simrad.AR +# V. Galligani, L. Fita (CIMA-IFAECI), Argentina. November 2025 +#ifdef SIMRADAR + +## L. Fita. This should work, but WRF is not prepared for 4D output variables like this + +# FROM Tiger's Registry documentation (https://www2.mmm.ucar.edu/wrf/WG2/software_2.0/registry_schaffer.pdf) +## Entry: The keyword “dimspec” +# dimspec DimName Order HowDefined CoordAxis DatName +dimspec ri 2 namelist=e_dist c e_dist +dimspec rj 2 namelist=e_angle c e_angle +dimspec rk 2 namelist=e_azimuth c e_azimuth +dimspec rn 2 namelist=nradar c nradar +dimpsec rng 2 namelist=ngrid_radar c ngrid_radar +dimspec dvh 2 namelist=Nqxbins c Nqxbins +dimspec cii 2 constant=2 z coord2 +dimspec civ 2 constant=4 z coord4 +dimspec iit 2 constant=2 z iit +## L. Fita, so, WRF output will provide a (nradar*e_radhgt) 3D variable instead +dimspec nrnn 2 namelist=nradar_ngrid c nradar_ngrid +dimspec nrnncoord 2 namelist=nradar_ngrid_coord2 c nradar_ngrid_coord2 +dimspec nrcoord 2 namelist=nradar_coord2 c nradar_coord2 +dimspec nr4coord 2 namelist=nradar_coord4 c nradar_coord4 +dimspec nncoord 2 namelist=ngrid_coord2 c ngrid_coord2 +#dimspec nrh 2 namelist=nradar_azimuth c nradar_azimuth +dimspec nrd 2 namelist=nradar_dist c nradar_dist +dimspec nbds 2 namelist=Nbands c Nbands +dimspec lqb 2 namelist=Nlutqxbins c Nlutqxbins +dimspec ddli 2 namelist=Nluticebins c Nluticebins +dimspec lnl 2 namelist=Nlutlevs c Nlutlevs +dimspec xiid - constant=12 c xiid + +# simrad.AR package namelist options +## From Tiger's Registry documentation +# rconfig [Type] [Sym] [How set] [Nentries] [Default] +rconfig integer nradar namelist,simradar 1 1 rh{11} "nradar" "amount of surface radar" "" +rconfig integer e_dist namelist,simradar 1 800 rh{11} "e_dist" "radar grid points along the radii" "" +rconfig integer e_angle namelist,simradar 1 376 rh{11} "e_angle" "radar amount of horizontal angles" "" +rconfig integer e_azimuth namelist,simradar 1 15 rh{11} "e_azimuth" "radar amount of azimuths" "" +rconfig integer ngrid_radar namelist,simradar 1 1500 rh{11} "ngrid_radar" "amount of radar cells into a given WRF grid point" "" +rconfig integer ngrid3_radar namelist,simradar 1 25000 rh{11} "ngrid3_radar" "amount of vertical radar cells into a given WRF grid point" "" +rconfig integer nradar_ngrid namelist,simradar 1 30000 rh{11} "nradar_ngrid" "nradar*ngrid_radar" "" +rconfig integer nradar_ngrid_coord2 namelist,simradar 1 60000 rh{11} "nradar_ngrid_coord2" "nradar*ngrid_radar*coord2" "" +rconfig integer nradar_coord2 namelist,simradar 1 40 rh{11} "nradar_coord2" "nradar*coord2" "" +rconfig integer nradar_coord4 namelist,simradar 1 160 rh{11} "nradar_coord4" "nradar*coord4" "" +rconfig integer ngrid_coord2 namelist,simradar 1 60000 rh{11} "ngrid_coord2" "ngrid_radar*coord2" "" +rconfig integer nradar_azimuth namelist,simradar 1 300 rh{11} "nradar_azimuth" "nradar*e_azimuth" "" +rconfig integer nradar_dist namelist,simradar 1 1600 rh{11} "nradar_dist" "nradar*e_dist" "" + +rconfig integer Nbands namelist,simradar 1 3 rh{11} "Nbands" "amount of radar bands" "" +rconfig integer Nluticebins namelist,simradar 1 200 rh{11} "Nluticebins" "amount of ice bins for LUT interpolations" "#" +rconfig integer Nlutqxbins namelist,simradar 1 50000 rh{11} "Nlutqxbins" "amount of water species bins for LUT interpolations" "#" +rconfig integer Nlutlevs namelist,simradar 1 15 rh{11} "Nlutlevs" "amount of azimuths for LUT interpolations" "" + +rconfig character radar_name namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_name" "name of the radar" "" +rconfig real radar_clon namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_clon" "longitude of the location of the radar" "" +rconfig real radar_clat namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_clat" "latitude of the location of the radar" "" +rconfig real radar_chgt namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_chgt" "height of the location of the radar" "" +rconfig character radar_band namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_band" "band of the radar: 'c', 's', or 'x'" "" +rconfig integer radar_scan namelist,simradar max_nradar (/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/) rh{11} "radar_scan" "type of scan 1:" "" +rconfig real radar_distmax namelist,simradar max_nradar (/250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250./) rh{11} "radar_distmax" "maximum distance from radar" "km" +rconfig real radar_hgtmax namelist,simradar max_nradar (/20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20./) rh{11} "radar_hgtmax" "maximum height from radar" "km" +rconfig integer em_ray namelist,simradar max_nradar (/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/) rh{11} "em_ray" "type of propagation of em ray. 1: 4/3 Ref" "" +rconfig real radards namelist,simradar max_nradar (/250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250./) rh{11} "radards" "nominal horizontal resolution of radar values in m" "" +rconfig integer radinterp namelist,simradar 1 1 rh{11} "radinterp" "method of interpolation to the radar space (1: lineal)" "" +rconfig real radii_dist namelist,simradar 1 15. rh{11} "radii_dist" "radii (km) at which arc resolution is the same as minimum ddist" "km" +rconfig integer eloption namelist,simradar 1 1 rh{11} "eloption" "Ellipsoide" "" +rconfig real missing_value namelist,simradar 1 1.e20 rh{11} "missing_value" "missing value" "" +rconfig integer doppler_scheme namelist,simradar 1 1 rh{11} "doppler_scheme" "method to compute the doppler velocities" "" + +rconfig integer timedbg namelist,simradar 1 720 rh{11} "timedbg" "frecuency of time-debugging" "minute" +rconfig real londbg namelist,simradar 1 -999. rh{11} "londbg" "longitude of grid-point debugging (-999., not used)" "degree" +rconfig real latdbg namelist,simradar 1 -999. rh{11} "latdbg" "latitude of grid-point debugging (-999., not used)" "degree" + +# LUT options (maybe from netCDF file!) +rconfig character lutfilen namelist,simradar 1 "TMATRIX_WRF_WSM6.nc" - "lutfilen" "name of the netCDF file with the T-Matrix lool-up-tables" "" +rconfig real q_rain_min namelist,simradar 1 1e-7 rh{11} "q_rain_min" "minimum rain mixing ratio to consider" "kgkg-1" +rconfig real q_rain_max namelist,simradar 1 0.02 rh{11} "q_rain_max" "maximum rain mixing ratio to consider" "kgkg-1" +rconfig real qn_rain_min namelist,simradar 1 1e-8 rh{11} "qn_rain_min" "minimum amount of nuclii for rain to consider" "#" +rconfig real qn_rain_max namelist,simradar 1 2.6e6 rh{11} "qn_rain_max" "maximum amount of nuclii for rain to consider" "#" +rconfig real q_snow_min namelist,simradar 1 1e-7 rh{11} "q_snow_min" "minimum snow mixing ratio to consider" "kgkg-1" +rconfig real q_snow_max namelist,simradar 1 0.018 rh{11} "q_snow_max" "maximum snow mixing ratio to consider" "kgkg-1" +rconfig real qn_snow_min namelist,simradar 1 1e-7 rh{11} "qn_snow_min" "minimum amount of nuclii for snow to consider" "#" +rconfig real qn_snow_max namelist,simradar 1 2.6e6 rh{11} "qn_snow_max" "maximum amount of nuclii for snow to consider" "#" +rconfig real q_grau_min namelist,simradar 1 1e-7 rh{11} "q_grau_min" "minimum grau mixing ratio to consider" "kgkg-1" +rconfig real q_grau_max namelist,simradar 1 0.018 rh{11} "q_grau_max" "maximum grau mixing ratio to consider" "kgkg-1" +rconfig real qn_grau_min namelist,simradar 1 1E-7 rh{11} "qn_grau_min" "minimum amount of nuclii for grau to consider" "#" +rconfig real qn_grau_max namelist,simradar 1 2.6e6 rh{11} "qn_grau_max" "maximum amount of nuclii for grau to consider" "#" +rconfig real q_ice_min namelist,simradar 1 1e-7 rh{11} "q_ice_min" "minimum cloud ice mixing ratio to consider" "kgkg-1" +rconfig real q_ice_max namelist,simradar 1 0.001 rh{11} "q_ice_max" "maximum cloud ice mixing ratio to consider" "kgkg-1" +rconfig real qn_ice_min namelist,simradar 1 1e-8 rh{11} "qn_ice_min" "minimum amount of nuclii for cloud ice to consider" "#" +rconfig real qn_ice_max namelist,simradar 1 2.6e6 rh{11} "qn_ice_max" "maximum amount of nuclii for cloud ice to consider" "#" +rconfig real q_hail_min namelist,simradar 1 -9999. rh{11} "q_hail_min" "minimum hail mixing ratio to consider" "kgkg-1" +rconfig real q_hail_max namelist,simradar 1 -9999. rh{11} "q_hail_max" "maximum hail mixing ratio to consider" "kgkg-1" +rconfig real qn_hail_min namelist,simradar 1 -9999. rh{11} "qn_hail_min" "minimum amount of nuclii for hail to consider" "#" +rconfig real qn_hail_max namelist,simradar 1 -9999. rh{11} "qn_hail_max" "maximum amount of nuclii for hail to consider" "#" +rconfig real q_clou_min namelist,simradar 1 1e-10 rh{11} "q_clou_min" "minimum clou mixing ratio to consider" "kgkg-1" +rconfig real q_clou_max namelist,simradar 1 0.013 rh{11} "q_clou_max" "maximum clou mixing ratio to consider" "kgkg-1" +rconfig integer dim_ice namelist,simradar 1 200 rh{11} "dim_ice" "dimension for ice bins" "-" +rconfig integer Nqxbins namelist,simradar 1 50000 rh{11} "Nqxbins" "dimension for water species bins" "-" + +## +# NOTE: WRF can not provide 5-dimensional output variables. Here amount of radars and vertical coordinate are collapsed into different dimensions, therefore outputs will be +# IJDISTGRIDRADAR(e_we,e_sn,nradar_ngrid,Time) therefore it will be needed to decompress the nradar, ngrid_radar dimensions IJDISTGRIDRADAR(e_we,e_sn,ngrid_radar,nradar,Time) +# IJGRIDRADAR(e_we,e_sn,nradar_ngrid_coord2,Time) --> IJGRIDRADAR(e_we,e_sn,ngrid_radar,coord2,nradar,Time) +# RADARLON_BNDS(e_angle,e_dist,nradar_coord4,Time) --> RADARLON_BNDS(e_angle,e_dist,coord4,nradar,Time) +# RADARLAT_BNDS(e_angle,e_dist,nradar_coord4,Time) --> RADARLAT_BNDS(e_angle,e_dist,coord4,nradar,Time) +# RADARHGT_BNDS(e_azimuth,e_dist,nradar_coord4,Time) --> RADARHGT_BNDS(e_azimuth,e_dist,coord4,nradar,Time) +# ZH_RAIN(e_angle,e_azimuth,nradar_dist,Time) --> RADARHGT_BNDS(e_angle,e_azimuth,e_dist,nradar,Time) +# +# For example: +# DO ia=1, e_angle +# DO id=1, e_dist +# DO ir=1, nradar +# DO j=1, coord4 +# irj = (ir-1)*coord4 + j +# radarlon_bnds(ia,id,irj) = lon_bnds(ia,id,j,ir) +# END DO +# END DO +# END DO +# END DO +## + +## Fromt Tiger's Registry +# state [Type] [Sym] [Dims] [Use] [NumTLev] [Stagger] [IO] [DName] [Descrip] [Units] +state integer ijraddbg {iit} misc - - rh{11} "ijraddbg" "i,j-grid points from londbg,latdbg for simrad.AR debugging" "-" +state real ellipse_a - misc - - - "ellipse_a" "semimajor axis equatorial of the ellipsoide" "m" +state real ellipse_f - misc - - - "ellipse_f" "flatteing of the ellipsoide" "-" +state character ellipse_name - misc - - - "ellipse_name" "name of the ellipsoide" "-" + +# Specific variables for the interpolation-integration of radar values retrieved from the T-MATRIX netCDF file 'lutfilen' +state real lutelevations {lnl} misc 1 - - "lutelevations" "azimuths used for T-MATRIX LUT interpolations" "degree" +state real lut_Zh_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_Zh_RAIN" "Zh rain values used for T-MATRIX LUT interpolations" "" +state real lut_Zdr_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_Zdr_RAIN" "Zdr rain values used for T-MATRIX LUT interpolations" "" +state real lut_LDR_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_LDR_RAIN" "LDR rain values used for T-MATRIX LUT interpolations" "" +state real lut_Aih_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_Aih_RAIN" "Aih rain values used for T-MATRIX LUT interpolations" "" +state real lut_Aiv_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_Aiv_RAIN" "Aiv rain values used for T-MATRIX LUT interpolations" "" +state real lut_KDP_RAIN {nbds}{lqb}{lnl} misc 1 - - "lut_KDP_RAIN" "KDP rain values used for T-MATRIX LUT interpolations" "" + +state real lut_Zh_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_Zh_SNOW" "Zh snow values used for T-MATRIX LUT interpolations" "" +state real lut_Zdr_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_Zdr_SNOW" "Zdr snow values used for T-MATRIX LUT interpolations" "" +state real lut_LDR_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_LDR_SNOW" "LDR snow values used for T-MATRIX LUT interpolations" "" +state real lut_Aih_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_Aih_SNOW" "Aih snow values used for T-MATRIX LUT interpolations" "" +state real lut_Aiv_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_Aiv_SNOW" "Aiv snow values used for T-MATRIX LUT interpolations" "" +state real lut_KDP_SNOW {nbds}{lqb}{lnl} misc 1 - - "lut_KDP_SNOW" "KDP snow values used for T-MATRIX LUT interpolations" "" + +state real lut_Zh_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_Zh_GRAU" "Zh graupel values used for T-MATRIX LUT interpolations" "" +state real lut_Zdr_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_Zdr_GRAU" "Zdr graupel values used for T-MATRIX LUT interpolations" "" +state real lut_LDR_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_LDR_GRAU" "LDR graupel values used for T-MATRIX LUT interpolations" "" +state real lut_Aih_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_Aih_GRAU" "Aih graupel values used for T-MATRIX LUT interpolations" "" +state real lut_Aiv_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_Aiv_GRAU" "Aiv graupel values used for T-MATRIX LUT interpolations" "" +state real lut_KDP_GRAU {nbds}{lqb}{lnl} misc 1 - - "lut_KDP_GRAU" "KDP graupel values used for T-MATRIX LUT interpolations" "" + +state real lut_Zh_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_Zh_ICE" "Zh cloud ice values used for T-MATRIX LUT interpolations" "" +state real lut_Zdr_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_Zdr_ICE" "Zdr cloud ice values used for T-MATRIX LUT interpolations" "" +state real lut_LDR_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_LDR_ICE" "LDR cloud ice values used for T-MATRIX LUT interpolations" "" +state real lut_Aih_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_Aih_ICE" "Aih cloud ice values used for T-MATRIX LUT interpolations" "" +state real lut_Aiv_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_Aiv_ICE" "Aiv cloud ice values used for T-MATRIX LUT interpolations" "" +state real lut_KDP_ICE {nbds}{lqb}{lnl} misc 1 - - "lut_KDP_ICE" "KDP cloud ice values used for T-MATRIX LUT interpolations" "" + +state real lut_Zh_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_Zh_HAIL" "Zh hail values used for T-MATRIX LUT interpolations" "" +state real lut_Zdr_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_Zdr_HAIL" "Zdr hail values used for T-MATRIX LUT interpolations" "" +state real lut_LDR_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_LDR_HAIL" "LDR hail values used for T-MATRIX LUT interpolations" "" +state real lut_Aih_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_Aih_HAIL" "Aih hail values used for T-MATRIX LUT interpolations" "" +state real lut_Aiv_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_Aiv_HAIL" "Aiv hail values used for T-MATRIX LUT interpolations" "" +state real lut_KDP_HAIL {nbds}{lqb}{lnl} misc 1 - - "lut_KDP_HAIL" "KDP hail values used for T-MATRIX LUT interpolations" "" + +state integer N_tocompute - misc - - - "N_tocompute" "Amount of variables to interpolate / use" "" +state integer tocompute {xiid} misc - - - "tocompute" "Indices of variables to interpolate / use" "" +state integer N_tcmoist - misc - - - "N_tcomist" "Amount of water species to interpolate / use" "" +state integer tcmoist {xiid} misc - - - "tcmoist" "Indices of water species to interpolate / use" "" + +state real q_rain_vec {lqb} misc 1 - - "q_rain_vec" "rain bin values used for T-MATRIX LUT interpolations" "kgkg-1" +state real q_snow_vec {lqb} misc 1 - - "q_snow_vec" "snow bin values used for T-MATRIX LUT interpolations" "kgkg-1" +state real q_grau_vec {lqb} misc 1 - - "q_grau_vec" "grau bin values used for T-MATRIX LUT interpolations" "kgkg-1" +state real q_ice_vec {lqb} misc 1 - - "q_ice_vec" "ice bin values used for T-MATRIX LUT interpolations" "kgkg-1" +state real q_hail_vec {lqb} misc 1 - - "q_hail_vec" "hail bin values used for T-MATRIX LUT interpolations" "kgkg-1" + +# microphysics Denisities +state real rhoo - misc 1 - - "rhoo" "microphysics basic density" "kgm-3" +state real rho_r - misc 1 - - "rho_r" "microphysics rain basic density" "kgm-3" +state real n0r - misc 1 - - "n0r" "microphysics rain basic number of particles" "#" +state real rho_s - misc 1 - - "rho_s" "microphysics snow basic density" "kgm-3" +state real n0s - misc 1 - - "n0s" "microphysics snow basic number of particles" "#" +state real rho_g - misc 1 - - "rho_g" "microphysics graupel basic density" "kgm-3" +state real n0g - misc 1 - - "n0g" "microphysics graupel basic number of particles" "#" +state real rho_h - misc 1 - - "rho_h" "microphysics hail basic density" "kgm-3" +state real n0h - misc 1 - - "n0h" "microphysics hail basic number of particles" "#" + +# All bands variables +state integer inradar ij misc 1 - rh{11} "inradar" "pertinence of the grid cell to any radar (1: yes)" "1" +state integer radarnum ij{rn} misc 1 - rh{11} "radarnum" "pertinence of the grid cell inside a given radar area" "1" +state real RADARS_DIST {ri}{rn} misc 1 - rh{11} "radars_dist" "ground radar gate distances to their center" "m" +state real RADARS_ANGLE {rj}{rn} misc 1 - rh{11} "radars_angle" "ground radar gate horizontal angles to emit beams" "m" +state real RADARS_AZIMUTH {rj}{rn} misc 1 - rh{11} "radars_azimuth" "ground radar gate vertical angles to emit beams" "m" +state real RADARLONS {rj}{ri}{rn} misc 1 - rh{11} "radarlons" "ground radar gate longitude" "degree_east" +state real RADARLATS {rj}{ri}{rn} misc 1 - rh{11} "radarlats" "ground radar gate latitude" "degree_north" +state real RADARHGTS {rk}{ri}{rn} misc 1 - rh{11} "radarhgts" "ground radar gate height" "m" +state real RADARLHGTS {rk}{ri}{rn} misc 1 - rh{11} "radarlhgts" "ground radar gate local height" "m" +state real RADARDIST {rj}{ri}{rn} misc 1 - rh{11} "radarhdist" "ground radar horizontal distance of the grid point" "m" +state real RADARLON_BNDS {rj}{ri}{nr4coord} misc 1 - rh{11} "radarlon_bnds" "ground radar gate longitude bounds" "degree_east" +state real RADARLAT_BNDS {rj}{ri}{nr4coord} misc 1 - rh{11} "radarlat_bnds" "ground radar gate latitude bounds" "degree_north" +state real RADARHGT_BNDS {rk}{ri}{nr4coord} misc 1 - rh{11} "radarhgt_bnds" "ground radar gate height bounds" "m" +state real GATEHDIST ij{rn} misc 1 - rh{11} "gatehdist" "distacnce of the grid points to the position of the radar" "m" +state integer NGRIDRADAR ij{rn} misc 1 - rh{11} "ngridradar" "amount of radar cells inside the grid point" "1" +state integer IJGRIDRADAR ij{nrnncoord} misc 1 - rh{11} "ijgridradar" "coordinates of radar cells inside the grid point" "-" +state real IJDISTGRIDRADAR ij{nrnn} misc 1 - rh{11} "ijdistgridradar" "distance of radar cells inside the grid point to its center" "m" + +# radar band values +state real ZH_RAIN {rj}{rk}{nrd} misc 1 - h{11} "ZH_RAIN" "Reflectivity for Rain" "dBZ" +state real ZDR_RAIN {rj}{rk}{nrd} misc 1 - h{11} "ZDR_RAIN" "ZDR for Rain" "dBz" +state real KDP_RAIN {rj}{rk}{nrd} misc 1 - h{11} "KDP_RAIN" "KDP for Rain" "degree/km" +state real ZH_SNOW {rj}{rk}{nrd} misc 1 - h{11} "ZH_SNOW" "Reflectivity for Snow" "dBZ" +state real ZDR_SNOW {rj}{rk}{nrd} misc 1 - h{11} "ZDR_SNOW" "ZDR for Snow" "dBz" +state real KDP_SNOW {rj}{rk}{nrd} misc 1 - h{11} "KDP_SNOW" "KDP for Snow" "degree/km" +state real ZH_GRAU {rj}{rk}{nrd} misc 1 - h{11} "ZH_GRAU" "Reflectivity for Graupel" "dBZ" +state real ZDR_GRAU {rj}{rk}{nrd} misc 1 - h{11} "ZDR_GRAU" "ZDR for Graupel" "dBz" +state real KDP_GRAU {rj}{rk}{nrd} misc 1 - h{11} "KDP_GRAU" "KDP for Graupel" "degree/km" +state real ZH_ICE {rj}{rk}{nrd} misc 1 - h{11} "ZH_ICE" "Reflectivity for cloud Ice" "dBZ" +state real ZDR_ICE {rj}{rk}{nrd} misc 1 - h{11} "ZDR_ICE" "ZDR for cloud Ice" "dBz" +state real KDP_ICE {rj}{rk}{nrd} misc 1 - h{11} "KDP_ICE" "KDP for cloud Ice" "degree/km" +state real ZH_HAIL {rj}{rk}{nrd} misc 1 - h{11} "ZH_HAIL" "Reflectivity for Hail" "dBZ" +state real ZDR_HAIL {rj}{rk}{nrd} misc 1 - h{11} "ZDR_HAIL" "ZDR for Hail" "dBz" +state real KDP_HAIL {rj}{rk}{nrd} misc 1 - h{11} "KDP_HAIL" "KDP for Hail" "degree/km" +state real ZH_TOTAL {rj}{rk}{nrd} misc 1 - h{11} "ZH_TOTAL" "Reflectivity for all species" "dBZ" +state real ZDR_TOTAL {rj}{rk}{nrd} misc 1 - h{11} "ZDR_TOTAL" "ZDR for all species" "dBz" +state real KDP_TOTAL {rj}{rk}{nrd} misc 1 - h{11} "KDP_TOTAL" "KDP for all species" "degree/km" + +state real TV_RAIN {rj}{rk}{nrd} misc 1 - h{11} "TV_RAIN" "terminal velocity for Rain" "ms-1" +state real TV_SNOW {rj}{rk}{nrd} misc 1 - h{11} "TV_SNOW" "terminal velocity for Snow" "ms-1" +state real TV_GRAU {rj}{rk}{nrd} misc 1 - h{11} "TV_GRAU" "terminal velocity for Graupel" "ms-1" +state real TV_ICE {rj}{rk}{nrd} misc 1 - h{11} "TV_ICE" "terminal velocity for cloud Ice" "ms-1" +state real TV_HAIL {rj}{rk}{nrd} misc 1 - h{11} "TV_HAIL" "terminal velocity for Hail" "ms-1" +state real TV_WGT {rj}{rk}{nrd} misc 1 - h{11} "TV_WGT" "weighted terminal velocity" "ms-1" +state real VR {rj}{rk}{nrd} misc 1 - h{11} "VR" "Doppler velocity" "ms-1" + +## OPTIONAL +## 3.dimensional variables interpolated at the radar space. +## 1. Uncomment ##SIMRADARop## +## 2. ./clean -a +## 3. Modify accordingly phys/module_simradar.F and phys/module_diagnostics_driver.F +## 4. Recompile +##SIMRADARop##state real qvradar {rj}{rk}{nrd} moist 1 - h{11} "QVRADAR" "water vapor mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qcradar {rj}{rk}{nrd} moist 1 - h{11} "QCRADAR" "cloud mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qrradar {rj}{rk}{nrd} moist 1 - h{11} "QRRADAR" "rain mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qsradar {rj}{rk}{nrd} moist 1 - h{11} "QSRADAR" "snow mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qiradar {rj}{rk}{nrd} moist 1 - h{11} "QIRADAR" "cloud ice mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qgradar {rj}{rk}{nrd} moist 1 - h{11} "QGRADAR" "graupel mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real qhradar {rj}{rk}{nrd} moist 1 - h{11} "QHRADAR" "hail mixing ratio at the radar space" "kg kg-1" +##SIMRADARop##state real uaradar {rj}{rk}{nrd} misc 1 - h{11} "UAERRADAR" "eastward Earth rotated wind speed at the radar space" "ms-1" +##SIMRADARop##state real varadar {rj}{rk}{nrd} misc 1 - h{11} "VAERRADAR" "northward Earth rotated wind speed at the radar space" "ms-1" +##SIMRADARop##state real waradar {rj}{rk}{nrd} misc 1 - h{11} "WAERRADAR" "vertical wind speed at the radar space" "ms-1" +##SIMRADARop##state real taradar {rj}{rk}{nrd} misc 1 - h{11} "TARADAR" "air temperature mixing ratio at the radar space" "K" +##SIMRADARop##state real pressradar {rj}{rk}{nrd} misc 1 - h{11} "PRESSRADAR" "air pressure at the radar space" "Pa" +#endif diff --git a/doc/README.simradAR b/doc/README.simradAR new file mode 100644 index 0000000000..5c8a0c06db --- /dev/null +++ b/doc/README.simradAR @@ -0,0 +1,124 @@ +simrad.AR: Implementation of a radar emulator in WRF + +L. Fita, V. Galligani (CIMA-IFAECI), Aregntina + +This effort is absolutely in shoulders of giants and is based in multiple different works: + + - Bessel, F. W., (1826) On the computation of geographical longitude and latitude from geodetic measurements, Astronomische Nachrichten (Astronomical Notes), Band 4, No. 86, Spalten 241-254, Altona 1826 + - Leinonen, J., High-level interface to T-matrix scattering calculations: architecture, capabilities and limitations, Opt. Express, vol. 22, issue 2, 1655-1660 (2014), doi: 10.1364/OE.22.001655. + - Mishchenko, M. I. and L. D. Travis, T-matrix computations of light scattering by large spheroidal particles, Opt. Commun., vol. 109, 16-21 (1994) + - Mishchenko, M. I., L. D. Travis, and A. Macke, Scattering of light by polydisperse, randomly oriented, finite circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996) + - Vincenty, T., (1975) Direct and inverse solutions on the ellipsoid with application of nested equations, Survey Review, Vol. 23, No. 176, pp. 88-93. + - Vincenty, T., (1976) Correspondence: solutions of geodesics, survey Review, Vol. 23, No. 180, p. 294. + - Wielaard, D. J., M. I. Mishchenko, A. Macke, and B. E. Carlson, Improved T-matrix computations for large, nonabsorbing and weakly absorbing nonspherical particles and comparison with geometrical optics approximation, Appl. Opt., vol. 36, 4305-4313 (1997) + + +There is a prototype off-line working version available from CIMA's GIT repository, from which the code is based on: +https://git.cima.fcen.uba.ar/lluis.fita/simradar/-/wikis/home + +This code is available from Lluis Fita's WRF fork (git@github.com:LluisFB/WRF.git) in branch 'simradAR' + +== Usage == +A new section into the WRF's namelist calles 'simradar' has to be set up + +=== namelist === +simradar& +/ + +integer nradar 1 amount of surface radar" "" +integer e_dist namelist,simradar 1 800 rh{11} "e_dist" "radar grid points along the radii" "" +integer e_angle namelist,simradar 1 376 rh{11} "e_angle" "radar amount of horizontal angles" "" +integer e_azimuth namelist,simradar 1 15 rh{11} "e_azimuth" "radar amount of azimuths" "" +integer ngrid_radar namelist,simradar 1 1500 rh{11} "ngrid_radar" "amount of radar cells into a given WRF grid point" "" +integer ngrid3_radar namelist,simradar 1 25000 rh{11} "ngrid3_radar" "amount of vertical radar cells into a given WRF grid point" "" +integer nradar_ngrid namelist,simradar 1 30000 rh{11} "nradar_ngrid" "nradar*ngrid_radar" "" +integer nradar_ngrid_coord2 namelist,simradar 1 60000 rh{11} "nradar_ngrid_coord2" "nradar*ngrid_radar*coord2" "" +integer nradar_coord2 namelist,simradar 1 40 rh{11} "nradar_coord2" "nradar*coord2" "" +integer nradar_coord4 namelist,simradar 1 160 rh{11} "nradar_coord4" "nradar*coord4" "" +integer ngrid_coord2 namelist,simradar 1 60000 rh{11} "ngrid_coord2" "ngrid_radar*coord2" "" +integer nradar_azimuth namelist,simradar 1 300 rh{11} "nradar_azimuth" "nradar*e_azimuth" "" +integer nradar_dist namelist,simradar 1 1600 rh{11} "nradar_dist" "nradar*e_dist" "" + +integer Nbands namelist,simradar 1 3 rh{11} "Nbands" "amount of radar bands" "" +integer Nluticebins namelist,simradar 1 200 rh{11} "Nluticebins" "amount of ice bins for LUT interpolations" "#" +integer Nlutqxbins namelist,simradar 1 50000 rh{11} "Nlutqxbins" "amount of water species bins for LUT interpolations" "#" +integer Nlutlevs namelist,simradar 1 15 rh{11} "Nlutlevs" "amount of azimuths for LUT interpolations" "" + +character radar_name namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_name" "name of the radar" "" +real radar_clon namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_clon" "longitude of the location of the radar" "" +real radar_clat namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_clat" "latitude of the location of the radar" "" +real radar_chgt namelist,simradar max_nradar (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) rh{11} "radar_chgt" "height of the location of the radar" "" +character radar_band namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_band" "band of the radar: 'c', 's', or 'x'" "" +integer radar_scan namelist,simradar max_nradar (/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/) rh{11} "radar_scan" "type of scan 1:" "" +real radar_distmax namelist,simradar max_nradar (/250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250./) rh{11} "radar_distmax" "maximum distance from radar" "km" +real radar_hgtmax namelist,simradar max_nradar (/20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20.,20./) rh{11} "radar_hgtmax" "maximum height from radar" "km" +integer em_ray namelist,simradar max_nradar (/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/) rh{11} "em_ray" "type of propagation of em ray. 1: 4/3 Ref" "" +real radards namelist,simradar max_nradar (/250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250.,250./) rh{11} "radards" "nominal horizontal resolution of radar values in m" "" +integer radinterp namelist,simradar 1 1 rh{11} "radinterp" "method of interpolation to the radar space (1: lineal)" "" +real radii_dist namelist,simradar 1 15. rh{11} "radii_dist" "radii (km) at which arc resolution is the same as minimum ddist" "km" +integer eloption namelist,simradar 1 1 rh{11} "eloption" "Ellipsoide" "" +real missing_value namelist,simradar 1 1.e20 rh{11} "missing_value" "missing value" "" +integer doppler_scheme namelist,simradar 1 1 rh{11} "doppler_scheme" "method to compute the doppler velocities" "" + +integer timedbg namelist,simradar 1 720 rh{11} "timedbg" "frecuency of time-debugging" "minute" +real londbg namelist,simradar 1 -999. rh{11} "londbg" "longitude of grid-point debugging (-999., not used)" "degree" +real latdbg namelist,simradar 1 -999. rh{11} "latdbg" "latitude of grid-point debugging (-999., not used)" "degree" + +# LUT options (maybe from netCDF file!) +character lutfilen namelist,simradar 1 "TMATRIX_WRF_WSM6.nc" - "lutfilen" "name of the netCDF file with the T-Matrix lool-up-tables" "" +real q_rain_min namelist,simradar 1 1e-7 rh{11} "q_rain_min" "minimum rain mixing ratio to consider" "kgkg-1" +real q_rain_max namelist,simradar 1 0.02 rh{11} "q_rain_max" "maximum rain mixing ratio to consider" "kgkg-1" +real qn_rain_min namelist,simradar 1 1e-8 rh{11} "qn_rain_min" "minimum amount of nuclii for rain to consider" "#" +real qn_rain_max namelist,simradar 1 2.6e6 rh{11} "qn_rain_max" "maximum amount of nuclii for rain to consider" "#" +real q_snow_min namelist,simradar 1 1e-7 rh{11} "q_snow_min" "minimum snow mixing ratio to consider" "kgkg-1" +real q_snow_max namelist,simradar 1 0.018 rh{11} "q_snow_max" "maximum snow mixing ratio to consider" "kgkg-1" +real qn_snow_min namelist,simradar 1 1e-7 rh{11} "qn_snow_min" "minimum amount of nuclii for snow to consider" "#" +real qn_snow_max namelist,simradar 1 2.6e6 rh{11} "qn_snow_max" "maximum amount of nuclii for snow to consider" "#" +real q_grau_min namelist,simradar 1 1e-7 rh{11} "q_grau_min" "minimum grau mixing ratio to consider" "kgkg-1" +real q_grau_max namelist,simradar 1 0.018 rh{11} "q_grau_max" "maximum grau mixing ratio to consider" "kgkg-1" +real qn_grau_min namelist,simradar 1 1E-7 rh{11} "qn_grau_min" "minimum amount of nuclii for grau to consider" "#" +real qn_grau_max namelist,simradar 1 2.6e6 rh{11} "qn_grau_max" "maximum amount of nuclii for grau to consider" "#" +real q_ice_min namelist,simradar 1 1e-7 rh{11} "q_ice_min" "minimum cloud ice mixing ratio to consider" "kgkg-1" +real q_ice_max namelist,simradar 1 0.001 rh{11} "q_ice_max" "maximum cloud ice mixing ratio to consider" "kgkg-1" +real qn_ice_min namelist,simradar 1 1e-8 rh{11} "qn_ice_min" "minimum amount of nuclii for cloud ice to consider" "#" +real qn_ice_max namelist,simradar 1 2.6e6 rh{11} "qn_ice_max" "maximum amount of nuclii for cloud ice to consider" "#" +real q_hail_min namelist,simradar 1 -9999. rh{11} "q_hail_min" "minimum hail mixing ratio to consider" "kgkg-1" +real q_hail_max namelist,simradar 1 -9999. rh{11} "q_hail_max" "maximum hail mixing ratio to consider" "kgkg-1" +real qn_hail_min namelist,simradar 1 -9999. rh{11} "qn_hail_min" "minimum amount of nuclii for hail to consider" "#" +real qn_hail_max namelist,simradar 1 -9999. rh{11} "qn_hail_max" "maximum amount of nuclii for hail to consider" "#" +real q_clou_min namelist,simradar 1 1e-10 rh{11} "q_clou_min" "minimum clou mixing ratio to consider" "kgkg-1" +real q_clou_max namelist,simradar 1 0.013 rh{11} "q_clou_max" "maximum clou mixing ratio to consider" "kgkg-1" +integer dim_ice namelist,simradar 1 200 rh{11} "dim_ice" "dimension for ice bins" "-" +integer Nqxbins namelist,simradar 1 50000 rh{11} "Nqxbins" "dimension for water species bins" "-" + + +== Installation == + +1.- Get the repository from a $WORKDIR +git clone git@github.com:LluisFB/WRF.git + +2.- getting submoules +git submodule update --init --recursive + +3.- Switching to the simrad.AR branch +git checkout simradAR + +4.- Preparing installation (remember to define WRFIO_NCD_LARGE_FILE_SUPPORT=1 and NETCDF) +./configure + +5.- Editing configure.wrf to add -DSIMRADAR pre-compilation flag +cp configure.wrf configure.simradAR.orig +cp configure.wrf configure.simradAR.wrf +vim configure.simradAR.wrf +diff configure.wrf.orig configure.simradAR.wrf +216c216 +< \ +--- +> -DSIMRADAR \ +cp configure.simradAR.wrf configure.wrf + +6.- Compile as usual +./compile em_real >& compile.log + + + diff --git a/dyn_em/start_em.F b/dyn_em/start_em.F index 96057ebdfb..8950b8602f 100644 --- a/dyn_em/start_em.F +++ b/dyn_em/start_em.F @@ -52,6 +52,11 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & USE module_model_constants USE module_avgflx_em, ONLY : zero_avgflx +#ifdef SIMRADAR +! simrad.AR + USE module_simradar_tools, ONLY: lut_vars +#endif + IMPLICIT NONE ! Input data. TYPE (domain) :: grid @@ -135,6 +140,9 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & REAL :: dt_s, dx_km REAL :: max_mf, max_rot_angle +#ifdef SIMRADAR + TYPE(lut_vars) :: lutgrid +#endif CALL get_ijk_from_grid ( grid , & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & @@ -1271,6 +1279,24 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & ,grid%NC_CU, grid%NI_CU, grid%NR_CU, grid%NS_CU,grid%CCN_CU & ,grid%alevsiz_cu,grid%num_months,grid%no_src_types_cu,grid%aeromcu & ,grid%aeropcu,grid%EFCG,grid%EFCS,grid%EFIG,grid%EFIS,grid%EFSG,grid%EFSS & +#ifdef SIMRADAR +! simrad.AR + ,grid%xlong_v, grid%xlat_u & + ,config_flags%nradar & + ,config_flags%e_angle, config_flags%e_azimuth, config_flags%e_dist & + ,config_flags%nradar_coord4, config_flags%nradar_ngrid_coord2 & + ,config_flags%nradar_ngrid, model_config_rec%ngrid_coord2 & + ,config_flags%londbg, config_flags%latdbg & + ,grid%ijraddbg & + ,grid%radars_dist, grid%radars_angle, grid%radars_azimuth & + ,grid%ellipse_a, grid%ellipse_f, grid%ellipse_name & + ,grid%gatehdist, grid%inradar, grid%radarnum & + ,grid%radarlons, grid%radarlats, grid%radardist & + ,grid%radarhgts, grid%radarlhgts & + ,grid%radarlon_bnds, grid%radarlat_bnds, grid%radarhgt_bnds & + ,grid%ngridradar, grid%ijgridradar, grid%ijdistgridradar & + ,lutgrid & +#endif ) ENDDO ! loop of tiles for phy_init @@ -1297,6 +1323,68 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & CALL wrf_debug ( 100 , 'start_domain_em: After call to phy_init' ) +#ifdef SIRADAR + !! Filling grid with lutgrid values + grid%lutelevations = grid%elevations + grid%lut_Zh_RAIN = lutgrid%Zh_RAIN + grid%lut_Zdr_RAIN = lutgrid%Zdr_RAIN + grid%lut_LDR_RAIN = lutgrid%LDR_RAIN + grid%lut_Aih_RAIN = lutgrid%Aih_RAIN + grid%lut_Aiv_RAIN = lutgrid%Aiv_RAIN + grid%lut_KDP_RAIN = lutgrid%KDP_RAIN + + grid%lut_Zh_SNOW = lutgrid%Zh_SNOW + grid%lut_Zdr_SNOW = lutgrid%Zdr_SNOW + grid%lut_LDR_SNOW = lutgrid%LDR_SNOW + grid%lut_Aih_SNOW = lutgrid%Aih_SNOW + grid%lut_Aiv_SNOW = lutgrid%Aiv_SNOW + grid%lut_KDP_SNOW = lutgrid%KDP_SNOW + + grid%lut_Zh_GRAU = lutgrid%Zh_GRAU + grid%lut_Zdr_GRAU = lutgrid%Zdr_GRAU + grid%lut_LDR_GRAU = lutgrid%LDR_GRAU + grid%lut_Aih_GRAU = lutgrid%Aih_GRAU + grid%lut_Aiv_GRAU = lutgrid%Aiv_GRAU + grid%lut_KDP_GRAU = lutgrid%KDP_GRAU + + grid%lut_Zh_ICE = lutgrid%Zh_ICE + grid%lut_Zdr_ICE = lutgrid%Zdr_ICE + grid%lut_LDR_ICE = lutgrid%LDR_ICE + grid%lut_Aih_ICE = lutgrid%Aih_ICE + grid%lut_Aiv_ICE = lutgrid%Aiv_ICE + grid%lut_KDP_ICE = lutgrid%KDP_ICE + + grid%lut_Zh_HAIL = lutgrid%Zh_HAIL + grid%lut_Zdr_HAIL = lutgrid%Zdr_HAIL + grid%lut_LDR_HAIL = lutgrid%LDR_HAIL + grid%lut_Aih_HAIL = lutgrid%Aih_HAIL + grid%lut_Aiv_HAIL = lutgrid%Aiv_HAIL + grid%lut_KDP_HAIL = lutgrid%KDP_HAIL + + grid%N_tocompute = lutgrid%N_tocompute + grid%tocompute = lutgrid%tocompute + grid%N_tcmoist = lutgrid%Nmoist + grid%tcmoist = lutgrid%imoist + + grid%q_rain_vec = lutgrid%q_rain_vec + grid%q_snow_vec = lutgrid%q_snow_vec + grid%q_grau_vec = lutgrid%q_grau_vec + grid%q_ice_vec = lutgrid%q_ice_vec + grid%q_hail_vec = lutgrid%q_rhail_vec + + ! microphysics-derived Denisities + grid%rhoo = lutgrid%rhoo + grid%rho_r = lutgrid%rho_r + grid%n0r = lutgrid%n0r + grid%rho_s = lutgrid%rho_s + grid%n0s = lutgrid%n0s + grid%rho_g = lutgrid%rho_g + grid%n0g = lutgrid%n0g + grid%rho_h = lutgrid%rho_h + grid%n0h = lutgrid%n0h + +#endif + IF (config_flags%do_avgflx_em .EQ. 1) THEN WRITE ( message , FMT = '("start_em: initializing avgflx on domain ",I3)' ) & & grid%id diff --git a/frame/module_driver_constants.F b/frame/module_driver_constants.F index 6c7d797d1a..1247ffa9cf 100644 --- a/frame/module_driver_constants.F +++ b/frame/module_driver_constants.F @@ -102,6 +102,12 @@ MODULE module_driver_constants ! The maximum number of domains used by the external model with which wrf is communicating in coupled mode INTEGER , PARAMETER :: max_extdomains = 5 + +#ifdef SIMRADAR +! L. Fita simrad.AR + ! The maximum amunt of surface radars + INTEGER , PARAMETER :: max_nradar = 20 +#endif ! 2. Following related to driver level data structures for DM_PARALLEL communications diff --git a/phys/Makefile b/phys/Makefile index 5eda61c111..2c185dd02e 100644 --- a/phys/Makefile +++ b/phys/Makefile @@ -209,7 +209,9 @@ MODULES = \ module_wind_mav.o \ module_sf_lake.o \ module_diagnostics_driver.o \ - module_irrigation.o + module_irrigation.o \ + module_simradar_tools.o \ + module_simradar.o FIRE_MODULES = \ module_fr_fire_driver.o \ diff --git a/phys/module_diagnostics_driver.F b/phys/module_diagnostics_driver.F index aa583b505f..b61ef0c9be 100644 --- a/phys/module_diagnostics_driver.F +++ b/phys/module_diagnostics_driver.F @@ -61,7 +61,14 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & model_config_rec USE module_streams +#ifdef SIMRADAR +! simrad.AR CIMA Feb. 26 + USE module_utility + USE module_simradar, ONLY: simradar_calc + USE module_simradar_tools, ONLY: lut_vars +#else USE module_utility, ONLY : WRFU_Time +#endif !============================================================= ! USE Association for the Diagnostic Packages @@ -81,7 +88,6 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & IMPLICIT NONE - !============================================================= ! Subroutine Arguments !============================================================= @@ -171,6 +177,16 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & TYPE(WRFU_Time) :: currentTime +#ifdef SIMRADAR + TYPE(lut_vars) :: lutgrid + LOGICAL :: is_output_t_simradar + TYPE(WRFU_Time) :: hist_time, aux11_time, CurrTime, StartTime + TYPE(WRFU_TimeInterval) :: dtint, histint, aux11int + CHARACTER(len=256) :: message, timestr + INTEGER :: i + CHARACTER(len=1024) :: msg +#endif + !============================================================= ! Start of executable code !============================================================= @@ -1212,6 +1228,142 @@ SUBROUTINE diagnostics_driver ( grid, config_flags, & kts=k_start, kte=min(k_end,kde-1), num_tiles=grid%num_tiles ) END IF SOLAR_FIELDS +#ifdef SIMRADAR + IF ( ( config_flags%auxhist11_interval == 0 ) ) THEN + WRITE (diag_message , * ) & + "SIMRAD.AR: ERROR -- error -- ERROR -- error : NO 'auxhist11_interval' has been " // & + "defined in 'namelist.input'" + CALL wrf_error_fatal ( diag_message ) + END IF + + CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: SIMRAD.AR DIAGNOSTICS' ) + ! Do not forget, SIMRADAR initialization in phys/module_physics_init.F at SUBROUTINE simradar_init + + ! From AFWA diagnostics ... + CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, & + ringinterval=histint) + CALL WRFU_ALARMGET( grid%alarms(AUXHIST11_ALARM), prevringtime=aux11_time, ringinterval=aux11int) + CALL domain_clock_get ( grid, current_time=CurrTime, simulationStartTime=StartTime, & + current_timestr=timestr, time_step=dtint ) + !is_output_t_simradar = ( Currtime .GE. hist_time + histint - dtint .OR. & + ! Currtime .ge. aux11_time + aux19int - dtint ) + is_output_t_simradar = ( Currtime .ge. aux11_time + aux11int - dtint ) + ! Complementary way... + !is_output_t_simradar=(MOD(NINT((curr_secs2+grid%dt)/grid%dt),NINT(config_flags%auxhist11_interval*60./grid%dt))==0 .OR. curr_secs2==0.) + + ! microphysics-derived Denisities. Scalars, can be duplicated + lutgrid%rhoo = grid%rhoo + lutgrid%rho_r = grid%rho_r + lutgrid%n0r = grid%n0r + lutgrid%rho_s = grid%rho_s + lutgrid%n0s = grid%n0s + lutgrid%rho_g= grid%rho_g + lutgrid%n0g = grid%n0g + lutgrid%rho_h = grid%rho_h + lutgrid%n0h = grid%n0h + + IF (is_output_t_simradar) THEN + CALL simradar_calc( & + IS_RESTART=config_flags%restart, IS_OUTPUT_TIME=is_output_t_simradar, & + NSTEPS=grid%nsteps, & + DT=grid%dt, XTIME=grid%xtime, CURR_SECS2=curr_secs2, & + + TIMEDBG=config_flags%timedbg, & + LONDBG=config_flags%londbg, LATDBG=config_flags%latdbg, & + IDBG=grid%ijraddbg(1), JDBG=grid%ijraddbg(2), & + + ! Input + N_MOIST=num_moist, & + MOIST=grid%moist, & + PRM_FIRST_SCALAR=param_first_scalar, & + M_QV=P_QV, M_QC=P_QC, M_QR=P_QR, M_QS=P_QS, M_QI=P_QI, M_QG=P_QG, M_QH=P_QH, & + M_QNG=P_QNG, M_QNH=P_QNH, M_QNR=P_QNR, M_QVOLG=P_QVOLG, & +! M_QG=F_QV, F_QC, F_QI, F_QS, & +! QV=grid%moist(:,:,:,m_qv), QC=grid%moist(:,:,:,m_qc), & +! QR=grid%moist(:,:,:,m_qr), QS=grid%moist(:,:,:,m_qs), & +! QI=grid%moist(:,:,:,m_qi), QG=grid%moist(:,:,:,m_qg), & +! QH=grid%moist(:,:,:,m_qh), & + LON=grid%xlong, LAT=grid%xlat, & + SINA=grid%sina, COSA=grid%cosa, & + U=grid%u_phy, V=grid%v_phy, W=grid%w_2, & + P=grid%p, PB=grid%pb, & + TA=grid%t_phy, PH=grid%ph_2, PHB=grid%phb, & + NRADAR=config_flags%nradar, NGRID_RADAR=config_flags%ngrid_radar, & + E_ANGLE=config_flags%e_angle, E_AZIMUTH=config_flags%e_azimuth, & + E_DIST=config_flags%e_dist, & + COORD2=2, COORD3=3, COORD4=4, AROUND2=9, AROUND3=27, & + NGRID3_RADAR=config_flags%ngrid3_radar, NRD=config_flags%nradar_dist, & + NRNN=config_flags%nradar_ngrid, NRNNCOORD=config_flags%nradar_ngrid_coord2, & + NRCOORD=config_flags%nradar_coord2, NR4COORD=config_flags%nradar_coord4, & + ELIPSN=grid%ellipse_name, A=grid%ellipse_a, F=grid%ellipse_f, & + RADARS_ANGLE=grid%radars_angle, RADARS_AZIMUTH=grid%radars_azimuth, & + RADARS_DIST=grid%radars_dist, & + IJRADARDBG=grid%ijraddbg, & + RADARLONS=grid%radarlons, RADARLATS=grid%radarlats, & + RADARHGTS=grid%radarhgts, RADARLHGTS=grid%radarlhgts, & + GATEHDIST=grid%gatehdist, INRADAR=grid%inradar, RADNUM=grid%radarnum, & + NGRIDRADAR=grid%ngridradar, IJGRIDRADAR=grid%ijgridradar, & + IJDISTGRIDRADAR=grid%ijdistgridradar, & + NBANDS=config_flags%Nbands, NLUTQXBINS=config_flags%Nlutqxbins, & + NLUTICEBINS=config_flags%Nluticebins, NLUTLEVS=config_flags%Nlutlevs, & + ELEVATIONS=grid%lutelevations, & + lutZH_RAIN=grid%lut_Zh_RAIN, lutZDR_RAIN=grid%lut_Zdr_RAIN, & + ! Not needed yet... +! lutLDR_RAIN=grid%lut_LDR_RAIN, lutAIH_RAIN=grid%lut_Aih_RAIN, & +! lutAIV_RAIN=grid%lut_Aiv_RAIN, lutKDP_RAIN=grid%lut_KDP_RAIN, & + lutKDP_RAIN=grid%lut_KDP_RAIN, & + lutZH_SNOW=grid%lut_Zh_SNOW, lutZDR_SNOW=grid%lut_Zdr_SNOW, & +! lutLDR_SNOW=grid%lut_LDR_SNOW, lutAIH_SNOW=grid%lut_Aih_SNOW, & +! lutAIV_SNOW=grid%lut_Aiv_SNOW, lutKDP_SNOW=grid%lut_KDP_SNOW, & + lutKDP_SNOW=grid%lut_KDP_SNOW, & + lutZH_GRAU=grid%lut_Zh_GRAU, lutZDR_GRAU=grid%lut_Zdr_GRAU, & +! lutLDR_GRAU=grid%lut_LDR_GRAU, lutAIH_GRAU=grid%lut_Aih_GRAU, & +! lutAIV_GRAU=grid%lut_Aiv_GRAU, lutKDP_GRAU=grid%lut_KDP_GRAU, & + lutKDP_GRAU=grid%lut_KDP_GRAU, & + lutZH_ICE=grid%lut_Zh_ICE, lutZDR_ICE=grid%lut_Zdr_ICE, & +! lutLDR_ICE=grid%lut_LDR_ICE, lutAIH_ICE=grid%lut_Aih_ICE, & +! lutAIV_ICE=grid%lut_Aiv_ICE, lutKDP_ICE=grid%lut_KDP_ICE, & + lutKDP_ICE=grid%lut_KDP_ICE, & + lutZH_HAIL=grid%lut_Zh_HAIL, lutZDR_HAIL=grid%lut_Zdr_HAIL, & +! lutLDR_HAIL=grid%lut_LDR_HAIL, lutAIH_HAIL=grid%lut_Aih_HAIL, & +! lutAIV_HAIL=grid%lut_Aiv_HAIL, lutKDP_HAIL=grid%lut_KDP_HAIL, & + lutKDP_HAIL=grid%lut_KDP_HAIL, & + N_TOCOMPUTE=grid%N_tocompute, TOCOMPUTE=grid%tocompute, & + NMOIST=grid%N_tcmoist, IMOIST=grid%tcmoist, & + ! These are too big to be duplicated + NQXBINS = config_flags%Nqxbins, & + Q_RAIN_VEC=grid%q_rain_vec, Q_SNOW_VEC=grid%q_snow_vec, Q_GRAU_VEC=grid%q_grau_vec, & + Q_ICE_VEC=grid%q_ice_vec, Q_HAIL_VEC=grid%q_hail_vec, & + LUTGRID=lutgrid, & + + ! Output + !!! Optional. Modify Registry/registry.simradar and recompile + !!SIMRADARop!!IQVM=grid%qvradar, IQCM=grid%qcradar, IQRM=grid%qrradar, & + !!SIMRADARop!!IQIM=grid%qiradar, IQSM=grid%qvsadar, IQGM=grid%qgradar, & + !!SIMRADARop!!IQHM=grid%qhradar, & + !!SIMRADARop!!IUAM=grid%uaradar, IVAM=grid%varadar, IWAM=grid%waradar, & + !!SIMRADARop!!ITAM=grid%taradar, IPRESSM=grid%pressradar, & + ZH_RAIN=grid%zh_rain, ZDR_RAIN=grid%zdr_rain, KDP_RAIN=grid%kdp_rain, & + ZH_SNOW=grid%zh_snow, ZDR_SNOW=grid%zdr_snow, KDP_SNOW=grid%kdp_snow, & + ZH_GRAU=grid%zh_grau, ZDR_GRAU=grid%zdr_grau, KDP_GRAU=grid%kdp_grau, & + ZH_ICE=grid%zh_ice, ZDR_ICE=grid%zdr_ice, KDP_ICE=grid%kdp_ice, & + ZH_HAIL=grid%zh_hail, ZDR_HAIL=grid%zdr_hail, KDP_HAIL=grid%kdp_hail, & + ZH_TOTAL=grid%zh_total, ZDR_TOTAL=grid%zdr_total, KDP_TOTAL=grid%kdp_total, & + TV_RAIN=grid%tv_rain, TV_SNOW=grid%tv_snow, TV_GRAU=grid%tv_grau, & + TV_ICE=grid%tv_ice, TV_HAIL=grid%tv_hail, & + TV_WGT=grid%tv_wgt, VR=grid%vr, & + + ! Dimension arguments + IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde, & + IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme, & + IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe, & + I_START=grid%i_start,I_END=min(grid%i_end, ide-1), & + J_START=grid%j_start,J_END=min(grid%j_end, jde-1), & + KTS=k_start, KTE=k_end, & + NUM_TILES=grid%num_tiles & + ) + END IF +#endif END SUBROUTINE diagnostics_driver diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 8f94e171e2..1169d16643 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -12,11 +12,16 @@ MODULE module_physics_init USE module_state_description USE module_model_constants USE module_configure, ONLY : grid_config_rec_type + #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) USE module_dm, ONLY : wrf_dm_max_real #endif USE module_ra_clWRF_support +#ifdef SIMRADAR + USE module_simradar_tools, ONLY: lut_vars +#endif + ! USE module_ssib_veg , ONLY : init_module_ssib_veg !fds (SSiB constants) !Local data for CAM's MG MP scheme integer :: ixcldliq, ixcldice, ixnumliq, ixnumice @@ -273,7 +278,25 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & ,alevsiz_cu,num_months,no_src_types_cu,aeromcu,aeropcu & ! PSH/TWG 06/10/16 ,EFCG,EFCS,EFIG,EFIS,EFSG,EFSS & ! TWG #endif - ) + +#ifdef SIMRADAR +! simrad.AR + ,xlong_v, xlat_u & + ,nradar & + ,e_angle, e_azimuth, e_dist & + ,nr4coord, nrnncoord & + ,nrnn, nncoord & + ,londbg, latdbg & + ,ijdbg & + ,radars_dist, radars_angle, radars_azimuth & + ,ra, rf, elipsn & + ,gatehdist, inradar, radnum & + ,radarlons, radarlats, radardist, radarhgts, radarlhgts & + ,radarlon_bnds, radarlat_bnds, radarhgt_bnds & + ,ngridradar, ijgridradar, ijdistgridradar & + ,lutgrid & +#endif + ) !----------------------------------------------------------------- USE module_domain @@ -910,6 +933,38 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & CHARACTER(LEN=8) :: name !----------------------------------------------------------------- +#ifdef SIMRADAR +! simrad.AR + INTEGER, INTENT(in) :: nradar, e_angle, e_dist, e_azimuth, & + nr4coord, nrnncoord, nrnn, nncoord + REAL, INTENT(in) :: londbg, latdbg + REAL, DIMENSION(ims:ime, jms:jme), INTENT(in) :: xlong_v, xlat_u + INTEGER, DIMENSION(2), INTENT(out) :: ijdbg + REAL, INTENT(out) :: ra, rf + CHARACTER(len=256),INTENT(out) :: elipsn + REAL, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(out) :: gatehdist + INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(out) :: inradar + INTEGER, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(out) :: radnum + REAL, DIMENSION(e_dist, nradar),INTENT(out) :: radars_dist + REAL, DIMENSION(e_angle, nradar),INTENT(out) :: radars_angle + REAL, DIMENSION(e_azimuth, nradar),INTENT(out) :: radars_azimuth + REAL, DIMENSION(e_angle, e_dist, nradar),INTENT(out) :: radarlons, radarlats, radardist + REAL, DIMENSION(e_azimuth, e_dist, nradar),INTENT(out) :: radarhgts, radarlhgts + REAL,DIMENSION(e_angle,e_dist,nr4coord), INTENT(out) :: radarlon_bnds, radarlat_bnds + REAL, DIMENSION(e_azimuth,e_dist,nr4coord), INTENT(out) & + :: radarhgt_bnds + INTEGER, DIMENSION(ims:ime,jms:jme, nradar), & + INTENT(out) :: ngridradar + INTEGER, DIMENSION(ims:ime, jms:jme, nrnncoord), & + INTENT(out) :: ijgridradar + REAL, DIMENSION(ims:ime, jms:jme, nrnn), INTENT(out) :: ijdistgridradar + + ! LUT space + TYPE(lut_vars), INTENT(out) :: lutgrid +#endif + #if ( EM_CORE == 1 ) ! Compute 2d grid distance and 2d grid cell area. For use with @@ -1739,6 +1794,37 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & end if #endif +#ifdef SIMRADAR +! CALL simradar_init(model_config_rec & +! ,xlong_u, xlat_u & + CALL simradar_init(xlong, xlat & + ,xlong_u, xlat_u & + ,xlong_v, xlat_v & + ,londbg, latdbg & + ,ids, ide, jds, jde, kds, kde & + ,ims, ime, jms, jme, kms, kme & + ,its, ite, jts, jte, kts, kte & + ,param_first_scalar & + ,p_qv, p_qc, p_qr, p_qs, p_qi, p_qg, p_qh & + + ,model_config_rec%nradar, model_config_rec%ngrid_radar & + ,model_config_rec%nradar_ngrid, model_config_rec%nradar_ngrid_coord2 & + ,model_config_rec%nradar_coord2, model_config_rec%nradar_coord4 & + ,model_config_rec%ngrid_coord2 & + ,model_config_rec%e_angle, model_config_rec%e_azimuth, model_config_rec%e_dist & + + ! output + ,ijdbg & + ,radars_dist, radars_angle, radars_azimuth & + ,ra, rf, elipsn & + ,gatehdist, inradar, radnum & + ,radarlons, radarlats, radardist, radarhgts, radarlhgts & + ,radarlon_bnds, radarlat_bnds, radarhgt_bnds & + ,ngridradar, ijgridradar, ijdistgridradar & + ,lutgrid & + ) +#endif + END SUBROUTINE phy_init !===================================================================== @@ -5747,4 +5833,620 @@ SUBROUTINE shalwater_init(ims,ime,jms,jme, & END SUBROUTINE shalwater_init +#ifdef SIMRADAR + SUBROUTINE simradar_init(xlong, xlat & + ,xlong_u, xlat_u & + ,xlong_v, xlat_v & + ,londbg, latdbg & + ,ids, ide, jds, jde, kds, kde & + ,ims, ime, jms, jme, kms, kme & + ,its, ite, jts, jte, kts, kte & + ,param_first_scalar & + ,p_qv, p_qc, p_qr, p_qs, p_qi, p_qg, p_qh & + ,nradar, ngrid_radar & + ,nrnn, nrnncoord & + ,nrcoord, nr4coord & + ,nncoord & + ,e_angle, e_azimuth, e_dist & +! ,gatelon, gatelat, gatehgt, distmax, heightmax, radards & + + ! output + ,ijdbg & + ,radars_dist, radars_angle, radars_azimuth & + ,ra, rf, elipsn & + ,gatehdist, inradar, radnum & + ,radarlons, radarlats, radardist, radarhgts, radarlhgts & + ,radarlon_bnds, radarlat_bnds, radarhgt_bnds & + ,ngridradar, ijgridradar, ijdistgridradar & + ,lutgrid & + ) + + ! Subroutine to initialize the simrad.AR syntehtic surface radar + ! NOTE radars are in polar space: + ! dx: distance dy: angle dz: azimuth + ! 3Dvariables (dy,dz,dx), lon,lats(dy,dx), hgt(dx,dz) + + ! There is no domain-related differences in radar configuration + USE module_configure, ONLY: model_config_rec + USE module_simradar_tools + + IMPLICIT NONE + +! TYPE (model_config_rec), INTENT(in) :: model_config_rec + INTEGER, INTENT(in) :: ids, ide, jds, jde, kds, kde + INTEGER, INTENT(in) :: ims, ime, jms, jme, kms, kme + INTEGER, INTENT(in) :: its, ite, jts, jte, kts, kte + INTEGER, INTENT(in) :: param_first_scalar + INTEGER, INTENT(in) :: p_qv, p_qc, p_qr, p_qs, p_qi, p_qg, p_qh + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(in) :: xlong, xlat + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(in) :: xlong_u, xlat_u + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(in) :: xlong_v, xlat_v + REAL, INTENT(in) :: londbg, latdbg + + INTEGER, INTENT(in) :: nradar, ngrid_radar + INTEGER, INTENT(in) :: e_azimuth, e_angle, e_dist +! REAL, DIMENSION(nradar), INTENT(in) :: gatelon, gatelat, gatehgt, distmax, & +! heightmax, radards +! INTEGER, DIMENSION(nradar), INTENT(in) :: gatetype, gateband, radinterp, em_ray +! REAL, INTENT(in) :: radii_ddist +! REAL, INTENT(in) :: missing_value + ! Compressed dimensions + INTEGER, INTENT(in) :: nrnn, nrnncoord, nrcoord, nr4coord, nncoord + + ! Output + INTEGER, DIMENSION(2), INTENT(out) :: ijdbg + REAL, DIMENSION(e_angle, nradar), INTENT(inout) :: radars_angle + REAL, DIMENSION(e_azimuth, nradar), INTENT(inout) :: radars_azimuth + REAL, DIMENSION(e_dist, nradar), INTENT(inout) :: radars_dist + REAL, INTENT(out) :: ra, rf + CHARACTER(len=256), INTENT(out) :: elipsn + REAL, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(out) :: gatehdist + INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(out) :: inradar + INTEGER, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(out) :: radnum + REAL, DIMENSION(e_angle, e_dist, nradar),INTENT(out) :: radarlons, radarlats, radardist + REAL, DIMENSION(e_azimuth,e_dist,nradar),INTENT(out) :: radarhgts, radarlhgts + REAL,DIMENSION(e_angle,e_dist,nr4coord), INTENT(out) :: radarlon_bnds, radarlat_bnds + REAL, DIMENSION(e_azimuth,e_dist,nr4coord), INTENT(out) & + :: radarhgt_bnds + INTEGER, DIMENSION(ims:ime,jms:jme, nradar), & + INTENT(out) :: ngridradar + INTEGER, DIMENSION(ims:ime, jms:jme, nrnncoord), & + INTENT(out) :: ijgridradar + REAL, DIMENSION(ims:ime, jms:jme, nrnn), INTENT(out) :: ijdistgridradar + + ! LUT space + TYPE(lut_vars), INTENT(out) :: lutgrid + + ! Local + INTEGER :: ir + INTEGER :: ierr + INTEGER :: inrnn, inrnncoord, inncoord, inrcoord, & + inr4coord + INTEGER :: i,j,k,l,n + INTEGER :: idbg, jdbg + INTEGER :: itile_rad, jtile_rad + INTEGER :: e_angle2, e_azimuth2, e_dist2 + INTEGER :: dbg_level + INTEGER, DIMENSION(1) :: klonmin, klatmin, klonmax, klatmax + REAL :: mindist + REAL :: ddlon, ddlat, mindd + CHARACTER(len=256) :: rband + REAL :: rlon, rlat, rhgt, rtype, rddist, rdistx, & + rhgtx + REAL(kind=rr8) :: a, f + REAL :: maxrdistmax, maxrhgtmax, minrds + REAL :: lonSW, lonNW, lonNE, lonSE + REAL :: latSW, latNW, latNE, latSE + REAL :: lonmin, lonmax, latmin, latmax + REAL :: lonrad, latrad, distrad + REAL, DIMENSION(4) :: lonvertices, latvertices + REAL, DIMENSION( ims:ime , jms:jme ) :: dist, ddist, dlon, dlat + REAL, DIMENSION( ims:ime , jms:jme, 4 ) :: bnds_lon, bnds_lat + REAL, DIMENSION( ims:ime ) :: distlon + REAL, DIMENSION( jms:jme ) :: distlat + REAL,DIMENSION(e_angle,e_dist,4) :: radarlon_bnds_NOc, radarlat_bnds_NOc, & + radarhgt_bnds_NOc + INTEGER, DIMENSION(ims:ime, jms:jme, 2, ngrid_radar) :: ijgridradar_NOc + REAL, DIMENSION(ims:ime, jms:jme, ngrid_radar) :: ijdistgridradar_NOc + CHARACTER (LEN=1000) :: diag_message + LOGICAL :: found, inside, dbg + CHARACTER(len=50) :: fname + +!!!!!!!! Variables +! xlong, xlat: matrices of longitudes and latitudes +! xlong_u/v, xlat_u/v: u/vstaggered matrices of longitudes and latitudes +! londbg, latdbg: coordinates of the debugging grid point +! ids, ide, jds, jde, kds, kde: extension of the matrices +! ims, ime, jms, jme, kms, kme: extension of the matrices in memory +! its, ite, jts, jte, kts, kte: extension of the martices by tile +! param_first_scalar index to account for the first active water species (microphysics scheme +! dependend) +! p_qv, p_qc, p_qr, p_qs, p_qi, p_qg, p_qh: moist indices for each specie: vapour, cloud, rain, +! snow, ice, graupel, ice and hail +! idbg, jdbg: point coordinates of the debugging point +!! Radar variables +! nradar: amount of radars +! gatelon: longitude of the radar location +! gatelat: latitude of the radar location +! gatehgt: height of the radar location [m] +! distmax: horizontal range of measuremnt of the radar [km] +! heightmax: vertical range of measuremnt of the radar [km] +! radards: radial resolution of the radar space [m] +! radii_ddist: radii (km) at which arc resolution is the same as minimum ddist (same for all radars) +! em_ray: type of propagation of em ray +! 1: 4/3 Ref +! eloption: type of ellipsoide to use +! 1: 'GRS80 / WGS84 (NAD83)' +! a = 6378137.d0 +! f = 1.d0/298.25722210088d0 +! 2: 'Clarke 1866 (NAD27)' +! a = 6378206.4d0 +! f = 1.d0/294.9786982138d0 +! e_angle: amount of angle beams (horizontal) in the radar space [computed as +! 2*pi*radii_ddist / (min(ddist)) +! e_azimuth: amount of azimuth angles (vertical) in the radar space (same for all radars) +! azimuths: azimuth angles +! IF azimuths(1) == -1; aximuths are internally computed from Nazimuth, by regular intervals +! 30. / (e_azimuth-1) +! e_dist: amount of distances (horizontal) in the radar space [computed as +! 2.*max(distmax) / min(ddist)] +! lutfilen: name of the netCDF file with the look-up-tables for a given miscrophysics of the +! polymetries obntained via the T-matrix derived +! doppler_scheme: mtehod to use to retrieve the radial velocities after LUT +! 0: Reflectivity weighted terminal velocity +! 1: Get PSD integrated backscattering cross section int = v S/S [not implemented] +! ngrid_radar: maximum amount of radar cells in a WRF grid +!! Compressed dimensions +! nrnn = nradar*ngrid_radar +! nrcoord = nradar*2 +! nrnncoord = nradar*ngrid_radar*2 +! nncoord = ngrid_radar*2 +! nr4coord = nradar*4 + + fname = 'simradar_init' + + ! Initializing outputs + inradar = 0 + radnum = 0 + gatehdist = 0. + + ! Consistency checks + !! Collapsed dimensions + + ! nrnn = nradar*ngrid_radar + IF ( nrnn /= nradar * ngrid_radar) THEN + PRINT *, ' nrnn: ', nrnn, ' nradar: ', nradar, ' ngrid_radar: ', ngrid_radar + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, ngrid_radar " // & + "and nrnn do not coincide! nrnn = nradar*ngrid_radar" + CALL wrf_error_fatal ( diag_message ) + END IF + ! nrcoord = nradar*2 + IF ( nrcoord /= nradar * 2) THEN + PRINT *, ' nrcoord: ', nrcoord, ' nradar: ', nradar + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar " // & + "and nrcoord do not coincide! nrcoord = nradar*2" + CALL wrf_error_fatal ( diag_message ) + END IF + ! nrnncoord = nradar*ngrid_radar*2 + IF ( nrnncoord /= nradar * ngrid_radar * 2) THEN + PRINT *, ' nrrnn: ', nrnn, ' nradar: ', nradar, ' ngrid_radar: ', ngrid_radar + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, ngrid_radar " // & + "and nrnncoord do not coincide! nrnn = nradar*ngrid_radar*2" + CALL wrf_error_fatal ( diag_message ) + END IF + ! nncoord = ngrid_radar*2 + IF ( nncoord /= ngrid_radar * 2) THEN + PRINT *, ' nrrnncoord: ', nrnncoord, ' ngrid_radar: ', ngrid_radar + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges ngrid_radar " // & + "and nrnncoord do not coincide! nrnncoord = ngrid_radar*2" + CALL wrf_error_fatal ( diag_message ) + END IF + ! nr4coord = nradar*4 + IF ( nr4coord /= nradar * 4) THEN + PRINT *, ' nrcoord: ', nrcoord, ' nradar: ', nradar + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar " // & + "and nr4coord do not coincide! nrcoord = nradar*4" + CALL wrf_error_fatal ( diag_message ) + END IF + +! IF ( e_rna /= nradar * e_angle) THEN +! PRINT *, ' e_rna: ', e_rna, ' nradar: ', nradar, ' e_angle: ', e_rna +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_angle " // & +! "and e_rna do not coincide! e_rna = nradar*e_angle" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rnh /= nradar * e_azimuth) THEN +! PRINT *, ' e_rnh: ', e_rnh, ' nradar: ', nradar, ' e_azimuth: ', e_azimuth +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_azimuth " // & +! "and e_rnh do not coincide! e_rnh = nradar*e_azimuth" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rnd /= nradar * e_dist) THEN +! PRINT *, ' e_rnd: ', e_rnd, ' nradar: ', nradar, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_dist " // & +! "and e_rnd do not coincide! e_rnd = nradar*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rnad /= nradar * e_angle * e_dist) THEN +! PRINT *, ' e_rnad: ', e_rnad, ' nradar: ', nradar, ' e_angle: ', e_angle, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_angle, e_dist " // & +! "and e_rnad do not coincide! e_rnad = nradar*e_angle*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rnhd /= nradar * e_azimuth * e_dist) THEN +! PRINT *, ' e_rnhd: ', e_rnhd,' nradar: ', nradar,' e_azimuth: ', e_azimuth, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_azimuth, e_dist " //& +! "and e_rnhd do not coincide! e_rnhd = nradar*e_azimuth*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rnahd /= nradar * e_angle * e_azimuth * e_dist) THEN +! PRINT *, ' e_rnhd: ', e_rnahd, ' nradar: ', nradar, ' e_angle: ', e_angle, ' e_azimuth: ', & +! e_azimuth, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges nradar, e_angle, " // & +! "e_azimuth, e_dist and e_rnahd do not coincide! e_rnahd = nradar*e_angle*e_azimuth*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rad /= e_angle * e_dist) THEN +! PRINT *, ' e_rad: ', e_rahd, ' e_angle: ', e_angle, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges e_angle, e_dist " // & +! "and e_rad do not coincide! e_rad = e_angle*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rhd /= e_azimuth * e_dist) THEN +! PRINT *, ' e_rhd: ', e_rhd, ' e_azimuth: ', e_azimuth, ' e_dist: ', e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges e_azimuth, e_dist " // & +! "and e_rhd do not coincide! e_rhd = e_azimuth*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + +! IF ( e_rahd /= e_angle * e_azimuth * e_dist) THEN +! PRINT *, ' e_rahd: ', e_rahd, ' e_angle: ', e_angle, ' e_azimuth: ', e_azimuth, ' e_dist: ', & +! e_dist +! WRITE (diag_message , * ) & +! "simrad.AR: ERROR -- error -- ERROR -- error: individual ranges e_angle, e_azimuth, " // & +! "e_dist and e_rahd do not coincide! e_rahd = e_angle*e_azimuth*e_dist" +! CALL wrf_error_fatal ( diag_message ) +! END IF + + !! Maximum radar coverage + maxrdistmax = MAXVAL(model_config_rec%radar_distmax(1:nradar)) + minrds = MINVAL(model_config_rec%radards(1:nradar)) + + IF ( e_dist /= INT(2*maxrdistmax / minrds)) THEN + PRINT *, ' e_dist: ', e_dist, ' maximum distmax: ', maxrdistmax, ' minimum radards: ', minrds + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: radar range along radial-axis does not " // & + "coincide with the maximum maximum spatial radar range (maximum distmax) and the minimal "// & + "resolution (minimum radards)! e_dist = INT(2*maxrdistmax / minradards)" + CALL wrf_error_fatal ( diag_message ) + END IF + + IF ( e_angle /= INT(2.*pict*model_config_rec%radii_dist*1000./minrds)) THEN + PRINT *, ' e_angle: ', e_angle, ' radii_ddist: ', model_config_rec%radii_dist, & + ' minimum radards: ', minrds + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: amount of range along radial-axis does not " // & + "coincide with the radii at the same arc size as radial resolution (radii_dist) and the " // & + "minimal resolution (minimum radards)! e_angle = INT(2*pi*radii_dist*1000. / minradards)" + CALL wrf_error_fatal ( diag_message ) + END IF + + ! lon, lat debugging point + ! L. Fita, UBA-CIMA-IFAECI. This should be done as in share/wrf_timeseries.F, but it uses + ! a lot of aditional modules, that might not be worth to mess with in the phys/ folder ¿? + ! Minimum grid distance + ! Filling matrices with absurd numbers + DO i=ims, ime + DO j=jms, jme + dlon(i,j) = 10000.*i + 10.*j + dlat(i,j) = 20000.*i + 10.*j + END DO + END DO + dlon(its:ite-1,jts:jte-1) = xlong(its:ite-1,jts:jte-1)-xlong(its+1:ite,jts+1:jte) + dlat(its:ite-1,jts:jte-1) = xlat(its:ite-1,jts:jte-1)-xlat(its+1:ite,jts+1:jte) + + ddist = 1000. + ddist = SQRT(dlon*dlon + dlat*dlat) + mindd = MINVAL(ddist) + + ! Minimum distance to the debugging point + dist = SQRT((xlong-londbg)**2+(xlat-latdbg)**2) + mindist = MINVAL(dist) + + ! Set-up simrad.AR debug + CALL get_wrf_debug_level( dbg_level ) + IF (dbg_level > 75) dbg = .TRUE. + + found = .FALSE. + WRITE(diag_message,*) 'its,ite:', its,ite, ' jts,jte:', jts,jte, ' mindist:', mindist, ' mindd:', & + mindd + CALL wrf_debug(75,diag_message) + IF (mindist < mindd) THEN + DO i=its, ite + DO j=jts, jte + + IF (dist(i,j) == mindist) THEN + idbg = i + jdbg = j + found = .TRUE. + WRITE(diag_message,*)' found debugging grid point _______' + CALL wrf_debug(75,diag_message) + WRITE(diag_message,*)' lon,lat dbg: ', londbg, ' ,', latdbg,' mindist:', mindist + CALL wrf_debug(75,diag_message) + WRITE(diag_message,*)' i,j dbg: ', idbg, jdbg + CALL wrf_debug(75,diag_message) + ijdbg(1) = idbg + ijdbg(2) = jdbg + EXIT + END IF + END DO + IF (found) EXIT + END DO + END IF + + ! Boundaries of the cells + ! SW + bnds_lon(ims:ime,jms:jme,1) = xlong_u(ims:ime,jms:jme) + bnds_lat(ims:ime,jms:jme,1) = xlat_u(ims:ime,jms:jme) + ! NW + bnds_lon(ims:ime,jms:jme,2) = xlong_v(ims:ime,jms+1:jme+1) + bnds_lat(ims:ime,jms:jme,2) = xlat_v(ims:ime,jms+1:jme+1) + ! NE + bnds_lon(ims:ime,jms:jme,3) = xlong_u(ims+1:ime+1,jms:jme) + bnds_lat(ims:ime,jms:jme,3) = xlat_u(ims+1:ime+1,jms:jme) + ! SW + bnds_lon(ims:ime,jms:jme,4) = xlong_v(ims:ime,jms:jme) + bnds_lat(ims:ime,jms:jme,4) = xlat_v(ims:ime,jms:jme) + + ! Getting borders of the tile + lonSW = xlong(its,jts) + lonNW = xlong(its,jte) + lonNE = xlong(ite,jte) + lonSE = xlong(ite,jts) + latSW = xlat(its,jts) + latNW = xlat(its,jte) + latNE = xlat(ite,jte) + latSE = xlat(ite,jts) + + ! WR lon,lat matrices can be rotated, thus matrix vertices might not coincide... + lonvertices = (/ lonSW, lonNW, lonNE, lonSE /) + latvertices = (/ latSW, latNW, latNE, latSE /) + lonmin = MINVAL(lonvertices) + lonmax = MAXVAL(lonvertices) + latmin = MINVAL(latvertices) + latmax = MAXVAL(latvertices) + + klonmin = MINLOC(lonvertices) + klatmin = MINLOC(latvertices) + klonmax = MAXLOC(lonvertices) + klatmax = MAXLOC(latvertices) + + e_angle2 = e_angle / 2 + e_azimuth2 = e_azimuth / 2 + e_dist2 = e_dist / 2 + + ! Getting projection details for the radar space + ! Ellipsoide + IF (model_config_rec%eloption == 1) THEN + a = 6378137.d0 + f = 1.d0/298.25722210088d0 + elipsn = 'GRS80 / WGS84 (NAD83)' + ELSE IF (model_config_rec%eloption == 2) THEN + a = 6378206.4d0 + f = 1.d0/294.9786982138d0 + elipsn = 'Clarke 1866 (NAD27)' + !ELSE IF (config%eloption == 3) THEN + ! CALL elipss (elips) + ELSE + WRITE (diag_message , * ) & + "simrad.AR: ERROR -- error -- ERROR -- error: ellipsoide option eloption = " // & + TRIM(ItoS(model_config_rec%eloption)) // " not ready !! Available ones: 1, 2" + CALL wrf_error_fatal ( diag_message ) + END IF + + ! To real as member of grid + ra = REAL(a) + rf = REAL(f) + + ! Getting actual areas encompassed by a radar + look_radar: DO ir=1, nradar + + rlon = model_config_rec%radar_clon(ir) + rlat = model_config_rec%radar_clat(ir) + rhgt = model_config_rec%radar_chgt(ir) + rband = model_config_rec%radar_band(ir) + rtype = model_config_rec%radar_scan(ir) + rddist = model_config_rec%radards(ir) + rdistx = model_config_rec%radar_distmax(ir) + rhgtx = model_config_rec%radar_hgtmax(ir) + ! Polar coordinates of the radar space + !! + CALL compute_radar_lonlat_radial(e_angle, e_azimuth, e_dist, rlon, rlat, rhgt, rddist, & + rdistx*1000., rhgtx*1000., radars_azimuth(:,ir), a, f, elipsn, radars_angle(:,ir), & + radars_dist(:,ir), radars_azimuth(:,ir), radarlons(:,:,ir), radarlats(:,:,ir), & + radarhgts(:,:,ir), radarlhgts(:,:,ir), radarlon_bnds_NOc, radarlat_bnds_NOc, & + radarhgt_bnds_NOc, dbg) + + ! Checking pertinence of the radar into this tile + CALL pertinence_polar(ims, ime, jms, jme, xlong, xlat, bnds_lon, bnds_lat, 1, e_dist, 1, & + e_angle, radarlons(:,:,ir), radarlats(:,:,ir), rlon, rlat, rdistx, & + model_config_rec%missing_value, ngrid_radar, gatehdist(:,:,ir), inradar, radnum(:,:,ir), & + ngridradar(:,:,ir), 2, ijgridradar_NOc, ijdistgridradar_NOc, inside, dbg) + + ! Compressing results + DO i=ims, ime + DO j=jms, jme + DO n=1, ngridradar(i,j,ir) + inrnn = (ir-1)*ngrid_radar + n + ijdistgridradar(i,j,inrnn) = ijdistgridradar_NOc(i,j,n) + DO l=1, 2 + inrnncoord = (ir-1)*ngrid_radar*2 + (n-1)*2 + l + ijgridradar(i,j,inrnncoord) = ijgridradar_NOc(i,j,n,l) + END DO + END DO + + DO l=1, 4 + inr4coord = (ir-1)*4 + l + radarlon_bnds(i,j,inr4coord) = radarlon_bnds_NOc(i,j,l) + radarlat_bnds(i,j,inr4coord) = radarlat_bnds_NOc(i,j,l) + radarhgt_bnds(i,j,inr4coord) = radarhgt_bnds_NOc(i,j,l) + END DO + END DO + END DO + + END DO look_radar + + !! Looking for the specees to interpolate and use + !! When simulating radar, there is an array that holds the variable to interpolate + !! cvals(:,1:12) = (/ qv, qc, qr, qi, qs, qg, qh, ua, va, wa, ta, press /) + !! + !! Not all micrphysics schemes have all the water species. This is managed via PARAM_FIRST_SCALAR + !! Extracted from phys/module_radiation_driver.F + ! Among the arguments is moist, a four-dimensional scalar array storing + ! a variable number of moisture tracers, depending on the physics + ! configuration for the WRF run, as determined in the namelist. The + ! highest numbered index of active moisture tracers the integer argument + ! n_moist (note: the number of tracers at run time is the quantity + ! n_moist - PARAM_FIRST_SCALAR + 1 , not n_moist. Individual tracers + ! may be indexed from moist by the Registry name of the tracer prepended + ! with P_; for example P_QC is the index of cloud water. An index + ! represents a valid, active field only if the index is greater than + ! or equal to PARAM_FIRST_SCALAR. PARAM_FIRST_SCALAR and the individual + ! indices for each tracer is defined in module_state_description and + ! set in set_scalar_indices_from_config defined in + ! frame/module_configure.F. + !! + ! + ! keeping in a vector which indices have to be computed + lutgrid%N_tocompute = 0 + lutgrid%Nmoist = 0 + lutgrid%tocompute = 0 + lutgrid%imoist = 0 + IF (p_qv .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 1 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qv + END IF + IF (p_qc .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 2 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qc + END IF + IF (p_qr .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 3 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qr + END IF + IF (p_qi .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 4 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qi + END IF + IF (p_qs .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 5 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qs + END IF + IF (p_qg .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 6 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qg + END IF + IF (p_qh .ge. param_first_scalar) THEN + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 7 + lutgrid%Nmoist = lutgrid%Nmoist + 1 + lutgrid%imoist(lutgrid%Nmoist) = p_qh + END IF + ! These will be always interpolated / used + !! ua, va, wa, ta, press + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 8 + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 9 + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 10 + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 11 + lutgrid%N_tocompute = lutgrid%N_tocompute + 1 + lutgrid%tocompute(lutgrid%N_tocompute) = 12 + + !! Look-up tables space for T-matrix + CALL read_lut(model_config_rec%lutfilen, lutgrid, dbg) + + ! From namelist or from file? + ! Let's put all into lutgrid for simplicity (but, duplicity) + lutgrid%q_rain_min = model_config_rec%q_rain_min + lutgrid%q_rain_max = model_config_rec%q_rain_max + lutgrid%qn_rain_min = model_config_rec%qn_rain_min + lutgrid%qn_rain_max = model_config_rec%qn_rain_min + lutgrid%q_snow_min = model_config_rec%q_snow_min + lutgrid%q_snow_max = model_config_rec%q_snow_min + lutgrid%qn_snow_min = model_config_rec%q_snow_min + lutgrid%qn_snow_max = model_config_rec%q_snow_min + lutgrid%q_grau_min = model_config_rec%q_grau_min + lutgrid%q_grau_max = model_config_rec%q_grau_min + lutgrid%qn_grau_min = model_config_rec%q_grau_min + lutgrid%qn_grau_max = model_config_rec%q_grau_min + lutgrid%q_ice_min = model_config_rec%q_ice_min + lutgrid%q_ice_max = model_config_rec%q_ice_min + lutgrid%qn_ice_min = model_config_rec%q_ice_min + lutgrid%qn_ice_max = model_config_rec%q_ice_min + lutgrid%q_clou_min = model_config_rec%q_clou_min + lutgrid%q_clou_max = model_config_rec%q_clou_min + lutgrid%dim_ice = model_config_rec%dim_ice + lutgrid%dimv = model_config_rec%Nqxbins + + ALLOCATE(lutgrid%q_rain_vec(model_config_rec%Nqxbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'q_rain_vec', 1, (/model_config_rec%Nqxbins/), fname, ierr) + ALLOCATE(lutgrid%q_snow_vec(model_config_rec%Nqxbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'q_snow_vec', 1, (/model_config_rec%Nqxbins/), fname, ierr) + ALLOCATE(lutgrid%q_grau_vec(model_config_rec%Nqxbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'q_grau_vec', 1, (/model_config_rec%Nqxbins/), fname, ierr) + ALLOCATE(lutgrid%q_ice_vec(model_config_rec%Nqxbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'q_ice_vec', 1, (/model_config_rec%Nqxbins/), fname, ierr) + + CALL create_logarithmic_scaled_vectorsD(model_config_rec%q_rain_min, model_config_rec%q_rain_max, & + model_config_rec%Nqxbins, lutgrid%q_rain_vec, dbg) + CALL create_logarithmic_scaled_vectorsD(model_config_rec%q_snow_min, model_config_rec%q_snow_max, & + model_config_rec%Nqxbins, lutgrid%q_snow_vec, .FALSE.) + CALL create_logarithmic_scaled_vectorsD(model_config_rec%q_grau_min, model_config_rec%q_grau_max, & + model_config_rec%Nqxbins, lutgrid%q_grau_vec, .FALSE.) + CALL create_logarithmic_scaled_vectorsD(model_config_rec%q_ice_min, model_config_rec%q_ice_max, & + model_config_rec%Nqxbins, lutgrid%q_ice_vec, .FALSE.) + + RETURN + + END SUBROUTINE simradar_init +#endif + END MODULE module_physics_init + diff --git a/phys/module_simradar.F b/phys/module_simradar.F new file mode 100644 index 0000000000..ec1641f494 --- /dev/null +++ b/phys/module_simradar.F @@ -0,0 +1,1043 @@ +MODULE module_simradar +#ifdef SIMRADAR +! Module for the simulation of surface polimeric radar following simrad.AR, Galligani et al (CONGREMET XIII), 2018 + + INTEGER, PARAMETER :: Ntotvarsinterp = 12 + + CONTAINS + + SUBROUTINE simradar_calc( & + is_restart, is_output_time, & + nsteps, & + dt, xtime, curr_secs2, & + + timedbg, & + londbg, latdbg, & + idbg, jdbg, & + + ! Input + n_moist, & + moist, & + prm_first_scalar, & + m_qv, m_qc, m_qr, m_qs, m_qi, m_qg, m_qh, & + m_qng, m_qnh, m_qnr, m_qvolg, & +! qv, qc, qr, qs, qi, qg, qh, & + lon, lat, & + sina, cosa, & + u, v, w, & + p, pb, & + ta, ph, phb, & + nradar, ngrid_radar, & + e_angle, e_azimuth, e_dist, & + coord2, coord3, coord4, around2, around3, & + ngrid3_radar, nrd, & + nrnn, nrnncoord, nrcoord, nr4coord, & + elipsn, a, f, & + radars_angle, radars_azimuth, radars_dist, & + ijradardbg, & + radarlons, radarlats, & + radarhgts, radarlhgts, & + gatehdist, inradar, radnum, & + ngridradar, ijgridradar, ijdistgridradar, & + nbands, nlutqxbins, & + nluticebins, nlutlevs, & + elevations, & + lutzh_rain, lutzdr_rain, & +! lutldr_rain, lutaih_rain, & +! lutaiv_rain, lutkdp_rain, & + lutkdp_rain, & + lutzh_snow, lutzdr_snow, & +! lutldr_snow, lutaih_snow, & +! lutaiv_snow, lutkdp_snow, & + lutkdp_snow, & + lutzh_grau, lutzdr_grau, & +! lutldr_grau, lutaih_grau, & +! lutaiv_grau, lutkdp_grau, & + lutkdp_grau, & + lutzh_ice, lutzdr_ice, & +! lutldr_ice, lutaih_ice, & +! lutaiv_ice, lutkdp_ice, & + lutkdp_ice, & + lutzh_hail, lutzdr_hail, & +! lutldr_hail, lutaih_hail, & +! lutaiv_hail, lutkdp_hail, & + lutkdp_hail, & + n_tocompute, tocompute, & + nmoist, imoist, & + Nqxbins, & + q_rain_vec, q_snow_vec, q_grau_vec, & + q_ice_vec, q_hail_vec, & + lutgrid, & + + ! Output + !!! Optional. Modify Registry/registry.simradar and recompile + !!SIMRADARop!! iqvm, iqcm, iqrm, iqim, & + !!SIMRADARop!! iqsm, iqgm, iqhm, & + !!SIMRADARop!! iuam, ivam, iwam, itam, ipressm, & + zh_rain, zdr_rain, kdp_rain, & + zh_snow, zdr_snow, kdp_snow, & + zh_grau, zdr_grau, kdp_grau, & + zh_ice, zdr_ice, kdp_ice, & + zh_hail, zdr_hail, kdp_hail, & + zh_total, zdr_total, kdp_total, & + tv_rain, tv_snow, tv_grau, & + tv_ice, tv_hail, & + tv_wgt, vr, & + + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + ips,ipe, jps,jpe, kps,kpe, & ! patch dims + i_start,i_end,j_start,j_end,kts,kte,num_tiles & + ) + +!---------------------------------------------------------------------- + ! Main subroutine to compute simrad.AR variables + ! NOTE radars are in polar space: + ! dx: distance dy: angle dz: azimuth + ! 3Dvariables (dy,dz,dx), lon,lats(dy,dx), hgt(dx,dz) + + ! There is no domain-related differences in radar configuration + USE module_configure, ONLY: model_config_rec + USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval + USE module_configure + USE module_simradar_tools + USE module_model_constants +! USE module_diagvar_simradar + + IMPLICIT NONE + +!====================================================================== +! Definitions +!----------- +!-- DT time step (second) +!-- XTIME forecast time +!-- curr_secs2 current time in seconds since simulation restart +!-- radt frequency of radiative calls (minutes) +! +!-- ids start index for i in domain +!-- ide end index for i in domain +!-- jds start index for j in domain +!-- jde end index for j in domain +!-- kds start index for k in domain +!-- kde end index for k in domain +!-- ims start index for i in memory +!-- ime end index for i in memory +!-- jms start index for j in memory +!-- jme end index for j in memory +!-- ips start index for i in patch +!-- ipe end index for i in patch +!-- jps start index for j in patch +!-- jpe end index for j in patch +!-- kms start index for k in memory +!-- kme end index for k in memory +!-- i_start start indices for i in tile +!-- i_end end indices for i in tile +!-- j_start start indices for j in tile +!-- j_end end indices for j in tile +!-- kts start index for k in tile +!-- kte end index for k in tile +!-- num_tiles number of tiles +!-- n_moist number of water species +!-- prm_first_scalar index to account for the first active water species (microphysics scheme +! dependend) +!-- m_qv, m_qc, m_qr, m_qs, m_qi, m_qg, m_qh: moist indices for each specie: vapour, cloud, rain, +! snow, ice, graupel, ice and hail +!! Extracted from phys/module_radiation_driver.F +! Among the arguments is moist, a four-dimensional scalar array storing +! a variable number of moisture tracers, depending on the physics +! configuration for the WRF run, as determined in the namelist. The +! highest numbered index of active moisture tracers the integer argument +! n_moist (note: the number of tracers at run time is the quantity +! n_moist - prm_first_scalar + 1 , not n_moist. Individual tracers +! may be indexed from moist by the Registry name of the tracer prepended +! with P_; for example P_QC is the index of cloud water. An index +! represents a valid, active field only if the index is greater than +! or equal to prm_first_scalar. prm_first_scalar and the individual +! indices for each tracer is defined in module_state_description and +! set in set_scalar_indices_from_config defined in +! frame/module_configure.F. +!! +! +!----------- +! is_restart: whether if simulation is a restart +! is_output_time: whether if this iteration is an output time-step + +! ------- INPUT vars ------- +! lon: longitude +! lat: latitude +! moist: 4-dimensional array with all water species [kgkg-1] +!! qv: water vapour mixing ratio [kgkg-1] +!! qc: cloud mixing ratio [kgkg-1] +!! qr: rain mixing ratio [kgkg-1] +!! qs: snow mixing ratio [kgkg-1] +!! qi: ice mixing ratio [kgkg-1] +!! qg: graupel mixing ratio [kgkg-1] +!! qh: hail mixing ratio [kgkg-1] +! sina: local sine of map rotation [1] +! cosa: local cosine of map rotation [1] +! u: eastward wind speed [ms-1] +! v: northward wind speed [ms-1] +! wa: vertical wind speed [ms-1] +! p: 3D pressure perturbation [Pa] +! pb: 3D pressure base state [Pa] +! ta: air temperature [K] +! ph: 3D gepotential perturbation [m2s-2] +! phb: 3D gepotential base state [m2s-2] + +!! Radar variables +! nradar: amount of radars +! e_angle: amount of angle beams (horizontal) in the radar space +! e_azimuth: amount of azimuth angles (vertical) in the radar space (same for all radars) +! azimuths: azimuth angles +! e_dist: amount of distances (horizontal) in the radar space +! coord2: amount of 2-dimensional coordinates (2) +! coord3: amount of 3-dimensional coordinates (3) +! coord4: amount of coordinates for grid/cell 2-dimensional boundaries (4) +! around2: amoount of points of 3-dimensional surroundings (9) +! around3: amoount of points of 3-dimensional surroundings (27) +! ngrid_radar: maximum amount of radar cells in a WRF grid +! ngrid3_radar: maximum amount of 3-dimensional radar cells in a WRF grid +! inradar: pertinence of the grid cell to any radar (1: yes) +! radarnum: pertinence of the grid cell inside a given radar area [1] +! radars_dist: ground radar gate distances to their center [m] +! radars_angle: ground radar gate horizontal angles to emit beams [m] +! radars_azimuth: ground radar gate vertical angles to emit beams [m] +! radarlons: ground radar gate longitude [degree_east] +! radarlats: ground radar gate latitude [degree_north] +! radarhgts: ground radar gate height [m] +! radarlhgts: ground radar gate local height [m] +! radarhdist: ground radar horizontal distance of the grid point [m] +! gatehdist: distacnce of the grid points to the position of the radar [m] +! ngridradar: amount of radar cells inside the grid point [1] +! ijgridradar: coordinates of radar cells inside the grid point [-] +! ijdistgridradar: distance of radar cells inside the grid point to its center [m] + +!! T-MATRIX LUT variables +! nbands: amount of radar bans +! nlutqxbins: mnount of mixing ratio bins +! nluticebins: mnount of ice mixing ratio bins +! nlutlevs: amount of elevations +! elevations: elevations [degree] +! lutzh_rain: Zh for rain by band, elevation and bind accordig to a microphysics scheme +! lutzdr_rain: Zdr for rain by band, elevation and bind accordig to a microphysics scheme +! lutldr_rain: LDR for rain by band, elevation and bind accordig to a microphysics scheme +! lutaih_rain: Aih for rain by band, elevation and bind accordig to a microphysics scheme +! lutaiv_rain: Aiv for rain by band, elevation and bind accordig to a microphysics scheme +! lutkdp_rain: KDP for rain by band, elevation and bind accordig to a microphysics scheme +! lutzh_snow: Zh for snow by band, elevation and bind accordig to a microphysics scheme +! lutzdr_snow: Zdr for snow by band, elevation and bind accordig to a microphysics scheme +! lutldr_snow: LDR for snow by band, elevation and bind accordig to a microphysics scheme +! lutaih_snow: Aih for snow by band, elevation and bind accordig to a microphysics scheme +! lutaiv_snow: Aiv for snow by band, elevation and bind accordig to a microphysics scheme +! lutkdp_snow: KDP for snow by band, elevation and bind accordig to a microphysics scheme +! lutzh_grau: Zh for graupel by band, elevation and bind accordig to a microphysics scheme +! lutzdr_grau: Zdr for graupel by band, elevation and bind accordig to a microphysics scheme +! lutldr_grau: LDR for graupel by band, elevation and bind accordig to a microphysics scheme +! lutaih_grau: Aih for graupel by band, elevation and bind accordig to a microphysics scheme +! lutaiv_grau: Aiv for graupel by band, elevation and bind accordig to a microphysics scheme +! lutkdp_grau: KDP for graupel by band, elevation and bind accordig to a microphysics scheme +! lutzh_ice: Zh for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutzdr_ice: Zdr for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutldr_ice: LDR for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutaih_ice: Aih for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutaiv_ice: Aiv for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutkdp_ice: KDP for cloud ice by band, elevation and bind accordig to a microphysics scheme +! lutzh_hail: Zh for hail by band, elevation and bind accordig to a microphysics scheme +! lutzdr_hail: Zdr for hail by band, elevation and bind accordig to a microphysics scheme +! lutldr_hail: LDR for hail by band, elevation and bind accordig to a microphysics scheme +! lutaih_hail: Aih for hail by band, elevation and bind accordig to a microphysics scheme +! lutaiv_hail: Aiv for hail by band, elevation and bind accordig to a microphysics scheme +! lutkdp_hail: KDP for hail by band, elevation and bind accordig to a microphysics scheme +! n_tocompute: amount of variables to interpolate / use in the radar space +! tocompute: indices of the variables to interpolate / use in the radar space +! nmoist: amount of water species to interpolate / use in the radar space +! imoist: indices of the of water species to interpolate / use in the radar space + +!! Microphysics depending variables +! Nqxbins: amount of bins for the water species for the interpolation into LUT space +! q_rain_vec: bins for rain for the interpolation into LUT space +! q_snow_vec: bins for snow for the interpolation into LUT space +! q_grau_vec: bins for graupel for the interpolation into LUT space +! q_ice_vec, q_hail_vec: bins for cloud ice for the interpolation into LUT space +! q_ice_vec, q_hail_vec: bins for hail for the interpolation into LUT space + +!! Compressed dimensions +! nrnn = nradar*ngrid_radar +! nrcoord = nradar*coord2 +! nrnncoord = nradar*ngrid_radar*coord2 +! nncoord = ngrid_radar*coord2 +! nr4coord = nradar*coord4 +! nrd: nradar*e_dist + +! timedbg: frecuency in minutes of time-debugging +! londbg: longitude in degrees of grid-point debugging (-999., not used) +! latdbg: latitude in degrees of grid-point debugging (-999., not used) + +! Output +!------- +!! Optional. !!SIMRADARop!! +!! Modify Registry/registry.simradar and recompile +! iqvm: 3-dimensional radar interpolated variable of water vapour mixing ratio [kgkg-1] +! iqcm: 3-dimensional radar interpolated variable of cloud mixing ratio [kgkg-1] +! iqrm: 3-dimensional radar interpolated variable of rain mixing ratio [kgkg-1] +! iqim: 3-dimensional radar interpolated variable of ice mixing ratio [kgkg-1] +! iqsm: 3-dimensional radar interpolated variable of snow mixing ratio [kgkg-1] +! iqgm: 3-dimensional radar interpolated variable of graupel mixing ratio [kgkg-1] +! iqhm: 3-dimensional radar interpolated variable of hail mixing ratio [kgkg-1] +! iuam: 3-dimensional radar interpolated variable of eastward wind speed [ms-1] +! ivam: 3-dimensional radar interpolated variable of northward wind speed [ms-1] +! iwam: 3-dimensional radar interpolated variable of vertical wind speed [ms-1] +! itam: 3-dimensional radar interpolated variable of air temperature [K] +! ipressm: 3-dimensional radar interpolated variable of air pressure [Pa] +!! END Optional +! zh_rain: Reflectivity for Rain [dBZ] +! zdr_rain: ZDR for Rain [dBz] +! kdp_rain: KDP for Rain [degree/km] +! zh_snow: Reflectivity for Snow [dBZ] +! zdr_snow: ZDR for Snow [dBz] +! kdp_snow: KDP for Snow [degree/km] +! zh_grau: Reflectivity for Graupel [dBZ] +! zdr_grau: ZDR for Graupel [dBz] +! kdp_grau: KDP for Graupel [degree/km] +! zh_ice: Reflectivity for cloud Ice [dBZ] +! zdr_ice: ZDR for cloud Ice [dBz] +! kdp_ice: KDP for cloud Ice [degree/km] +! zh_hail: Reflectivity for Hail [dBZ] +! zdr_hail: ZDR for Hail [dBz] +! kdp_hail: KDP for Hail [degree/km] +! zh_total: Reflectivity for all species [dBZ] +! zdr_total: ZDR for all species [dBz] +! kdp_total: KDP for all species [degree/km] +! tv_rain: terminal velocity for Rain [ms-1] +! tv_snow: terminal velocity for Snow [ms-1] +! tv_grau: terminal velocity for Graupel [ms-1] +! tv_ice: terminal velocity for cloud Ice [ms-1] +! tv_hail: terminal velocity for Hail [ms-1] +! tv_wgt: weighted terminal velocity [ms-1] +! vr: Doppler velocity [ms-1] + +!====================================================================== + + INTEGER, INTENT(in) :: & + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + ips,ipe, jps,jpe, kps,kpe, & + kts,kte, & + num_tiles, & + n_moist + + + INTEGER, DIMENSION(num_tiles), INTENT(in) :: i_start, i_end, j_start, j_end + INTEGER, INTENT(in) :: prm_first_scalar + INTEGER, INTENT(in) :: m_qv, m_qc, m_qr, m_qs, m_qi, m_qg, m_qh, & + m_qng, m_qnh, m_qnr, m_qvolg + REAL, DIMENSION(ims:ime,jms:jme), INTENT(in) :: lon, lat + REAL, DIMENSION(ims:ime,jms:jme), INTENT(in) :: sina, cosa + REAL, DIMENSION(ims:ime,kms:kme,jms:jme,n_moist), & + INTENT(in) :: moist +! REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, & +! INTENT(in) :: qv, qc, qr, qs, qi, qg, qh + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(in) :: u, v, w, p, pb, ta, ph, phb + + REAL, INTENT(in) :: DT, XTIME + REAL, INTENT(in) :: curr_secs2 + LOGICAL, INTENT(in) :: is_restart, is_output_time + INTEGER, INTENT(in) :: timedbg + REAL, INTENT(in) :: londbg, latdbg + INTEGER, INTENT(in) :: idbg, jdbg + + INTEGER, INTENT(inout) :: nsteps + INTEGER, INTENT(in) :: coord2, coord3, coord4, around2, around3 + +! Specific Radar variables + INTEGER, INTENT(in) :: nradar, ngrid_radar + INTEGER, INTENT(in) :: e_azimuth + INTEGER, INTENT(in) :: e_angle, e_dist +! REAL, DIMENSION(nradar), INTENT(in) :: gatelon, gatelat, gatehgt, distmax, & +! heightmax, radards +! INTEGER, DIMENSION(nradar), INTENT(in) :: gatetype, gateband, radinterp, em_ray +! INTEGER, DIMENSION(nradar), INTENT(in) :: em_ray +! REAL, INTENT(in) :: radii_ddist +! REAL, INTENT(in) :: missing_value + REAL, DIMENSION(e_angle, nradar), INTENT(in) :: radars_angle + REAL, DIMENSION(e_azimuth, nradar), INTENT(in) :: radars_azimuth + REAL, DIMENSION(e_dist, nradar), INTENT(in) :: radars_dist + REAL, INTENT(in) :: a, f + CHARACTER(len=30), INTENT(in) :: elipsn + ! Compressed dimensions + INTEGER, INTENT(in) :: nrnn, nrnncoord, nrcoord, nr4coord + INTEGER, INTENT(in) :: ngrid3_radar, nrd + INTEGER, DIMENSION(coord2), INTENT(in) :: ijradardbg + REAL, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(in) :: gatehdist + INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(in) :: inradar + INTEGER, DIMENSION(ims:ime, jms:jme, nradar), & + INTENT(in) :: radnum + REAL, DIMENSION(e_angle, e_dist, nradar), INTENT(in) :: radarlons, radarlats + REAL, DIMENSION(e_azimuth,e_dist,nradar), INTENT(in) :: radarhgts, radarlhgts +! REAL,DIMENSION(e_angle,e_dist,nr4coord), INTENT(in) :: radarlons_bnds, radarlat_bnds, & +! REAL, DIMENSION(e_azimuth,e_dist,nr4coord), INTENT(in) & +! :: radarhgt_bnds + INTEGER, DIMENSION(ims:ime,jms:jme, nradar), & + INTENT(in) :: ngridradar + INTEGER, DIMENSION(ims:ime, jms:jme, nrnncoord), & + INTENT(in) :: ijgridradar + REAL, DIMENSION(ims:ime, jms:jme, nrnn), INTENT(in) :: ijdistgridradar + + ! LUT space + INTEGER, INTENT(in) :: nbands, nlutqxbins, nluticebins, nlutlevs + REAL, DIMENSION(nlutlevs), INTENT(in) :: elevations + REAL, DIMENSION(nbands,nlutlevs,nlutqxbins), & + INTENT(in) :: lutzh_rain, lutzdr_rain, lutkdp_rain, & + lutzh_snow, lutzdr_snow, lutkdp_snow, & + lutzh_grau, lutzdr_grau, lutkdp_grau, & + lutzh_ice, lutzdr_ice, lutkdp_ice, & + lutzh_hail, lutzdr_hail, lutkdp_hail +! lutaih_rain, lutaiv_rain, lutkdp_rain, & +! lutzh_snow, lutzdr_snow, lutldr_snow, lutaih_snow, lutaiv_snow, lutkdp_snow, & +! lutzh_grau, lutzdr_grau, lutldr_grau, lutaih_grau, lutaiv_grau, lutkdp_grau, & +! lutzh_ice, lutzdr_ice, lutldr_ice, lutaih_ice, lutaiv_ice, lutkdp_ice, & +! lutzh_hail, lutzdr_hail, lutldr_hail, lutaih_hail, lutaiv_hail, lutkdp_hail + INTEGER, INTENT(in) :: n_tocompute, nmoist + INTEGER, DIMENSION(12), INTENT(in) :: tocompute, imoist + INTEGER, INTENT(in) :: Nqxbins + REAL, DIMENSION(Nqxbins), INTENT(in) :: q_rain_vec, q_snow_vec, q_grau_vec, & + q_ice_vec, q_hail_vec + TYPE(lut_vars), INTENT(in) :: lutgrid + + ! Output + !! optional !! SIMRADARop. Uncomment, modify Regisry/registry.simradar and recompile + !!REAL, DIMENSION(e_angle,e_azimiuth,nrd), INTENT(out) & + !! :: iqvm, iqcm, iqrm, iqim, iqsm, iqgm, iqhm, & + !! iuam, ivam, iwam, itam, ipressm + + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_rain, zdr_rain, kdp_rain + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_snow, zdr_snow, kdp_snow + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_grau, zdr_grau, kdp_grau + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_ice, zdr_ice, kdp_ice + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_hail, zdr_hail, kdp_hail + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: zh_total, zdr_total, kdp_total + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: tv_rain, tv_snow, tv_grau + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: tv_ice, tv_hail + REAL, DIMENSION(e_angle, e_azimuth, nrd),INTENT(out) :: tv_wgt, vr + + ! LOCAL VAR + INTEGER :: i, j, k, l, m, n, its, ite, jts, jte, ij + INTEGER :: idp, jdp, dimz, dimz1, dz1 + INTEGER :: i1, i_1, j1, j_1, k1, k2, tdbg + INTEGER :: ii1, ei1, ij1, ej1 + INTEGER, DIMENSION(num_tiles) :: i2, j2 + INTEGER :: inrnn, inrnncoord, inncoord, inrcoord, & + inr4coord, inrd + + REAL :: xtimep + REAL :: prv + + LOGICAL, EXTERNAL :: wrf_dm_on_monitor + CHARACTER(len=256) :: msg + LOGICAL :: ijdbg, cdone + INTEGER :: dbg_level + + INTEGER :: ir, ia, ih, id, iband + INTEGER :: ix, iy, iz + INTEGER, DIMENSION(ngrid_radar, coord2) :: ijgridradar_NOc + REAL, DIMENSION(ngrid_radar, nradar) :: ijdistgridradar_NOc + + ! radar 3-dimensional + INTEGER :: ijkNtot + INTEGER, DIMENSION(kms:kme) :: ijkngridradar + INTEGER, DIMENSION(kms:kme, coord3, ngrid3_radar) :: ijkgridradar + REAL, DIMENSION(kms:kme, around3, ngrid3_radar) :: wghts + + INTEGER :: ii2, ii3, iqm, iicv, ic + INTEGER :: ei2, ei3 + INTEGER :: i123 + INTEGER :: ijkfound + REAL, DIMENSION(around3,Ntotvarsinterp) :: cvals + CHARACTER(len=1), DIMENSION(Nbands) :: availrband + + REAL, DIMENSION(ims:ime,kms:kme,jms:jme) :: zg, press, unua, unva, ua, va + + REAL :: iqv, iqc, iqr, iqi, iqs, iqg, iqh, iua, & + iva, iwa, ita, ipress + REAL(kind=rr8) :: dqv, dqc, dqr, dqi, dqs, dqg, dqh, dua, & + dva, dwa, dta, dpress + REAL(kind=rr8) :: dZh_q, dZdr_q, dKDP_q + REAL :: zh_rain_NOc, zdr_rain_NOc, kdp_rain_NOc + REAL :: zh_snow_NOc, zdr_snow_NOc, kdp_snow_NOc + REAL :: zh_grau_NOc, zdr_grau_NOc, kdp_grau_NOc + REAL :: zh_ice_NOc, zdr_ice_NOc, kdp_ice_NOc + REAL :: zh_hail_NOc, zdr_hail_NOc, kdp_hail_NOc + REAL :: zh_total_NOc, zdr_total_NOc, kdp_total_NOc + REAL(kind=rr8) :: dtv_rain, dtv_snow, dtv_grau, dtv_ice, & + dtv_hail + REAL :: tv_rain_NOc, tv_snow_NOc, tv_grau_NOc + REAL :: tv_ice_NOc, tv_hail_NOc + REAL(kind=rr8) :: dangle, dazimuth, ddist + REAL(kind=rr8) :: cos_el, sin_el, cos_an, sin_an + REAL(kind=rr8) :: denom + REAL(kind=rr8) :: dwt, dvr + REAL :: tv_wgt_NOc, vr_NOc, wt_NOc + LOGICAL :: dbg, printdbg + CHARACTER(len=S32) :: fname + + fname = 'simradar_calc' + + ! Showing version of the module + IF (curr_secs2 == 0. ) THEN + WRITE(msg,*) ' simrad.AR module version: ' // TRIM(simradARversion) // ' ...' + PRINT *, TRIM(msg) + CALL wrf_debug(0,msg) + END IF + + ! timedbg + ! when time-debug in minutes + tdbg = MOD(xtime, timedbg*1.) + + ! Getting the overall debug level + CALL get_wrf_debug_level( dbg_level ) + IF (dbg_level > 75) dbg = .TRUE. + cdone = .FALSE. + + ! L. Fita: NOTE about efficiency + ! I am not sure if this is the more efficient way: + ! 1.- IF --> 2.- repeat DO + ! 1.- DO --> 2.- IF + + dimz = kme - kms + 1 + dimz1 = dimz + 1 + dz1 = dimz - 1 + + IF (ALL(inradar == 0)) THEN + !!! Optional. Modify Registry/registry.simradar and recompile + !!SIMRADARop!!iqvm = model_config_rec%missing_value + !!SIMRADARop!!iqcm = model_config_rec%missing_value + !!SIMRADARop!!iqrm = model_config_rec%missing_value + !!SIMRADARop!!iqim = model_config_rec%missing_value + !!SIMRADARop!!iqsm = model_config_rec%missing_value + !!SIMRADARop!!iqgm = model_config_rec%missing_value + !!SIMRADARop!!iqhm = model_config_rec%missing_value + !!SIMRADARop!!iuam = model_config_rec%missing_value + !!SIMRADARop!!ivam = model_config_rec%missing_value + !!SIMRADARop!!iwam = model_config_rec%missing_value + !!SIMRADARop!!itam = model_config_rec%missing_value + !!SIMRADARop!!ipressm = model_config_rec%missing_value + + zh_rain = model_config_rec%missing_value + zdr_rain = model_config_rec%missing_value + kdp_rain = model_config_rec%missing_value + zh_snow = model_config_rec%missing_value + zdr_snow = model_config_rec%missing_value + kdp_snow = model_config_rec%missing_value + zh_grau = model_config_rec%missing_value + zdr_grau = model_config_rec%missing_value + kdp_grau = model_config_rec%missing_value + zh_ice = model_config_rec%missing_value + zdr_ice = model_config_rec%missing_value + kdp_ice = model_config_rec%missing_value + zh_hail = model_config_rec%missing_value + zdr_hail = model_config_rec%missing_value + kdp_hail = model_config_rec%missing_value + zh_total = model_config_rec%missing_value + zdr_total = model_config_rec%missing_value + kdp_total = model_config_rec%missing_value + tv_rain = model_config_rec%missing_value + tv_snow = model_config_rec%missing_value + tv_grau = model_config_rec%missing_value + tv_ice = model_config_rec%missing_value + tv_hail = model_config_rec%missing_value + tv_wgt = model_config_rec%missing_value + vr = model_config_rec%missing_value + RETURN + ELSE + !! Getting standard atmospheric variables + DO ij = 1 , num_tiles + DO i = i_start(ij), i_end(ij) + DO j = j_start(ij), j_end(ij) + zg(i,:,j) = (ph(i,:,j) + phb(i,:,j))/9.81 + press(i,:,j) = p(i,:,j) + pb(i,:,j) + unua(i,:,j) = 0.5*(u(i+1,:,j) + u(i,:,j)) + unva(i,:,j) = 0.5*(u(i,:,j+1) + u(i,:,j)) + ua(i,:,j) = unua(i,:,j)*cosa(i,j) - unva(i,:,j)*sina(i,j) + va(i,:,j) = unua(i,:,j)*cosa(i,j) + unva(i,:,j)*sina(i,j) + END DO + END DO + END DO + END IF + + ! For test purposes + IF ((londbg == -999.) .OR. (latdbg == -999.)) THEN + !$OMP PARALLEL DO & + !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + i2(ij) = INT(0.5*(i_start(ij) + i_end(ij))) + j2(ij) = INT(0.5*(j_start(ij) + j_end(ij))) + END DO + !$OMP END PARALLEL DO + ELSE + !$OMP PARALLEL DO & + !$OMP PRIVATE ( ij ) + DO ij = 1 , num_tiles + i2(ij) = idbg + j2(ij) = jdbg + END DO + !$OMP END PARALLEL DO + END IF + k2 = 0.5*(kms+kme) + + ! Computing variables + ijdbg = .FALSE. + + !! Looping for all radars + radars_loop: DO ir=1, nradar + + ijkfound = 0 + + IF (TRIM(model_config_rec%radar_band(ir)) == 'S') THEN + iband = 1 + ELSE IF (TRIM(model_config_rec%radar_band(ir)) == 'C') THEN + iband = 2 + ELSE IF (TRIM(model_config_rec%radar_band(ir)) == 'X') THEN + iband = 3 + ELSE + availrband(1) = 'S' + availrband(2) = 'C' + availrband(3) = 'X' + CALL ErrMsgSAvail("radar band '" // TRIM(model_config_rec%radar_band(ir)) // "' not ready !!",& + 3, availrband, fname) + END IF + + !$OMP PARALLEL DO & + !$OMP PRIVATE ( ij ) + tiles_loop: DO ij = 1 , num_tiles + DO i = i_start(ij), i_end(ij) + DO j = j_start(ij), j_end(ij) + IF (ij == 1 .AND. i == i2(ij) .AND. j == j2(ij)) ijdbg = .TRUE. + i1 = MIN(i+1,i_end(ij)) + i_1 = MAX(i-1,i_start(ij)) + j1 = MIN(j+1,j_end(ij)) + j_1 = MAX(j-1,j_start(ij)) + + IF (ngridradar(i,j,ir) == 0) CYCLE + + ii1 = i1-1 + ei1 = i1+1 + ij1 = j1-1 + ej1 = j1+1 + IF (i == ims) ii1 = i + IF (i == ime) ei1 = i + IF (j == jms) ij1 = j + IF (j == jme) ej1 = j + + ! Decompressing + DO n=1, ngridradar(i,j,ir) + inrnn = (ir-1)*ngrid_radar + n + ijdistgridradar_NOc(n,ir) = ijdistgridradar(i,j,inrnn) + DO l=1, coord2 + inrnncoord = (ir-1)*ngrid_radar*2 + (n-1)*2 + l + ijgridradar_NOc(n,l) = ijgridradar(i,j,inrnncoord) + END DO + END DO + + ! Computing 27 weights around each WRF i,j,k for all radar cells + CALL weights_hcoordsR27around_polar_ellipsoide_1D(i, j, ii1, ei1, ij1, ej1, kms, kme, a, & + f, elipsn, lon(ii1:ei1,ij1:ej1), lat(ii1:ei1,ij1:ej1), zg(ii1:ei1,kms:kme,ij1:ej1), & + model_config_rec%missing_value, ngrid_radar, ngridradar(i,j,ir), coord2, & + ijgridradar_NOc, around2, 1, e_dist, 1, e_angle, 1, e_azimuth, radarlons(:,:,ir), & + radarlats(:,:,ir), radarhgts(:,:,ir), model_config_rec%missing_value, ngrid3_radar, & + ijkngridradar, coord3, ijkgridradar, ijkNtot, around3, wghts, dbg) + + loop_z: DO k=kms, kme + IF (ijkngridradar(k) == 0) CYCLE + + printdbg = .FALSE. + IF (dbg) THEN + IF (MOD(i+j+k,10) == 0 .AND. ijkfound < 10) THEN + IF (ijkfound < 10) THEN + printdbg = .TRUE. + ijkfound = ijkfound + 1 + END IF + END IF + END IF + + IF (printdbg) PRINT *,' ' // TRIM(fname) // ": ", ijkfound, ' i,j,k ', i,j,k, & + ' amount of radar cells:', ijkngridradar(k) + + ! 27 values around reference space at i,j,k + i123 = 1 + cvals = 0. + IF (k==kms) THEN + DO ii1=-1,1 + DO ii2=-1,1 + DO ii3=0,1 + DO iqm=1, Nmoist + iicv = tocompute(iqm) + cvals(i123,iicv) = moist(i+ii1,k+ii3,j+ii2,imoist(iqm)) + END DO + cvals(i123,8) = ua(i+ii1,k+ii3,j+ii2) + cvals(i123,9) = va(i+ii1,k+ii3,j+ii2) + cvals(i123,10) = w(i+ii1,k+ii3,j+ii2) + cvals(i123,11) = ta(i+ii1,k+ii3,j+ii2) + cvals(i123,12) = press(i+ii1,k+ii3,j+ii2) + i123 = i123 + 1 + END DO + END DO + END DO + ELSE IF (k1==kme) THEN + DO ii1=-1,1 + DO ii2=-1,1 + DO ii3=-1,0 + DO iqm=1, Nmoist + iicv = tocompute(iqm) + cvals(i123,iicv) = moist(i+ii1,k+ii3,j+ii2,imoist(iqm)) + END DO + cvals(i123,8) = ua(i+ii1,k+ii3,j+ii2) + cvals(i123,9) = va(i+ii1,k+ii3,j+ii2) + cvals(i123,10) = w(i+ii1,k+ii3,j+ii2) + cvals(i123,11) = ta(i+ii1,k+ii3,j+ii2) + cvals(i123,12) = press(i+ii1,k+ii3,j+ii2) + i123 = i123 + 1 + END DO + END DO + END DO + ELSE + DO ii1=-1,1 + DO ii2=-1,1 + DO ii3=-1,1 + DO iqm=1, Nmoist + iicv = tocompute(iqm) + cvals(i123,iicv) = moist(i+ii1,k+ii3,j+ii2,imoist(iqm)) + END DO + cvals(i123,8) = ua(i+ii1,k+ii3,j+ii2) + cvals(i123,9) = va(i+ii1,k+ii3,j+ii2) + cvals(i123,10) = w(i+ii1,k+ii3,j+ii2) + cvals(i123,11) = ta(i+ii1,k+ii3,j+ii2) + cvals(i123,12) = press(i+ii1,k+ii3,j+ii2) + i123 = i123 + 1 + END DO + END DO + END DO + END IF + + ! Removing too small values... + DO ii1=1, around3 + DO ii2=1, 12 + IF (cvals(ii1,ii2) < eps) cvals(ii1,ii2) = 0. + END DO + END DO + + loop_ptsijk: DO ic= 1, ijkngridradar(k) + !! Angle + iy = ijkgridradar(k,ic,1) + !! Dist + ix = ijkgridradar(k,ic,2) + !! Azimuth + iz = ijkgridradar(k,ic,3) + + IF (printdbg) PRINT *,' cell: ', ic, ' coords; ', ijkgridradar(k,ic,:), & + ' || lon: ', radarlons(iy,ix,ir), ' lat: ', radarlats(iy,ix,ir), ' hgt: ', & + radarhgts(iz,ix,ir), ' amount of values: ', i123 + ! Need to filter Nan and missing values ? + !ijkvals(i,j,k,ic) = 0. + !DO ijk=1,acoords3 + ! ijkvals(i,j,k,ic) = ijkvals(i,j,k,ic) + cvals(ijk)*wghts(i,j,k,ic,ijk) + !END DO + + ! Let's be memory constrained (might be ofered for output checks?) + ! Assuming we will always have, vapour, clouds and rain + iqv = SUM(cvals(:,1)*wghts(k,ic,:)) + + !! Computing radar diagnostics + ! Avoiding masked values + IF (iqv > 0.9*fillValueR) CYCLE + + iqc = SUM(cvals(:,2)*wghts(k,ic,:)) + ! Filtering to 0. + IF (iqc .lt. model_config_rec%q_clou_min) THEN + iqc = 0. + END IF + + iqr = SUM(cvals(:,3)*wghts(k,ic,:)) + IF (iqr .lt. model_config_rec%q_rain_min) THEN + iqr = 0. + END IF + + IF (m_qi .GE. prm_first_scalar) THEN + iqi = SUM(cvals(:,4)*wghts(k,ic,:)) + ! Filtering to 0 + IF (iqi .LT. model_config_rec%q_ice_min) THEN + iqi = 0. + END IF + ELSE + iqi = 0. + END IF + + IF (m_qs .GE. prm_first_scalar) THEN + iqs = SUM(cvals(:,5)*wghts(k,ic,:)) + ! Filtering to 0 + IF (iqs .LT. model_config_rec%q_snow_min) THEN + iqs = 0.d0 + END IF + ELSE + iqs = 0.d0 + END IF + + IF (m_qg .GE. prm_first_scalar) THEN + iqg = SUM(cvals(:,6)*wghts(k,ic,:)) + ! Filtering to 0 + IF (iqg .LT. model_config_rec%q_grau_min) THEN + iqg = 0. + END IF + ELSE + iqg = 0. + END IF + + IF (m_qh .GE. prm_first_scalar) THEN + iqh = SUM(cvals(:,7)*wghts(k,ic,:)) + ! Filtering to 0 + IF (iqh .LT. model_config_rec%q_hail_min) THEN + iqh = 0. + END IF + ELSE + iqh = 0. + END IF + + iua = SUM(cvals(:,8)*wghts(k,ic,:)) + iva = SUM(cvals(:,9)*wghts(k,ic,:)) + iwa = SUM(cvals(:,10)*wghts(k,ic,:)) + ita = SUM(cvals(:,11)*wghts(k,ic,:)) + ipress = SUM(cvals(:,12)*wghts(k,ic,:)) + + IF (printdbg ) PRINT *, ' 27values: ', cvals(:,1), ' 27wghts: ', wghts(k,ic,:), & + ' int. value: ', iqv + + ! Avoiding wasting time + IF (dqr+dqs+dqg+dqh == 0.d0) THEN + zh_rain_NOc = 0. + zdr_rain_NOc = 0. + kdp_rain_NOc = 0. + zh_snow_NOc = 0. + zdr_snow_NOc = 0. + kdp_snow_NOc = 0. + zh_grau_NOc = 0. + zdr_grau_NOc = 0. + kdp_grau_NOc = 0. + zh_ice_NOc = 0. + zdr_ice_NOc = 0. + kdp_ice_NOc = 0. + zh_hail_NOc = 0. + zdr_hail_NOc = 0. + kdp_hail_NOc = 0. + zh_total_NOc = 0. + zdr_total_NOc = 0. + kdp_total_NOc = 0. + tv_rain_NOc = 0. + tv_snow_NOc = 0. + tv_grau_NOc = 0. + wt_NOc = 0. + vr_NOc = 0. + CYCLE + END IF + dua = iua*1.d0 + dva = iva*1.d0 + dwa = iwa*1.d0 + dta = ita*1.d0 + dpress = ipress*1.d0 + + IF (printdbg) THEN + PRINT *,' ' // TRIM(fname) // " radar space: ia,ih,id:", iy,iz,ix + PRINT *,' qv: ', iqv, ' qc: ', iqc, ' qr: ', iqr, ' qci: ', iqi, ' qs: ', & + iqs, ' qg: ', iqg, ' qh: ', iqh + PRINT *,' ta: ',ita, ' press: ', ipress, ' ua: ', iua, ' va: ', iva, ' wa: ', & + iwa + END IF + + ! Polarimetrics_fromLUT_1R(qx, iz, dimv, q_log_vec, Nbnds, Nelevs, Nfrq, Zh_x, Zdr_x, & + ! KDP_x, Zh_out, Zdr_out, KDP_out) + ! Rain + CALL Polarimetrics_fromLUT_1R(iqr, iz, iband, model_config_rec%Nqxbins, & + q_rain_vec, Nlutqxbins, lutZh_RAIN(iband,iz,:), lutZdr_RAIN(iband,iz,:), & + lutKDP_RAIN(iband,iz,:), zh_rain_NOc, zdr_rain_NOc, kdp_rain_NOc, dbg) +! lutKDP_RAIN(iband,iz,:), dZh_q, dZdr_q, dKDP_q, dbg) +! zh_rain_NOc = REAL(dZh_q) +! zdr_rain_NOc = REAL(dZdr_q) +! kdp_rain_NOc = REAL(dKDP_q) + + + ! Snow + IF (m_qs .GE. prm_first_scalar) THEN + CALL Polarimetrics_fromLUT_1R(iqs, iz, iband, model_config_rec%Nqxbins, & + q_snow_vec, Nlutqxbins, lutZh_SNOW(iband,iz,:), lutZdr_SNOW(iband,iz,:), & + lutKDP_SNOW(iband,iz,:), zh_snow_NOc, zdr_snow_NOc, kdp_snow_NOc, dbg) +! lutKDP_SNOW(iband,iz,:), dZh_q, dZdr_q, dKDP_q, dbg) +! zh_snow_NOc = REAL(dZh_q) +! zdr_snow_NOc = REAL(dZdr_q) +! kdp_snow_NOc = REAL(dKDP_q) + END IF + + ! Graupel + IF (m_qg .GE. prm_first_scalar) THEN + CALL Polarimetrics_fromLUT_1R(iqg, iz, iband, model_config_rec%Nqxbins, & + q_grau_vec, Nlutqxbins, lutZh_GRAU(iband,iz,:), lutZdr_GRAU(iband,iz,:), & + lutKDP_GRAU(iband,iz,:), zh_grau_NOc, zdr_grau_NOc, kdp_grau_NOc, dbg) +! lutKDP_RAIN(iband,iz,:), dZh_q, dZdr_q, dKDP_q, dbg) +! zh_grau_NOc = REAL(dZh_q) +! zdr_grau_NOc = REAL(dZdr_q) +! kdp_grau_NOc = REAL(dKDP_q) + END IF + + ! Cloud ice + IF (m_qi .GE. prm_first_scalar) THEN + CALL Polarimetrics_fromLUT_1R(iqi, iz, iband, model_config_rec%Nqxbins, & + q_ice_vec, Nlutqxbins, lutZh_ICE(iband,iz,:), lutZdr_ICE(iband,iz,:), & + lutKDP_ICE(iband,iz,:), zh_ice_NOc, zdr_ice_NOc, kdp_ice_NOc, dbg) +! lutKDP_ICE(iband,iz,:), dZh_q, dZdr_q, dKDP_q, dbg) +! zh_ice_NOc = REAL(dZh_q) +! zdr_ice_NOc = REAL(dZdr_q) +! kdp_ice_NOc = REAL(dKDP_q) + + END IF + + ! Hail + IF (m_qh .GE. prm_first_scalar) THEN + CALL Polarimetrics_fromLUT_1R(iqh, iz, iband, model_config_rec%Nqxbins, & + q_hail_vec, Nlutqxbins, lutZh_HAIL(iband,iz,:), lutZdr_HAIL(iband,iz,:), & + lutKDP_HAIL(iband,iz,:), zh_hail_NOc, zdr_hail_NOc, kdp_hail_NOc, dbg) +! lutKDP_HAIL(iband,iz,:), dZh_q, dZdr_q, dKDP_q, dbg) +! zh_hail_NOc = REAL(dZh_q) +! zdr_hail_NOc = REAL(dZdr_q) +! kdp_hail_NOc = REAL(dKDP_q) + END IF + IF (printdbg) THEN + PRINT *, ' band:', TRIM(model_config_rec%radar_band(ir)), ' RAIN Zh:', & + zh_rain_NOc, ' Zdr: ', zdr_rain_NOc, ' KDP: ', kdp_rain_NOc + END IF + + ! Terminal velocities + CALL calc_tvelocities_1D(iqr, iqs, iqg, iqh, ipress, ita, lutgrid, dtv_rain, & + dtv_snow, dtv_grau, dtv_hail) + tv_rain_NOc = REAL(dtv_rain) + tv_snow_NOc = REAL(dtv_snow) + tv_grau_NOc = REAL(dtv_grau) + tv_hail_NOc = REAL(dtv_hail) + IF (printdbg) PRINT *, ' terminal velocities rain:', tv_rain_NOc, ' snow: ', & + tv_snow_NOc, ' graupel: ', tv_grau_NOc, ' hail: ', tv_hail_NOc + + ! Calcualte radial velocities ... (not full spectrum) + ! The average radial velocity vrad is the power-weighted sum of the projections of + ! U (eastward wind component), V (northward wind component), W (vertical wind component), + ! and vt, the hydrometeor terminal velocity, onto the axis of the radar beam defined by + ! angle and azimuth + dwt = 0.D0 + dvr = 0.D0 + IF (model_config_rec%doppler_scheme == 0) THEN + !doppler_scheme_name = 'massweightedmeanvt' + + ! Reflectivity weighted terminal velocity + ! 'Same approximation as wrf2radar' + ! Radial velocity (calculated geometrically) + dangle = radars_angle(iy,ir)*1.D0 / dRad + dazimuth = radars_azimuth(iz,ir)*1.D0 / dRad + cos_el = DCOS(dazimuth) + sin_el = DSIN(dazimuth) + sin_an = DSIN(dangle) + cos_an = DCOS(dangle) + denom = zh_rain_NOc + zh_snow_NOc + zh_grau_NOc + zh_hail_NOc + IF (denom /= 0.D0) THEN + dwt = (dtv_rain*zh_rain_NOc + dtv_snow*zh_snow_NOc + dtv_grau*zh_grau_NOc + & + dtv_hail*zh_hail_NOc) / denom + wt_NOc = REAL(dwt) + END IF + + ! Radial velocity + dvr = dua*cos_el*sin_an + dva*cos_el*cos_an + (dwa-wt_NOc)*sin_el + vr_NOc = REAL(dvr) + + IF (printdbg) THEN + PRINT *, ' wt:', wt_NOc + PRINT *, ' radial terminal velocity:', vr_NOc + END IF + + ELSE IF (model_config_rec%doppler_scheme == 1) THEN + ! Get PSD integrated backscattering cross section int = v S/S + CALL ErrMsg('Still need to add this functionality. Get the backscattering cross '// & + 'section for everygrid point', fname, -1) + END IF + + ! Compressing + inrd = (ir-1)*e_dist + ix + + !! Optional + !! SIMRADARop !! + !!iqvm(iy,iz,inrd) = iqv + !!iqcm(iy,iz,inrd) = iqc + !!iqrm(iy,iz,inrd) = iqr + !!iqim(iy,iz,inrd) = iqi + !!iqsm(iy,iz,inrd) = iqs + !!iqgm(iy,iz,inrd) = iqg + !!iqhm(iy,iz,inrd) = iqh + !!iuam(iy,iz,inrd) = iua + !!ivam(iy,iz,inrd) = iva + !!iwam(iy,iz,inrd) = iwa + !!itam(iy,iz,inrd) = ita + !!ipressm(iy,iz,inrd) = ipress + + zh_rain(iy,iz,inrd) = zh_rain_NOc + zdr_rain(iy,iz,inrd) = zdr_rain_NOc + kdp_rain(iy,iz,inrd) = kdp_rain_NOc + + zh_snow(iy,iz,inrd) = zh_snow_NOc + zdr_snow(iy,iz,inrd) = zdr_snow_NOc + kdp_snow(iy,iz,inrd) = kdp_snow_NOc + + zh_grau(iy,iz,inrd) = zh_grau_NOc + zdr_grau(iy,iz,inrd) = zdr_grau_NOc + kdp_grau(iy,iz,inrd) = kdp_grau_NOc + + zh_ice(iy,iz,inrd) = zh_ice_NOc + zdr_ice(iy,iz,inrd) = zdr_ice_NOc + kdp_ice(iy,iz,inrd) = kdp_ice_NOc + + zh_hail(iy,iz,inrd) = zh_hail_NOc + zdr_hail(iy,iz,inrd) = zdr_hail_NOc + kdp_hail(iy,iz,inrd) = kdp_hail_NOc + + zh_total(iy,iz,inrd) = zh_rain_NOc + zh_snow_NOc + zh_grau_NOc + zh_ice_NOc + & + zh_hail_NOc + zdr_total(iy,iz,inrd) = zdr_rain_NOc + zdr_snow_NOc + zdr_grau_NOc + zdr_ice_NOc + & + zdr_hail_NOc + kdp_total(iy,iz,inrd) = kdp_rain_NOc + kdp_snow_NOc + kdp_grau_NOc + kdp_ice_NOc + & + kdp_hail_NOc + + tv_rain(iy,iz,inrd) = tv_rain_NOc + tv_snow(iy,iz,inrd) = tv_snow_NOc + tv_grau(iy,iz,inrd) = tv_grau_NOc + tv_ice(iy,iz,inrd) = tv_ice_NOc + tv_hail(iy,iz,inrd) = tv_hail_NOc + + tv_wgt(iy,iz,inrd) = wt_NOc + vr(iy,iz,inrd) = vr_NOc + + END DO loop_ptsijk + END DO loop_z + + IF (dbg_level >= 75 .AND. ijdbg .AND. tdbg==0) THEN + WRITE(msg,*) ' ' // TRIM(fname) // ' ',i,' ,',j,' : final checks _______' + CALL wrf_debug(75,msg) + END IF + + ijdbg = .FALSE. + + END DO + END DO + END DO tiles_loop + !$OMP END PARALLEL DO + + END DO radars_loop + + RETURN + + END SUBROUTINE simradar_calc + +#endif +END MODULE module_simradar diff --git a/phys/module_simradar_tools.F b/phys/module_simradar_tools.F new file mode 100644 index 0000000000..9b108fd33d --- /dev/null +++ b/phys/module_simradar_tools.F @@ -0,0 +1,3303 @@ +MODULE module_simradar_tools +#ifdef SIMRADAR +! Module wiuth the tools for simrad.AR +! L. Fita, CIMA-IFAECI, November 2025 + + USE module_wrf_error + +!!!!!!! Subroutines / functions +! beam_height: Function to compute the height of the beam +! calc_tvelocities_1D: Subroutine to compute terminal velocities for water species 1D version +! calc_vt: Function to compute terminal velocity after a gamma function +! com_gamma: Gamma function +! compute_radar_lonlat_radial: Subroutine to copmute lon, lat, hgt coordinates in the radar matrix +! following polar coordinates +! create_logarithmic_scaled_vectorsD: Function to create logarithmic-based equi-spaced scaled double +! Real vectors +! distance_between_points: Subroutine tro compute the coordinates at a given distance and angle from a +! given point +! DIRCT1: Solution of the geodetic direct problem after T. Vincenty +! inverse: Subroutine to compute distance between 2 points +! is_inside: Subroutine to check if a tile gets a radar +! lin_interpolate_1RR: Function to provide linear interpolation a Real location on a given Real axis +! pertinence_polar: Subroutine to compute the 2-dimensional pertinence of a lonlat polar radial space +! (r) into another regular lonlat (WRF) +! Polarimetrics_fromLUT_1R: Subroutine to interpolate the polarimetic characteristics of a given water +! specie following T-matrix +! position_from_point: Subroutine tro compute the coordinates at a given distance and angle from a +! given point +! read_lut: Subroutine to read the T-MATRIX look-up-table (LUT) from an specific microphysics from a +! netCDF file +! radar_inside: Subroutine to determine if a radar fells inside the tile +! real_distances: Subroutine to compute the diferential of coordinattes of a distance from a given +! point +! real_distances_ellipsoide: Subroutine to compute the diferential of coordinattes of a distance from +! a given point on an ellipsoide +! weights_hcoordsR27around_polar_ellipsoide_1D: Subroutine to compute weights of the 27 around grid +! points to 3-dimensional interpolate to a new space using the list of cells on the horizontal of +! the targeted space that lay within each grid cell of the reference space with the values to +! interpolate version on an ellipsoide + +! basic: (from PyNCplot) +! ErrMsg: Subroutine to stop execution and provide an error message +! ErrMsg[I/R/S]Avail: Subroutine to stop execution by non exsitent Integer/Real/String value from +! available ones and provide an error message +! FileExist: Checks if a file exists +! freeunit: provides the number of a free unit in which open a file +! [I/R/D]toS: Function to transform an Integer/Real/double Real to String +! removeNONnum: Subroutine to remove non numeric characters from a string +! set_dbg: Function to setup debug value +! Sto[I/R/L]: Function to transform a String to an Integer/Real/Boolean + +! netCDF: (from PyNCplot) +! close_NCfile: Subroutine to close a netCDF file +! get_dimlength_ncid: Subroutine to get the length of a dimension from an existing to an open netCDF +! get_dims_ncid: Subroutine to get the dimensions from an open netCDF file +! get_Ndims_ncid: Subroutine to get the amount of dimensions from an open netCDF file +! get_Nvars_ncid: Subroutine to get the amount of variables from an open netCDF file +! get_vars_ncid: Subroutine to provide all the ids and names of the variables existing to an open +! netCDF +! get_varR[0/1]D_ncid: Subroutine to get a [0/1]D Real variable from a netCDF file unit passing +! starting and ending +! get_varD3D_ncid: Subroutine to get a 3D double Real variable from a netCDF file unit passing +! starting and ending +! get_varS0D_ncid: Subroutine to get a 0D String variable from a netCDF file unit passing starting +! and ending +! handle_errf: Subroutine to provide the error message when something with netCDF went wrong +! (including fname) +! isin_ncunit: Function to tell if a given variable is inside a netcdf file unit +! isin_dim_ncid: Function to tell if a given dimension is inside a netcdf file unit +! open_NCfile: Subroutine to open a basic netCDF file + + IMPLICIT NONE + + ! Integer 4 bytes + INTEGER, PARAMETER :: ii4 = SELECTED_INT_KIND(4) + ! REal 16 bytes + INTEGER, PARAMETER :: rr8 = SELECTED_REAL_KIND(8) + + REAL, PARAMETER :: pict = ACOS(-1.) + REAL, PARAMETER :: DegRad = 180./pict + ! double Real + REAL(kind=rr8), PARAMETER :: dpict = 4.d0*DATAN(1.d0) + REAL(kind=rr8), PARAMETER :: drad = 180.d0/dpict + + ! 16 charcater String + INTEGER, PARAMETER :: S16 = 16 + ! 32 character String + INTEGER, PARAMETER :: S32 = 32 + ! 50 character String + INTEGER, PARAMETER :: Sl = 50 + ! 64 character String + INTEGER, PARAMETER :: S64 = 64 + ! 128 character String + INTEGER, PARAMETER :: S128 = 128 + ! 256 character String + INTEGER, PARAMETER :: S256 = 256 + ! 512 character String + INTEGER, PARAMETER :: S512 = 512 + ! 1024 character String + INTEGER, PARAMETER :: S1024 = 1024 + + ! Fill values + ! !! 'B', NF_BYTE + ! BYTE, PARAMETER :: fillValueB = 255 + !! 'C', NF_CHAR + CHARACTER(len=1), PARAMETER :: fillValueC = '-' + !! 'I', NF_SHORT + INTEGER, PARAMETER :: fillValueI = -9999 + !!! 'I16', NF_INT + !INTEGER(kind=ii4), PARAMETER :: fillValueID = fI4 + !! 'R', NF_FLOAT + REAL, PARAMETER :: fillValueR = 1.e20 + !! 'D', NF_DOUBLE + REAL(kind=rr8), PARAMETER :: fillValueD = 1.d20 + + CHARACTER(len=33) :: errormsg = 'ERROR -- error -- ERROR -- error' + + ! error tolerance + REAL(kind=rr8), PARAMETER :: tol=1.d-14 + ! Epsilon + REAL(kind=rr8), PARAMETER :: eps=1.d-15 + + ! Maximum rank of a Frotran variable + INTEGER, PARAMETER :: maxfrank = 7 + + ! simrad.AR version + CHARACTER(len=S32) :: simradARversion = '1.0' + + ! 4 / 3 radar beam height + REAL(kind=rr8), PARAMETER :: ke = 4.d0 / 3.d0 + + ! Derived type with all the ranges for all water species for LUT + TYPE lut_vars + INTEGER :: Nbands + INTEGER :: Nlutqxbins + INTEGER :: Nlutelevs + CHARACTER(len=S16) :: mp_scheme + + REAL, DIMENSION(:), ALLOCATABLE :: elevations + + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zh_RAIN + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zdr_RAIN + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: LDR_RAIN + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aih_RAIN + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aiv_RAIN + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: KDP_RAIN + + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zh_SNOW + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zdr_SNOW + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: LDR_SNOW + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aih_SNOW + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aiv_SNOW + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: KDP_SNOW + + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zh_GRAU + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zdr_GRAU + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: LDR_GRAU + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aih_GRAU + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aiv_GRAU + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: KDP_GRAU + + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zh_ICE + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zdr_ICE + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: LDR_ICE + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aih_ICE + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aiv_ICE + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: KDP_ICE + + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zh_HAIL + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Zdr_HAIL + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: LDR_HAIL + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aih_HAIL + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: Aiv_HAIL + REAL(kind=rr8), DIMENSION(:,:,:), ALLOCATABLE :: KDP_HAIL + + ! Which variables need to be interpolated / used depending on the microphysics scheme + INTEGER :: N_tocompute + INTEGER, DIMENSION(12) :: tocompute + !! Indices of the water species to compute / use + INTEGER :: Nmoist + INTEGER, DIMENSION(12) :: imoist + + REAL :: q_rain_min + REAL :: q_rain_max + REAL :: qn_rain_min + REAL :: qn_rain_max + REAL :: q_snow_min + REAL :: q_snow_max + REAL :: qn_snow_min + REAL :: qn_snow_max + REAL :: q_grau_min + REAL :: q_grau_max + REAL :: qn_grau_min + REAL :: qn_grau_max + REAL :: q_ice_min + REAL :: q_ice_max + REAL :: qn_ice_min + REAL :: qn_ice_max + REAL :: q_hail_min + REAL :: q_hail_max + REAL :: qn_hail_min + REAL :: qn_hail_max + REAL :: q_clou_min + REAL :: q_clou_max + INTEGER :: dim_ice + INTEGER :: dimv + + REAL(kind=rr8), DIMENSION(:), ALLOCATABLE :: q_rain_vec + REAL(kind=rr8), DIMENSION(:), ALLOCATABLE :: q_snow_vec + REAL(kind=rr8), DIMENSION(:), ALLOCATABLE :: q_grau_vec + REAL(kind=rr8), DIMENSION(:), ALLOCATABLE :: q_ice_vec + REAL(kind=rr8), DIMENSION(:), ALLOCATABLE :: q_hail_vec + + ! Denisities + REAL :: rhoo + REAL :: rho_r + REAL :: n0r + REAL :: rho_s + REAL :: n0s + REAL :: rho_g + REAL :: n0g + REAL :: rho_h + REAL :: n0h + + END TYPE lut_vars + + CONTAINS + + SUBROUTINE create_logarithmic_scaled_vectorsD(vec_min, vec_max, Nvals, logvec, idbg) + ! Function to create logarithmic-based equi-spaced scaled double Real vectors + + IMPLICIT NONE + + INTEGER, INTENT(in) :: Nvals + REAL, INTENT(in) :: vec_min, vec_max + LOGICAL, OPTIONAL :: idbg + REAL(kind=rr8), DIMENSION(NVals), INTENT(out) :: logvec + + ! Local + INTEGER :: i + REAL(kind=rr8) :: log10n, log10x, ddvec + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! vec_min: minimum value + ! vec_max: maximum value + ! Nvals: amount of values + fname = 'create_logarithmic_scaled_vectorsD' + + dbg = set_dbg(idbg) + + log10n = DLOG10(vec_min*1.D0) + log10x = DLOG10(vec_max*1.D0) + + ddvec = (log10x - log10n)/(Nvals-1)*1.D0 + + DO i=1, Nvals + logvec(i) = (10.D0)**(log10n + (i-1)*ddvec) + END DO + + IF (dbg) THEN + PRINT *, ' ' // TRIM(fname) // ": log10 min: ", log10n, ' max: ', log10x, & + ' estimated distance: ', ddvec + PRINT *, ' LOG10 distance between consecutive log10(1), log10(2), dist:', & + DLOG10(logvec(1)), ', ', DLOG10(logvec(2)), ', ', DLOG10(logvec(2)) - DLOG10(logvec(1)) + PRINT *, ' LOG10 distance between consecutive log10(Ndim-1), log10(Ndim), dist:', & + DLOG10(logvec(Nvals-1)), ', ', DLOG10(logvec(Nvals)), ', ', DLOG10(logvec(Nvals)) - & + DLOG10(logvec(Nvals-1)) + END IF + + RETURN + + END SUBROUTINE create_logarithmic_scaled_vectorsD + + REAL(kind=rr8) FUNCTION lin_interpolate_1RR(pt, Ndim, coords, vals) +! If willing to manage error as function +! REAL(kind=rr8) FUNCTION lin_interpolate_1RD(pt, Ndim, coords, vals, statv) + ! Function to provide linear interpolation a Real location on a given real axis + + IMPLICIT NONE + + INTEGER, INTENT(in) :: Ndim + REAL, INTENT(in) :: pt + REAL, DIMENSION(Ndim), INTENT(in) :: coords, vals + ! 0: no error; 1: error +! INTEGER, INTENT(out) :: statv + + ! Local + INTEGER :: icor + INTEGER :: refi + REAL :: minc, maxc + REAL :: signv, a, b + REAL, DIMENSION(Ndim) :: origsing + CHARACTER (LEN=1000) :: diag_message + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! pt: pt to interpolate + ! coords: coordinates along the axis (incremental) + ! vals: values at the coordinates + fname = 'lin_interpolate_1RR' + + lin_interpolate_1RR = fillValueD + ! Sign of the axis + origsing = ABS(coords(Ndim)-coords(1))/(coords(Ndim)-coords(1)) + + minc = MINVAL(coords) + maxc = MAXVAL(coords) + + a = fillValueR + refi = -1 + IF (pt < minc) THEN + a = vals(1) + refi = 1 + signv = -1. + ELSE IF (pt > maxc) THEN + a = vals(Ndim) + refi = Ndim-1 + signv = 1. + ELSE IF (pt == minc) THEN + lin_interpolate_1RR = vals(1) + refi = 1 + signv = 1. + + RETURN + ELSE IF (pt == maxc) THEN + lin_interpolate_1RR = vals(Ndim) + refi = Ndim-1 + signv = 1. + + RETURN + ELSE + DO icor=1, Ndim-1 + IF (pt == coords(icor)) THEN + lin_interpolate_1RR = vals(icor) + RETURN + ELSE IF (pt >= coords(icor) .AND. pt < coords(icor+1)) THEN + a = vals(icor) + refi = icor + signv = 1. + EXIT + END IF + END DO + END IF + +! statv = 0 + IF (refi == -1) THEN + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " //TRIM(fname)// ": No interpolation found for: " // TRIM(RtoS(pt)) // & + ' within range: [' // TRIM(RtoS(minc)) // ', ' // TRIM(RtoS(maxc)) + CALL wrf_error_fatal ( diag_message ) + END IF + + b = ABS(pt-coords(refi))/(coords(refi+1)-coords(refi)) + + lin_interpolate_1RR = a + signv*(vals(refi+1)-a)*b + + RETURN + + END FUNCTION lin_interpolate_1RR + + SUBROUTINE com_gamma(x, ga) + ! GAMMA FUNCTION + !----------------------------------------------------------------------- + + !================================================== + ! Purpose: Compute the gamma function â(x) + ! Input : x --- Argument of â(x) + ! ( x is not equal to 0,-1,-2,úúú ) + ! Output: GA --- â(x) + !================================================== + ! Proff. Jianming Jin + ! Department of Electrical and Computer Engineering + ! University of Illinois + + IMPLICIT NONE + + REAL(kind=rr8), INTENT(in) :: X + REAL(kind=rr8), INTENT(out) :: GA + + ! Local + REAL(kind=rr8), DIMENSION(26) :: G + REAL(kind=rr8) :: Z , R , GR + INTEGER :: M1, K , M + CHARACTER(len=Sl) :: fname + + DATA G /1.0D0, 0.5772156649015329D0, -0.6558780715202538D0, -0.420026350340952D-1, & + 0.1665386113822915D0,-.421977345555443D-1, -.96219715278770D-2, .72189432466630D-2, & + -.11651675918591D-2, -.2152416741149D-3, .1280502823882D-3, -.201348547807D-4, & + -.12504934821D-5, .11330272320D-5, -.2056338417D-6, .61160950D-8, .50020075D-8, -.11812746D-8, & + .1043427D-9, .77823D-11, -.36968D-11, .51D-12, -.206D-13, -.54D-14, .14D-14, .1D-15 / + + !!!!!!! Variables + ! x: argument + fname = 'com_gamma' + + IF (X.EQ.INT(X)) THEN + IF (X.GT.0.0D0) THEN + GA=1.0D0 + M1=X-1 + DO K=2,M1 + GA=GA*K + END DO + ELSE + GA=1.0D+300 + ENDIF + + ELSE + IF (DABS(X).GT.1.0D0) THEN + Z = DABS(X) + M = INT(Z) + R = 1.0D0 + DO K=1, M + R = R*(Z-K) + END DO + Z=Z-M + ELSE + Z=X + ENDIF + GR=G(26) + DO K = 25, 1, -1 + GR = GR*Z + G(K) + END DO + + GA=1.0D0/(GR*Z) + IF (DABS(X).GT.1.0D0) THEN + GA=GA*R + IF (X.LT.0.0D0) GA=-DPICT/(X*GA*DSIN(DPICT*X)) + ENDIF + ENDIF + + RETURN + + END SUBROUTINE com_gamma + + SUBROUTINE Polarimetrics_fromLUT_1R(qx, iz, ibnd, dimv, q_log_vec, Nbins, Zh_x, Zdr_x, KDP_x, & + Zh_out, Zdr_out, KDP_out, idbg) + ! Subroutine to interpolate the polarimetic characteristics of a given water specie following T-matrix + + IMPLICIT NONE + + INTEGER, INTENT(in) :: iz, ibnd, dimv, Nbins + REAL, INTENT(in) :: qx + REAL, DIMENSION(dimv), INTENT(in) :: q_log_vec + REAL, DIMENSION(Nbins), INTENT(in) :: Zh_x, Zdr_x, KDP_x + LOGICAL, OPTIONAL :: idbg + REAL, INTENT(out) :: Zh_out, Zdr_out, KDP_out + + ! Local + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! qx: value of the water specie + ! iz: lut azimuth index + ! ibnd: band + ! 1: 'C' + ! 2: 'S' + ! 3: 'X' + ! dimv: length of the logarithmic scale + ! q_log_vec: logarithmic scale + ! Nbins: amount of bins + ! Zh_x: values of 'Zh' for the specific water specie + ! Zdr_x: values of 'Zdr' for the specific water specie + ! KDP_x: values of 'KDP' for the specific water specie + fname = 'Polarimetrics_fromLUT_1R' + + IF (dbg) THEN + PRINT *, ' ' // TRIM(fname) // ": q_log_vec 1 a 5:", q_log_vec(1:5) + PRINT *, ' ' // TRIM(fname) // ": Zh_x at ibnd= ", ibnd ," iz= ", iz, " 1 a 5:", Zh_x(1:5) + END IF + + dbg = set_dbg(idbg) + + Zh_out = fillValueR + Zdr_out = fillValueR + KDP_out = fillValueR + + Zh_out = lin_interpolate_1RR(qx, dimv, q_log_vec, Zh_x) + Zdr_out = lin_interpolate_1RR(qx, dimv, q_log_vec, Zdr_x) + KDP_out = lin_interpolate_1RR(qx, dimv, q_log_vec, KDP_x) + + RETURN + + END SUBROUTINE Polarimetrics_fromLUT_1R + + REAL(kind=rr8) FUNCTION calc_vt(a, b, lambdax) + ! Function to compute terminal velocity after a gamma function + + IMPLICIT NONE + + REAL, INTENT(in) :: a, b + REAL(kind=rr8), INTENT(in) :: lambdax + + ! Local + REAL(kind=rr8) :: db, gammav, lambda6 + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! a, b: gamma parameters + ! lambdax: lambda of the water specie + + fname = 'calc_vt' + + calc_vt = fillValueD + + db = b*1.D0 + CALL com_gamma(4.D0+db, gammav) + + calc_vt = a*gammav/6 + lambda6 = lambdax**db + + calc_vt = calc_vt / lambda6 + + RETURN + + END FUNCTION calc_vt + + SUBROUTINE calc_tvelocities_1D(qr, qs, qg, qh, p, T, lutgrid, v_rain, v_snow, v_grau, v_hail) + ! Subroutine to compute terminal velocities for water species 1D version + ! + ! Calculate the MASS-WEIGHTED MEAN terminal velocities following Srivastava (1967). + ! int vt l(D) dD/l, where vt is the terminal velocity of diamater D, l(D) is the mixing + ! ratio of the precipitation particle of diameter D, and l is the mixing ratio of a precip. field + ! see eq. 11-13 in Lin et al. 1983. (same as wrf2radar) + + IMPLICIT NONE + + REAL, INTENT(in) :: qr, qs, qg, qh, p, T + TYPE(lut_vars), INTENT(in) :: lutgrid + REAL(kind=rr8), INTENT(out) :: v_rain, v_snow, v_grau, v_hail + + ! Local + REAL :: rhoo, rho_air, rhofac + REAL :: rho_r, n0r, rho_s, n0s, rho_g, n0g, & + rho_h, n0h + REAL(kind=rr8) :: lambdar, lambdas, lambdag, lambdah + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! qr: mixing ratio for rain [kgkg-1] + ! qs: mixing ratio for snow [kgkg-1] + ! qg: mixing ratio for graupel [kgkg-1] + ! qh: mixing ratio for hail [kgkg-1] + ! p: pressure [Pa] + ! T: air temperature [K] + ! lutgrid: derived type with all the T-MATRIX look-up-table values + fname = 'calc_tvelocities_1D' + + !! Attention bring herer rl rhofactor! + !!!! THIS must be generalized from lutgrid !!!! + !!! reference air density used in some WSM6 docs, kg m^-3 + rhoo = lutgrid%rhoo + !!! kg/m3 + rho_air = p / (287.058 * T) + rhofac = (rhoo/rho_air)**0.5 + + rho_r = lutgrid%rho_r + n0r = lutgrid%n0r + rho_s = lutgrid%rho_s + n0s = lutgrid%n0s + rho_g = lutgrid%rho_g + n0g = lutgrid%n0g + rho_h = lutgrid%rho_h + n0h = lutgrid%n0g + + v_rain = fillValueD + v_snow = fillValueD + v_grau = fillValueD + + lambdar = fillValueR + lambdas = fillValueR + lambdag = fillValueR + + ! Rain + IF (qr > 0.) THEN + lambdar = (dpict*rho_r*n0r / (qr))**(1./4.) + v_rain = calc_vt(841.9, 0.8, lambdar)*rhofac + END IF + + ! Snow + IF (qs > 0.) THEN + lambdas = (dpict*rho_s*n0s / (qs))**(1./4.) + v_snow = calc_vt(11.72, 0.41, lambdas)*rhofac + END IF + + ! Grau + IF (qg > 0.) THEN + lambdag = (dpict*rho_g*n0g / (qg))**(1./4.) + v_grau = calc_vt(330., 0.8, lambdag)*rhofac + END IF + + ! Hail (LLUIS not know range !!) + IF (qh > 0.) THEN + lambdah = (dpict*rho_h*n0h / (qh))**(1./4.) + v_hail = calc_vt(-9999., -9999., lambdah)*rhofac + END IF + + RETURN + + END SUBROUTINE calc_tvelocities_1D + + SUBROUTINE weights_hcoordsR27around_polar_ellipsoide_1D(i1, j1, id1x, ed1x, id1y, ed1y, id1z, & + ed1z, a, f, elips, lon1, lat1, vert1, miss_v1, maxNcoords, ij21Ncoords, Nc2, ij21coords, & + Na2, id2x, ed2x, id2y, ed2y, id2z, ed2z, lon2, lat2, vert2, miss_v2, maxN3coords, ijk21Ncoords, & + Nc3, ijk21coords, totngridradar, Na3, wghts, idbg) + ! Subroutine to compute weights of the 27 around grid points to 3-dimensional interpolate to a new + ! space using the list of cells on the horizontal of the targeted space that lay within each grid + ! cell of the reference space with the values to interpolate version on an ellipsoide + ! lon1,lat1,vert1: space with data to interpolate (reference) + ! lon2,lat2,vert2: polar space towards which interpolate (target) + ! id1x= i1-1, edx1= i1+1; id1y= j1-1, edy1= j1+1; id1z=1, ed1z=dimz + ! NOTE in polar space: + ! dx: distance dy: angle dz: azimuth + ! 3Dvariables (dy,dz,dx), lon,lats(dy,dx), hgt(dx,dz) + + IMPLICIT NONE + + INTEGER, INTENT(in) :: i1, j1, id1x, ed1x, id1y, ed1y, id1z, & + ed1z, MaxNcoords, MaxN3coords, Nc2, Na2, Nc3, Na3 + INTEGER, INTENT(in) :: id2x, ed2x, id2y, ed2y, id2z, ed2z + REAL, INTENT(in) :: miss_v1, miss_v2 + REAL, INTENT(in) :: a, f + CHARACTER(len=30), INTENT(in) :: elips + REAL, DIMENSION(id1x:ed1x,id1y:ed1y), INTENT(in) :: lon1, lat1 + REAL, DIMENSION(id1x:ed1x,id1y:ed1y,id1z:ed1z), & + INTENT(in) :: vert1 + INTEGER, INTENT(in) :: ij21Ncoords + INTEGER, DIMENSION(MaxNcoords,Nc2), INTENT(in) :: ij21coords + REAL, DIMENSION(id2y:ed2y,id2x:ed2x), INTENT(in) :: lon2, lat2 + REAL, DIMENSION(id2z:ed2z,id2x:ed2x), INTENT(in) :: vert2 + LOGICAL, OPTIONAL :: idbg + INTEGER, INTENT(out) :: totngridradar + INTEGER, DIMENSION(id1z:ed1z), INTENT(out) :: ijk21Ncoords + REAL, DIMENSION(id1z:ed1z,maxN3coords,Na3), & + INTENT(out) :: wghts + INTEGER, DIMENSION(id1z:ed1z,maxN3coords,Nc3), & + INTENT(out) :: ijk21coords + + ! Local + INTEGER :: k1, ic, i2, j2, k2, kk2, ij, ijk + INTEGER :: ici, icj, ick + INTEGER :: ijkfound + REAL :: ijkvert1, sumdists, totwght + REAL(kind=rr8) :: da, df + INTEGER, DIMENSION(1) :: kmindist + INTEGER, DIMENSION(1) :: imaxwght, izero + REAL, DIMENSION(Na2) :: ijlon1, ijlat1, ijdists + REAL, DIMENSION(Na3) :: zdists, dists, wgs + REAL, DIMENSION(id1z:ed1z) :: vert1mask, zdist1 + LOGICAL, DIMENSION(Na3) :: ijkmask + LOGICAL :: printdbg + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! id1x, ed1x, id1y, ed1y, id1z, ed1z: Extent of the reference space + ! a: semimajor axis equatorial (in meters) + ! f: flattening + ! elips: name of the ellipsis being used + ! lon1, lat1, vert1: coordinates of the reference space [degrees_east, degrees_north, m] + ! MaxNcoords: maximum amount of target space grid points in a single reference grid point + ! ij21Ncoords: amount of target grid points in a single reference grid pointv + ! Nc2: amount of coordinates of 2-dimensional space + ! ij21coords: coordinates of the target grid points in a single reference grid point + ! ij21dists: horizonatal distancess of the target grid points in a single reference grid point + ! id2x, ed2x, id2y, ed2y, id2z, ed2z: Extent of the target space + ! lon2, lat2, vert2: coordinates of the target space [degrees_east, degrees_north, m] + ! ijk21Ncoords: amount of target grid 3-dimensional points in a single reference grid pointv + ! Nc3: amount of coordinates of 3-dimensional space + ! ijk21coords: 3-dimensional coordinates of the target grid points in a single reference grid point + ! Na3: amount of points around in a 3-dimensional space + ! wghts: weights to use to interpolate values of a target grid point to its closest 3x3x3 grid + ! points in the reference space + ! Rdists: distances of a target grid point to its closest 3x3x3 grid points in the reference space + fname = 'weights_hcoordsR27around_polar_ellipsoide_1D' + + dbg = set_dbg(idbg) + + ijk21Ncoords = 0 + ijk21coords = fillValueI + wghts = 0. + + da = a*1.D0 + df = f*1.D0 + + IF (dbg) PRINT *,' ' // TRIM(fname) // ": Compute 3-dimensional weights at ", i1, ', ', j1 , & + ' with ', ij21Ncoords, " radar cells horizontally" + + totngridradar = 0 + ijkfound = 0 + + gridpts: DO ic=1, ij21Ncoords + ! Indices of target space within reference grid + !! angle + ici = ij21coords(ic,1) + !! dist + icj = ij21coords(ic,2) + + ! computing distances from 3x3 horizontal grid (the already computed is only for i1,j1) + ! horizontal 3x3 array as + ! -1,-1,(1,1)[1]: SW; 0,-1,(2,1)[2]: SC, 1,-1,(3,1)[3]: SE; + ! -1,0,(1,2)[4]: CW; 0,0,(2,2)[5]: CC; 1,0,(3,2)[6]: CE; + ! -1,1,(1,3)[7]: NW; 0,1,(2,3)[8]: NC; 1,1,(3,3)[9]: NE; + ij = 1 + DO i2=-1,1 + DO j2=-1,1 + ijlon1(ij) = lon1(i1+i2,j1+j2) + ijlat1(ij) = lat1(i1+i2,j1+j2) + ij = ij + 1 + END DO + END DO + + CALL real_distances(1, 3, 1, 3, ijlon1, ijlat1, lon2(ici,icj), lat2(ici,icj), ijdists) + + ! Computing vertical distances (from closest i,j) + zloop: DO k2=id2z, ed2z + IF (vert2(k2,icj) /= miss_v2) THEN + + printdbg = .FALSE. + IF (dbg) THEN + IF (MOD(ic+k2, 10) == 0 .AND. ijkfound < 10) THEN + printdbg = .TRUE. + ijkfound = ijkfound + 1 + END IF + END IF + + ! vertical level at center of 3x3 + DO k1=id1z, ed1z + IF (vert1(i1,j1,k1) == miss_v1) THEN + vert1mask(k1) = fillValueR + ELSE + vert1mask(k1) = vert1(i1,j1,k1) + END IF + END DO + + zdist1 = ABS( vert1mask - vert2(k2,icj) ) + kmindist = MINLOC(zdist1) + ick = kmindist(1) + + IF (printdbg) PRINT *,' ' // TRIM(fname) // ": i1,j1:", i1, j1, " i2, j2, k2: ", & + ici, icj, k2, ' vert1 masked:', vert1mask, ' vert2: ', vert2(k2,icj), ' zdist1: ', & + zdist1, ' ick:', ick, ' found: ', totngridradar+1 + + ijk21Ncoords(ick) = ijk21Ncoords(ick) + 1 + IF (ijk21Ncoords(ick) > maxN3coords) THEN + CALL ErrMsg("Too much radar cells: " // TRIM(ItoS(ijk21Ncoords(ick))) // & + " in a WRF grid. Increase 'ngrid_radar3' in namelist current value: " // & + TRIM(ItoS(maxN3coords)), fname, -1) + END IF + ijk21coords(ick,ijk21Ncoords(ick),1) = ici + ijk21coords(ick,ijk21Ncoords(ick),2) = icj + ijk21coords(ick,ijk21Ncoords(ick),3) = k2 + totngridradar = totngridradar + 1 + + ! vertical 27 array as + ! SW + ! -1,-1,-1,(1,1,1),[1]: bottom; -1,-1,0,(1,1,2),[2]: middle, -1,-1,1,(1,1,3),[3]: top; + ! (...) + ! NE + ! 1,1,-1,(3,3,1),[25]: bottom; 1,1,0,(3,3,2),[26]: middle, 1,1,1,(3,3,3),[27]: top; + ij = 1 + ijk = 1 + dists = fillValueR + ijkmask = .FALSE. + IF (ick == 1) THEN + IF (printdbg) PRINT *, ' ' // TRIM(fname) // ": radar height ", vert2(k2,icj), & + " at first WRF level ", ick, ' hgt= ', vert1(i1,j1,ick), ' m' + DO i2=-1,1 + DO j2=-1,1 + DO kk2=0,1 + ijkvert1 = vert1(i1+i2,j1+j2,ick+kk2) + IF (ijkvert1 /= miss_v1) THEN + ! 3-dimensional distance + ! dist = SQRT(hdist**2 + zdist**2) + dists(ijk) = SQRT(ijdists(ij)**2 + (ijkvert1-vert2(k2,icj))**2) + ijkmask(ijk) = .TRUE. + END IF + ijk = ijk + 1 + END DO + ij = ij + 1 + END DO + END DO + ELSE IF (ick == ed1z) THEN + IF (printdbg) PRINT *, ' ' // TRIM(fname) // ": radar height ", vert2(k2,icj), & + " at last WRF level ", ed1z, ' hgt= ', vert1(i1,j1,ed1z), ' m' + DO i2=-1,1 + DO j2=-1,1 + DO kk2=-1,0 + ijkvert1 = vert1(i1+i2,j1+j2,ick+kk2) + IF (ijkvert1 /= miss_v1) THEN + ! 3-dimensional distance + ! dist = SQRT(hdist**2 + zdist**2) + dists(ijk) = SQRT(ijdists(ij)**2 + (ijkvert1-vert2(k2,icj))**2) + ijkmask(ijk) = .TRUE. + END IF + ijk = ijk + 1 + END DO + ij = ij + 1 + END DO + END DO + ELSE + DO i2=-1,1 + DO j2=-1,1 + DO kk2=-1,1 + ijkvert1 = vert1(i1+i2,j1+j2,ick+kk2) + IF (ijkvert1 /= miss_v1) THEN + ! 3-dimensional distance + ! dist = SQRT(hdist**2 + zdist**2) + dists(ijk) = SQRT(ijdists(ij)**2 + (ijkvert1-vert2(k2,icj))**2) + ijkmask(ijk) = .TRUE. + END IF + ijk = ijk + 1 + END DO + ij = ij + 1 + END DO + END DO + END IF + + ! Filtering for dist == 0!! + wgs = 0. + IF (ANY(dists == 0.)) THEN + izero = MINLOC(dists) + wgs = 0. + wgs(izero(1)) = 1. + sumdists = 1. + ELSE + sumdists = SUM(1./dists, mask=ijkmask) + ! Applying mask in weights + DO k1=1, Na3 + IF (.NOT.ijkmask(k1)) THEN + wgs(k1) = 0. + ELSE + wgs(k1) = 1./dists(k1)/sumdists + END IF + END DO + totwght = SUM(wgs) + ! Fixing + IF (totwght < 1.) THEN + imaxwght = MAXLOC(wgs) + wgs(imaxwght(1)) = wgs(imaxwght(1)) + 1. - totwght + END IF + END IF + wghts(ick,ijk21Ncoords(ick),:) = wgs + + IF (printdbg) THEN + PRINT *, ' '//TRIM(fname) // ": i,j,k radar cell: ", & + ijk21coords(ick,ijk21Ncoords(ick),:), ' hgt: ', vert2(k2,icj), ' i,j,k reference: ', & + i1,j1,ick, ' dists: ', dists, ' mask: ', ijkmask, ' wghts: ', & + wghts(ick,ijk21Ncoords(ick),:) + END IF + END IF + END DO zloop + + END DO gridpts + + IF (totngridradar > ij21Ncoords*(ed2z-id2z+1) ) THEN + CALL ErrMsg('Too many radar cells were found! Theoretical maximum: ' // & + TRIM(ItoS(ij21Ncoords*(ed2z-id2z+1))) // ', but I found: ' // TRIM(ItoS(totngridradar)), & + fname, -1) + END IF + + + IF (dbg) THEN + PRINT *,' ' // TRIM(fname) // ": amount of vertical radar cells in WRF grid: ", i1, j1, & + ' = ', totngridradar + PRINT *, ' found along the reference column [amount (k_lev)]:', & + (TRIM(ItoS(ijk21Ncoords(k1))) // ' (' // TRIM(ItoS(k1)) // '), ', k1=id1z, ed1z) + END IF + + RETURN + + END SUBROUTINE weights_hcoordsR27around_polar_ellipsoide_1D + + SUBROUTINE compute_radar_lonlat_radial(Nangle, Nazimuth, Ndist, lon_pt, lat_pt, hgt_pt, ddist, & + max_dist, max_hgt, azimuths, a, f, elips, angle1d, dist1d, azimuth1d, lons, lats, hgts, lhgts, & + lon_bnds, lat_bnds, hgt_bnds, idbg) + ! Subroutine to copmute lon, lat, hgt coordinates in the radar matrix following polar coordinates + + IMPLICIT NONE + + INTEGER, INTENT(in) :: Nangle, Nazimuth, Ndist + REAL, INTENT(in) :: lon_pt, lat_pt, hgt_pt, ddist, max_dist, & + max_hgt + REAL, DIMENSION(Nazimuth), INTENT(in) :: azimuths + REAL(kind=rr8), INTENT(in) :: a, f + CHARACTER(len=S32), INTENT(in) :: elips + LOGICAL, OPTIONAL :: idbg + REAL, DIMENSION(Nangle), INTENT(out) :: angle1d + REAL, DIMENSION(Ndist), INTENT(out) :: dist1d + REAL, DIMENSION(Nazimuth), INTENT(out) :: azimuth1d + REAL, DIMENSION(Nangle,Ndist), INTENT(out) :: lons, lats + REAL, DIMENSION(Nazimuth,Ndist), INTENT(out) :: hgts + REAL, DIMENSION(Nazimuth,Ndist), INTENT(out) :: lhgts + REAL, DIMENSION(Nangle,Ndist,4),INTENT(out) :: lon_bnds, lat_bnds + REAL, DIMENSION(Nazimuth,Ndist,4),INTENT(out) :: hgt_bnds + + ! Local + INTEGER :: i, j, k, ijfound + REAL :: dangle, ddist2 + REAL :: angle, dist, azimuth, hgt, lhgt + REAL :: lonv, latv + REAL, DIMENSION(Nazimuth) :: radazimuths, dhgt, dhgt2 + LOGICAL :: printdbg + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!! Test inverse distance + REAL :: faz, s + + !!!!!!! Variables + ! Nangle, Nazimuth, Ndist: dimension of the radar space + ! lon_pt, lat_pt, hgt_pt: location of the radar + ! ddist: horizontal resolution [m] + ! max_dist: maximum distance of the radar [km] + ! max_hgt: maximum height of the radar [km] + ! azimuths: angles of the azimuths [degree] + ! a, f: axis / semi-axis of the ellipsoide + ! offset: offset to apply to the radar array + ! iis, jjs, kks: indices of the radar array + ! hdists: horizontal distances of the radar + ! zdists: vertical distances of the radar + ! lons, lats, hgts: coordinates of the radar array + ! lon_bnds, lat_bnds, hgt_bnds: vertices of the grids of the radar array + fname = 'compute_radar_lonlat_radial' + + dbg = set_dbg(idbg) + + ! Initializing + angle1d = fillValueR + dist1d = fillValueR + azimuth1d = azimuths + lhgts = fillValueR + lons = fillValueR + lats = fillValueR + hgts = fillValueR + lat_bnds = fillValueR + lon_bnds = fillValueR + hgt_bnds = fillValueR + + ijfound = 0 + ! Angle at the center of the ray + dangle = 360. / (Nangle) + ! Angle at the center of the ray + radazimuths = azimuths / DegRad + dhgt = 0. + dhgt(1:Nazimuth-1) = radazimuths(2:Nazimuth) - radazimuths(1:Nazimuth-1) + dhgt(Nazimuth) = dhgt(Nazimuth-1) + + ddist2 = ddist / 2. + dhgt2 = dhgt / 2. + + angl: DO i=1, Nangle + angle = dangle*(i-1) + angle1d(i) = angle + + DO j=1, Ndist + printdbg = .FALSE. + IF (dbg) THEN + IF (MOD(i+j,10) == 0) THEN + ijfound = ijfound + 1 + IF (ijfound < 10) printdbg = .TRUE. + END IF + END IF + + ! Centering at half distance + dist = ddist*(j-1) + ddist2 + IF (i == 1) dist1d(j) = dist + + ! Vertices + ! Cell boundaries + ! SW + dist = ddist*(j-1.) + angle = dangle*(i-1.-0.5) + CALL position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lonv, latv) + lon_bnds(i,j,1) = lonv + lat_bnds(i,j,1) = latv + + ! SE + dist = ddist*(j+1.) + angle = dangle*(i-1.-0.5) + CALL position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lonv, latv) + lon_bnds(i,j,2) = lonv + lat_bnds(i,j,2) = latv + + ! NE + dist = ddist*(j+1.) + angle = dangle*(i-1.+0.5) + CALL position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lonv, latv) + lon_bnds(i,j,3) = lonv + lat_bnds(i,j,3) = latv + + ! NW + dist = ddist*(j-1.) + angle = dangle*(i-1.+0.5) + CALL position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lonv, latv) + lon_bnds(i,j,4) = lonv + lat_bnds(i,j,4) = latv + + ! Getting the distance at the center of the radar cell for height computations + dist = ddist*(j-1) + ddist2 + CALL position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lonv, latv) + + lons(i,j) = lonv + lats(i,j) = latv + + IF (printdbg) PRINT *, ' ' // TRIM(fname), ' i,j ', i, j, ' <> distance ', dist, & + ' angle ', angle, ' lon: ', lons(i,j), ' lat: ', lats(i,j) + + DO k=1, Nazimuth + ! Centered at half angle + azimuth = radazimuths(k) + + CALL beam_height(hgt_pt, a, angle, azimuth, dist, hgt, lhgt) + IF (hgt > max_hgt) CYCLE + + hgts(k,j) = hgt + lhgts(k,j) = lhgt + + IF (printdbg) PRINT *, ' ', i, j, k, ' angle:', angle, ' dist:', dist, ' azimuth: ', & + azimuth, ' height: ', hgt + + ! Only to be computed at first angle + IF (i == 1 ) THEN + + ! Vertices + ! Cell boundaries + ! bottomW + CALL beam_height(hgt_pt, a, angle, azimuth-dhgt2(k), dist-ddist2, hgt, lhgt) + hgt_bnds(k,j,1) = hgt + + ! bottomE + CALL beam_height(hgt_pt, a, angle, azimuth-dhgt2(k), dist+ddist2, hgt, lhgt) + hgt_bnds(k,j,2) = hgt + + ! topE + CALL beam_height(hgt_pt, a, angle, azimuth+dhgt2(k), dist+ddist2, hgt, lhgt) + hgt_bnds(k,j,3) = hgt + + ! topW + CALL beam_height(hgt_pt, a, angle, azimuth+dhgt2(k), dist-ddist2, hgt, lhgt) + hgt_bnds(k,j,4) = hgt + + END IF + + END DO + END DO + + ! debugging at N, W, S, E axes + IF (i== 1 .OR. i == Nangle/4 .OR. i == Nangle/2 .OR. i == 3*Nangle/4 .OR. i == Nangle) THEN + PRINT *, ' ' // TRIM(fname) // ": values at major AXIS N;E;S;W angle = ", angle, ' _______' + PRINT *, ' lon first:', lons(i,1), ' last: ', lons(i,Ndist) + PRINT *, ' lat first:', lats(i,1), ' last: ', lats(i,Ndist) + PRINT *, ' hgt at first:', hgts(1,1), ' top: ', hgts(Nazimuth,1) + PRINT *, ' hgt at last:', hgts(1,Ndist), ' top: ', hgts(Nazimuth,Ndist) + PRINT *, ' values idist dist lon lat invangle invdist _______' + DO j=1, Ndist, Ndist/10 + CALL distance_between_points(lon_pt, lat_pt, lons(i,j), lats(i,j), a, f, elips, s, faz) + PRINT *,' ', j, ddist*(j-1) + ddist2, lons(i,j), lats(i,j), faz, s + + END DO + END IF + + END DO angl + + RETURN + + END SUBROUTINE compute_radar_lonlat_radial + + SUBROUTINE pertinence_polar(iwdx, ewdx, iwdy, ewdy, WRFlon, WRFlat, WRFlon_bnds, WRFlat_bnds, & + irdx, erdx, irdy, erdy, rlon, rlat, crx, cry, rmax, miss_v, maxNcoords, distradar, inradar, & + radarnum, ngridradar, Nc2, ijradar, ijdistradar, inside, idbg) + ! Subroutine to compute the 2-dimensional pertinence of a lonlat polar radial space (r) into + ! another regular lonlat (WRF) + ! NOTE: + ! - In radial space + ! dx: distance dy: angle + ! lon(dy,dx), lat(dy,dx) + + IMPLICIT NONE + + INTEGER, INTENT(in) :: iwdx, ewdx, iwdy, ewdy, maxNcoords, Nc2 + INTEGER, INTENT(in) :: irdx, erdx, irdy, erdy + REAL, INTENT(in) :: crx, cry, rmax, miss_v + REAL, DIMENSION(iwdx:ewdx,iwdy:ewdy), INTENT(in) :: WRFlon, WRFlat + REAL, DIMENSION(irdy:erdy,irdx:erdx), INTENT(in) :: rlon, rlat + REAL, DIMENSION(iwdx:ewdx,iwdy:ewdy,4), INTENT(in) :: WRFlon_bnds, WRFlat_bnds + LOGICAL, OPTIONAL :: idbg + INTEGER, DIMENSION(iwdx:ewdx,iwdy:ewdy), INTENT(out) :: radarnum + INTEGER, DIMENSION(iwdx:ewdx,iwdy:ewdy), INTENT(out) :: ngridradar + INTEGER, DIMENSION(iwdx:ewdx,iwdy:ewdy), & + INTENT(inout) :: inradar + REAL, DIMENSION(iwdx:ewdx,iwdy:ewdy), INTENT(out) :: distradar + INTEGER, DIMENSION(iwdx:ewdx,iwdy:ewdy,maxNcoords,Nc2), & + INTENT(out) :: ijradar + REAL, DIMENSION(iwdx:ewdx,iwdy:ewdy,maxNcoords), & + INTENT(out) :: ijdistradar + LOGICAL, INTENT(out) :: inside + + ! Local + INTEGER :: i, j, ir + INTEGER :: ic, jc, ijfound + INTEGER, DIMENSION(1) :: klonmin, klatmin, klonmax, klatmax, kmindist + INTEGER :: rdx2, rdy2, rdx + INTEGER, DIMENSION(2) :: ijtile_rad + REAL :: lonmin, lonmax, latmin, latmax, mindist + REAL :: lonrad, latrad, distrad + INTEGER :: totngridradar + REAL :: lonSW, lonNW, lonNE, lonSE + REAL :: latSW, latNW, latNE, latSE + REAL, DIMENSION(4) :: lonvertices, latvertices + REAL, DIMENSION(iwdx:ewdx,iwdy:ewdy) :: ijdist + LOGICAL :: centerin + LOGICAL :: dbg, printdbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! iwdx, ewdx, iwdy, ewdy: size of the WRF lonlathgt + ! rnum: number of radar from the list of radars to process + ! WRFlon, WRFlat: coordinates of the WRF lonlathgt [degrees_east, degrees_North] + ! WRFlon_bnds, WRFlat_bnds: coordinates of the vertices of the grids of the WRF lonlat + ! irdx, erdx, irdy, erdy: size of the polar r lonlathgt that we want to know how much fits inside + ! WRF lonlat + ! rlon, rlat: coordinates of the r lonlathgt [degrees_east, degrees_North] + ! crx, cry: coordinates of the center of the radial space + ! rmax: maximum radii from the center of the radial space + ! miss_v: missing value + ! rmax: maximum possible amount of r lonlathgt cells in a single WRFlonlat grid + ! + ! OUTPUT + ! distradar: distance in WRF lonlat grids to the center of r lonlat + ! inradar: whether a WRF lonlat grid holds inside any r lonlat cell (horizontal only) + ! radarnum: whether a WRF lonlat grid holds inside the r lonlat cell (horizontal only) + ! ngridradar: amount of r lonlat cells wihtin a WRF lonlat grid + ! ijkradar: indices of the polar r lonlat cells wihtin a WRF lonlat grid + ! ijdistradar: distance to the WRF lonlat grid of the r lonlat cells wihtin a WRF lonlat grid + fname = 'pertinence_polar' + + dbg = set_dbg(idbg) + + IF (dbg) THEN + PRINT *,' ' // TRIM(fname) // ": looking for radar at ", crx, cry, ' max dist: ', rmax + PRINT *,' size of WRF lonlat:', iwdx, ewdx, iwdy, ewdy + PRINT *,' size of r lonlat:', irdx, erdx, irdy, erdy + END IF + + ! Getting borders of the tile + lonSW = WRFlon(iwdx,iwdy) + lonNW = WRFlon(iwdx,ewdy) + lonNE = WRFlon(ewdx,ewdy) + lonSE = WRFlon(ewdx,iwdy) + latSW = WRFlat(iwdx,iwdy) + latNW = WRFlat(iwdx,ewdy) + latNE = WRFlat(ewdx,ewdy) + latSE = WRFlat(ewdx,iwdy) + + ! WR lon,lat matrices can be rotated, thus matrix vertices might not coincide... + lonvertices = (/ lonSW, lonNW, lonNE, lonSE /) + latvertices = (/ latSW, latNW, latNE, latSE /) + lonmin = MINVAL(lonvertices) + lonmax = MAXVAL(lonvertices) + latmin = MINVAL(latvertices) + latmax = MAXVAL(latvertices) + + klonmin = MINLOC(lonvertices) + klatmin = MINLOC(latvertices) + klonmax = MAXLOC(lonvertices) + klatmax = MAXLOC(latvertices) + + rdx2 = (erdx - irdx + 1) / 2 + rdy2 = (erdy - irdy + 1) / 2 + + ! Getting actual areas encompassed by a radar + distradar = fillValueR + inradar = fillValueI + radarnum = fillValueI + ngridradar = 0 + ijradar = fillValueI + ijdistradar = fillValueR + + !irng = 0 + centerin = .FALSE. + CALL radar_inside (iwdx, ewdx, iwdy, ewdy, WRFlon, WRFlat, WRFlon_bnds, WRFlat_bnds, crx, cry, & + rmax, miss_v, distradar, inradar, radarnum, inside, dbg) + + ! Getting lon,lat for the radar matrix + ! Can we know from the radar perspective which of its grid points are inside this tile? + ! Even if I could know that, radar matrix is not (by now) parallelized ! + ! Therefore, I need to compute the whole radar matrix at each tile + + totngridradar = 0 + ijfound = 1 + !IF (rlon > lonvertices(klonmin(1)) .AND. rlon < lonvertices(klonmax(1)) .AND. & + ! rlat > latvertices(klonmin(1)) .AND. rlat < latvertices(klatmax(1))) THEN + insideradar: IF (inside) THEN + + centerin = .TRUE. + + ! Looping all over the radar matrix + lonlat_radar: DO i=irdy, erdy + DO j=irdx, erdx +! lonlat_radar: DO i=1,erdy +! DO j=irdx, 10 + printdbg = .FALSE. + IF (dbg) THEN + IF (MOD(i+j,10) == 0 .AND. ijfound < 10) THEN + printdbg = .TRUE. + ijfound = ijfound + 1 + END IF + END IF + + ! Avoid values outside radar range + IF (rlon(i,j) >= fillValueR*0.95 .OR. rlat(i,j) >= fillValueR*0.95) CYCLE + + ins: IF (rlon(i,j) > lonvertices(klonmin(1)) .AND. rlon(i,j) < lonvertices(klonmax(1)) .AND. & + rlat(i,j) > latvertices(klonmin(1)) .AND. rlat(i,j) < latvertices(klatmax(1))) THEN + + totngridradar = totngridradar + 1 + ! Closest WRF grid point to the radar cell + ijdist=SQRT((WRFlon-rlon(i,j))**2+(WRFlat-rlat(i,j))**2) + ijtile_rad = MINLOC(ijdist) + ic = ijtile_rad(1) + jc = ijtile_rad(2) + + inradar(ic,jc) = 1 + radarnum(ic,jc) = 1 + + ! Amount of r cells in this WRF grid point + ngridradar(ic,jc) = ngridradar(ic,jc) + 1 + IF (ngridradar(ic,jc) > maxNcoords) THEN + CALL ErrMsg("For WRF's i,j: " // TRIM(ItoS(ic)) //', ' // TRIM(ItoS(jc)) // & + " too much radar cells: " // TRIM(ItoS(ngridradar(ic,jc))) // " in a WRF grid. " // & + " distance to the radar: " // TRIM(RtoS(MINVAL(ijdist))) // & + "'ngrid_radar' in namelist current value: " // TRIM(ItoS(maxNcoords)), fname, -1) + END IF + + ! Coordinates of the radar cells within this WRF grid + ijradar(ic,jc,ngridradar(ic,jc), 1) = i + ijradar(ic,jc,ngridradar(ic,jc), 2) = j + ijdistradar(ic,jc,ngridradar(ic,jc)) = MINVAL(ijdist) + + IF (printdbg) THEN + PRINT *,' ' // TRIM(fname) // ": check total amount of radar cells wihtin " // & + "WRF grid: ", totngridradar + PRINT *, ' > WRF grid:', ijtile_rad + PRINT *, ' > amount radar cells in WRF grid:', ngridradar(ic,jc) + PRINT *, ' > radar cell coordinates of last found:', i,j + END IF + + END IF ins + END DO + END DO lonlat_radar + + rdx = erdx-irdx+1 + IF (dbg) THEN + PRINT *,' WRF lonlat box ...' + PRINT *,' SW:', lonvertices(klonmin(1)), latvertices(klatmin(1)) + PRINT *,' NE:', lonvertices(klonmax(1)), latvertices(klatmax(1)) + PRINT *,' ' // TRIM(fname) // ": with radar at ", crx,', ',cry, "; amount of " // & + "radar cells in WRF grid: ", totngridradar, ' theoretical value:', & + INT(pict*((rdx)/2)**2), INT(SUM(ngridradar)) + ijtile_rad = MAXLOC(ngridradar) + ic = ijtile_rad(1) + jc = ijtile_rad(2) + ijdist=SQRT((WRFlon(ic,jc)-crx)**2+(WRFlat(ic,jc)-cry)**2) + PRINT *,' largest amount radar cells in a WRF grid:', MAXVAL(ngridradar), ' found at: ', & + ijtile_rad, ' distance to the radar:', ijdist(1,1), ' degree' + END IF + + IF (totngridradar < INT((pict*(rdx/2)**2)/10)) THEN + CALL ErrMSg(' found radar cells in WRF grid: ' // TRIM(ItoS(totngridradar)) // & + ' theoretical: ' // TRIM(ItoS(INT(pict*(rdx/2)**2))) // ' Not enough ' // & + 'radar cells were found',fname,-1) + END IF + + ELSE + IF (dbg) THEN + PRINT *,' ' // TRIM(fname) // ": radar at ", crx,', ', cry, "; outside WRF grid !!" + PRINT *,' WRF lonlat box ...' + PRINT *,' SW:', lonvertices(klonmin(1)), latvertices(klatmin(1)) + PRINT *,' NE:', lonvertices(klonmax(1)), latvertices(klatmax(1)) + END IF + + END IF insideradar + + RETURN + + END SUBROUTINE pertinence_polar + + SUBROUTINE read_lut(lutfilen, lutgrid, idbg) + ! Subroutine to read the T-MATRIX look-up-table (LUT) from an specific microphysics from a netCDF + ! file + + IMPLICIT NONE + + CHARACTER(len=S256), INTENT(in) :: lutfilen + LOGICAL, OPTIONAL :: idbg + TYPE(lut_vars), INTENT(out) :: lutgrid + + ! Local + INTEGER :: lutid, ierr + INTEGER :: Nlutdims, Nlutvars + INTEGER :: Nbnds, Nlevs, Nqbins + CHARACTER(len=S16) :: S16a + CHARACTER(len=S32) :: S32a + CHARACTER(len=S256) :: S256a + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lutfilen: name of the netCDF file with the values + ! lutgrid: derived type where to kept the variables + + fname = 'read_lut' + + dbg = set_dbg(idbg) + + CALL open_NCfile(lutfilen, lutid, dbg) + + CALL get_Ndims_ncid(lutid, Nlutdims) + CALL get_Nvars_ncid(lutid, Nlutvars) + IF (dbg) PRINT *,' ' // TRIM(fname) // ": getting from file '" // TRIM(lutfilen) // & + "' Ndims: ", Nlutdims, ' Nvars: ', Nlutvars + + S32a = 'WRF_mp_physics' + CALL get_varS0D_ncid(lutid, 256, S32a, S256a, dbg) + S16a = TRIM(S256a) + lutgrid%mp_scheme = S16a + + ! LUT space + S32a = 'elevation' + CALL get_dimlength_ncid(lutid, S32a, Nlevs) + S32a = 'radar_band' + CALL get_dimlength_ncid(lutid, S32a, Nbnds) + S32a = 'qxbin' + CALL get_dimlength_ncid(lutid, S32a, Nqbins) + + ! I do think it makes the code more readable + lutgrid%Nlutelevs = Nlevs + lutgrid%Nbands = Nbnds + lutgrid%Nlutqxbins = Nqbins + + !! Variables + !!! RAIN + IF (ALLOCATED(lutgrid%Zh_RAIN)) DEALLOCATE(lutgrid%Zh_RAIN) + ALLOCATE(lutgrid%Zh_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zh_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zh_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zh_RAIN, dbg) + IF (ALLOCATED(lutgrid%Zdr_RAIN)) DEALLOCATE(lutgrid%Zdr_RAIN) + ALLOCATE(lutgrid%Zdr_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zdr_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zdr_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zdr_RAIN, dbg) + IF (ALLOCATED(lutgrid%LDR_RAIN)) DEALLOCATE(lutgrid%LDR_RAIN) + ALLOCATE(lutgrid%LDR_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'LDR_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'LDR_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%LDR_RAIN, dbg) + IF (ALLOCATED(lutgrid%Aih_RAIN)) DEALLOCATE(lutgrid%Aih_RAIN) + ALLOCATE(lutgrid%Aih_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aih_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aih_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aih_RAIN, dbg) + IF (ALLOCATED(lutgrid%Aiv_RAIN)) DEALLOCATE(lutgrid%Aiv_RAIN) + ALLOCATE(lutgrid%Aiv_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aiv_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aiv_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aiv_RAIN, dbg) + IF (ALLOCATED(lutgrid%KDP_RAIN)) DEALLOCATE(lutgrid%KDP_RAIN) + ALLOCATE(lutgrid%KDP_RAIN(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'KDP_RAIN', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'KDP_RAIN' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%KDP_RAIN, dbg) + + !!! SNOW + IF (ALLOCATED(lutgrid%Zh_SNOW)) DEALLOCATE(lutgrid%Zh_SNOW) + ALLOCATE(lutgrid%Zh_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zh_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zh_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zh_SNOW, dbg) + IF (ALLOCATED(lutgrid%Zdr_SNOW)) DEALLOCATE(lutgrid%Zdr_SNOW) + ALLOCATE(lutgrid%Zdr_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zdr_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zdr_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zdr_SNOW, dbg) + IF (ALLOCATED(lutgrid%LDR_SNOW)) DEALLOCATE(lutgrid%LDR_SNOW) + ALLOCATE(lutgrid%LDR_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'LDR_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'LDR_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%LDR_SNOW, dbg) + IF (ALLOCATED(lutgrid%Aih_SNOW)) DEALLOCATE(lutgrid%Aih_SNOW) + ALLOCATE(lutgrid%Aih_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aih_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aih_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aih_SNOW, dbg) + IF (ALLOCATED(lutgrid%Aiv_SNOW)) DEALLOCATE(lutgrid%Aiv_SNOW) + ALLOCATE(lutgrid%Aiv_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aiv_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aiv_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aiv_SNOW, dbg) + IF (ALLOCATED(lutgrid%KDP_SNOW)) DEALLOCATE(lutgrid%KDP_SNOW) + ALLOCATE(lutgrid%KDP_SNOW(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'KDP_SNOW', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'KDP_SNOW' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%KDP_SNOW, dbg) + + !!! GRAU + IF (ALLOCATED(lutgrid%Zh_GRAU)) DEALLOCATE(lutgrid%Zh_GRAU) + ALLOCATE(lutgrid%Zh_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zh_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zh_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zh_GRAU, dbg) + IF (ALLOCATED(lutgrid%Zdr_GRAU)) DEALLOCATE(lutgrid%Zdr_GRAU) + ALLOCATE(lutgrid%Zdr_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zdr_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zdr_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zdr_GRAU, dbg) + IF (ALLOCATED(lutgrid%LDR_GRAU)) DEALLOCATE(lutgrid%LDR_GRAU) + ALLOCATE(lutgrid%LDR_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'LDR_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'LDR_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%LDR_GRAU, dbg) + IF (ALLOCATED(lutgrid%Aih_GRAU)) DEALLOCATE(lutgrid%Aih_GRAU) + ALLOCATE(lutgrid%Aih_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aih_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aih_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aih_GRAU, dbg) + IF (ALLOCATED(lutgrid%Aiv_GRAU)) DEALLOCATE(lutgrid%Aiv_GRAU) + ALLOCATE(lutgrid%Aiv_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aiv_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aiv_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aiv_GRAU, dbg) + IF (ALLOCATED(lutgrid%KDP_GRAU)) DEALLOCATE(lutgrid%KDP_GRAU) + ALLOCATE(lutgrid%KDP_GRAU(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'KDP_GRAU', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'KDP_GRAU' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%KDP_GRAU, dbg) + + !!! ICE + IF (ALLOCATED(lutgrid%Zh_ICE)) DEALLOCATE(lutgrid%Zh_ICE) + ALLOCATE(lutgrid%Zh_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zh_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zh_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zh_ICE, dbg) + IF (ALLOCATED(lutgrid%Zdr_ICE)) DEALLOCATE(lutgrid%Zdr_ICE) + ALLOCATE(lutgrid%Zdr_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Zdr_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Zdr_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Zdr_ICE, dbg) + IF (ALLOCATED(lutgrid%LDR_ICE)) DEALLOCATE(lutgrid%LDR_ICE) + ALLOCATE(lutgrid%LDR_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'LDR_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'LDR_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%LDR_ICE, dbg) + IF (ALLOCATED(lutgrid%Aih_ICE)) DEALLOCATE(lutgrid%Aih_ICE) + ALLOCATE(lutgrid%Aih_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aih_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aih_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aih_ICE, dbg) + IF (ALLOCATED(lutgrid%Aiv_ICE)) DEALLOCATE(lutgrid%Aiv_ICE) + ALLOCATE(lutgrid%Aiv_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'Aiv_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'Aiv_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%Aiv_ICE, dbg) + IF (ALLOCATED(lutgrid%KDP_ICE)) DEALLOCATE(lutgrid%KDP_ICE) + ALLOCATE(lutgrid%KDP_ICE(Nbnds,Nlevs,Nqbins), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'KDP_ICE', 3, (/Nbnds,Nlevs,Nqbins/), fname, ierr) + S32a = 'KDP_ICE' + CALL get_varD3D_ncid(lutid, 1, 1, 1, Nbnds, Nlevs, Nqbins, S32a, lutgrid%KDP_ICE, dbg) + +! From namelist or from file? +! S32a = 'q_rain_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_rain_min, dbg) +! S32a = 'q_rain_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_rain_max, dbg) +! S32a = 'qn_rain_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_rain_min, dbg) +! S32a = 'qn_rain_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_rain_max, dbg) +! S32a = 'q_snow_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_snow_min, dbg) +! S32a = 'q_snow_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_snow_max, dbg) +! S32a = 'qn_snow_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_snow_min, dbg) +! S32a = 'qn_snow_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_snow_max, dbg) +! S32a = 'q_grau_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_grau_min, dbg) +! S32a = 'q_grau_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_grau_max, dbg) +! S32a = 'qn_grau_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_grau_min, dbg) +! S32a = 'qn_grau_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_grau_max, dbg) +! S32a = 'q_ice_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_ice_min, dbg) +! S32a = 'q_ice_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_ice_max, dbg) +! S32a = 'qn_ice_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_ice_min, dbg) +! S32a = 'qn_ice_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%qn_ice_max, dbg) +! S32a = 'q_clou_min' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_clou_min, dbg) +! S32a = 'q_clou_max' +! CALL get_varR0D_ncid(lutid, S32a, lutgrid%q_clou_max, dbg) +! S32a = 'dim_ice' +! CALL get_varN0D_ncid(lutid, S32a, lutgrid%dim_ice, dbg) +! S32a = 'dimv' +! CALL get_varN0D_ncid(lutid, S32a, lutgrid%dimv, dbg) + + !! Getting densities + S32a = 'rhoo' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%rhoo, dbg) + S32a = 'rho_r' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%rho_r, dbg) + S32a = 'n0r' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%n0r, dbg) + S32a = 'rho_s' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%rho_s, dbg) + S32a = 'n0s' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%n0s, dbg) + S32a = 'rho_g' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%rho_g, dbg) + S32a = 'n0g' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%n0g, dbg) + S32a = 'rho_h' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%rho_h, dbg) + S32a = 'n0h' + CALL get_varR0D_ncid(lutid, S32a, lutgrid%n0h, dbg) + + IF (ALLOCATED(lutgrid%elevations)) DEALLOCATE(lutgrid%elevations) + ALLOCATE(lutgrid%elevations(Nlevs), STAT=ierr) + IF (ierr /= 0) CALL ErrMsgAlloc('', 'elevations', 1, (/Nlevs/), fname, ierr) + + S32a = 'elevation' + CALL get_varR1D_ncid(lutid, 1, Nlevs, S32a, lutgrid%elevations, dbg) + + CALL close_NCfile(lutfilen, lutid, dbg) + + RETURN + + END SUBROUTINE read_lut + + SUBROUTINE distance_between_points(lon1_pt, lat1_pt, lon2_pt, lat2_pt, a, f, elips, dist, angle) + ! Subroutine tro compute the coordinates at a given distance and angle from a given point + + IMPLICIT NONE + + REAL, INTENT(in) :: lon1_pt, lat1_pt, lon2_pt, lat2_pt + REAL(kind=rr8), INTENT(in) :: a, f + CHARACTER(len=30), INTENT(in) :: elips + REAL, INTENT(out) :: dist, angle + + ! Local + REAL(kind=rr8) :: ddist + INTEGER(kind=ii4) :: lon_sign, lat_sign, angle_sign + REAL(kind=rr8) :: lon_deg, lat_deg, angle_deg, lon_min, & + lat_min, angle_min, lon_sec, lat_sec, angle_sec + REAL(kind=rr8) :: glat1, glon1, faz, distd + REAL(kind=rr8) :: glat2, glon2, baz + INTEGER :: Niter, kindv + REAL(kind=rr8) :: signv, lamv + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon1_pt, lat1_pt: coordinates of the first point [deg] + ! lon2_pt, lat2_pt: coordinates of the second point [deg] + ! a: semimajor axis equatorial (in meters) + ! f: flattening + ! elips: name of the ellipsis being used + ! dist: distance [m] + ! angle: angle [deg] + + fname = 'distance_between_points' + + glat1 = lat1_pt / drad + glon1 = lon1_pt / drad + glat2 = lat2_pt / drad + glon2 = lon2_pt / drad + + CALL inverse(a, 1.d0/f, glat1, glon1, glat2, glon2, faz, baz, distd, Niter, signv, lamv, kindv) + + angle = faz*DegRad + dist = distd*1. + + RETURN + + END SUBROUTINE distance_between_points + + SUBROUTINE inverse(a, rf, b1, l1, b2, l2, faz, baz, s, it, sig, lam, kindv) + ! Subroutine to compute distance between 2 points + ! + ! Following software by NOAA + ! https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html + ! + ! Bessel, F. W., (1826) On the computation of geographical longitude and latitude from geodetic + ! measurements, Astronomische Nachrichten (Astronomical Notes), Band 4, No. 86, Spalten 241-254, + ! Altona 1826 + ! Vincenty, T., (1975) Direct and inverse solutions on the ellipsoid with application of nested + ! equations, Survey Review, Vol. 23, No. 176, pp. 88-93. + ! Vincenty, T., (1976) Correspondence: solutions of geodesics, survey Review, Vol. 23, No. 180, + ! p. 294. + + ! + !*** inverse for long-line and antipodal cases. + !*** latitudes may be 90 degrees exactly. + !*** latitude positive north, longitude positive east, radians. + !*** azimuth clockwise from north, radians. + !*** original programmed by thaddeus vincenty, 1975, 1976 + !*** removed back side solution option, debugged, revised -- 2011may01 -- dgm + !*** this version of code is interim -- antipodal boundary needs work + + !*** output (besides faz,baz, and s): + !*** it, iteration count + !*** sig, spherical distance on auxiliary sphere + !*** lam, longitude dIFference on auxiliary sphere + !*** kind, solution flag: kind=1, long-line; kind=2, antipodal + + IMPLICIT NONE + + REAL(kind=rr8), INTENT(in) :: a, rf + REAL(kind=rr8), INTENT(in) :: b1, b2, l1, l2 + REAL(kind=rr8), INTENT(out) :: faz, baz, s + INTEGER, INTENT(out) :: it, kindv + REAL(kind=rr8), INTENT(out) :: sig, lam + + ! Local + INTEGER :: i, iter, Niter + REAL(kind=rr8) :: boa, l + REAL(kind=rr8) :: biga, bigb, bige, bigf, c, d, z + REAL(kind=rr8) :: cosal2, coslam, cossig, costm, costm2 + REAL(kind=rr8) :: dsig, ep2, eps + REAL(kind=rr8) :: prev, tem1, tem2, temp, test + REAL(kind=rr8) :: sinal, sinlam, sinsig + REAL(kind=rr8) :: beta1, sinu1, cosu1 + REAL(kind=rr8) :: beta2, sinu2, cosu2 + LOGICAL :: converged, keeplooping + CHARACTER(len=S32) :: fname + + fname = 'inverse' + + boa = 1.d0-1.d0/rf + !***** sinu1 = boa*dsin(b1)/dsqrt((boa*dsin(b1))**2+dcos(b1)**2) + !***** cosu1 = dsqrt(-sinu1**2+1.d0) !*** rounDOff + !***** sinu2 = boa*dsin(b2)/dsqrt((boa*dsin(b2))**2+dcos(b2)**2) + !***** cosu2 = dsqrt(-sinu2**2+1.d0) !*** rounDOff + + !*** better reduced latitude + beta1 = datan(boa*dtan(b1)) + sinu1 = dsin(beta1) + cosu1 = dcos(beta1) + + !*** better reduced latitude + beta2 = datan(boa*dtan(b2)) + sinu2 = dsin(beta2) + cosu2 = dcos(beta2) + + !*** longitude dIFference [-pi,pi] + l = l2-l1 + IF (l .gt. dpict) l = l - dpict - dpict + IF (l .lt. -dpict) l = l + dpict + dpict + prev = l + test = l + kindv = 1 + !*** v13 (rapp part II) + lam = l + + !*** top of the long-line loop (kind=1) + it = 0 + Niter = 100 + + DO i=1, Niter + + sinlam = dsin(lam) + + !*** no--troublesome + !***** IF(dabs(pi-dabs(l)).lt.0.2d-10) sinlam=0.d0 + coslam = dcos(lam) + temp = cosu1*sinu2-sinu1*cosu2*coslam + !*** v14 (rapp part II) + sinsig = dsqrt((cosu2*sinlam)**2+temp**2) + !*** v15 (rapp part II) + cossig = sinu1*sinu2 + cosu1*cosu2*coslam + sig = datan2(sinsig,cossig) + + !*** v17 (rapp part II) + !***** sinal = cosu1*cosu2*sinlam/sinsig + IF(dabs(sinsig).lt.eps) THEN + !*** avoid div 0 + sinal = cosu1*cosu2*sinlam/dsign(eps,sinsig) + ELSE + sinal = cosu1*cosu2*sinlam/sinsig + ENDIF + + cosal2 = -sinal**2+1.d0 + !*** v18 (rapp part II) + !***** costm = -2.d0*sinu1*sinu2/(cosal2+tol)+cossig + IF(dabs(cosal2).lt.eps) THEN + !*** avoid div 0 + costm = -2.d0*(sinu1*sinu2/dsign(eps,cosal2))+cossig + ELSE + costm = -2.d0*(sinu1*sinu2/cosal2)+cossig + ENDIF + costm2=costm*costm + !*** v10 (rapp part II) + c=((-3.d0*cosal2+4.d0)/rf+4.d0)*cosal2/rf/16.d0 + + converged = .FALSE. + DO it=1, Niter + + !*** v11 + d=(((2.d0*costm2-1.d0)*cossig*c+costm)*sinsig*c+sig)*(1.d0-c)/rf + + IF (kindv .eq. 1) THEN + lam = l+d*sinal + IF (DABS(lam-test) .lt. tol) THEN + converged = .TRUE. + EXIT + + END IF + + IF (DABS(lam) .gt. dpict) THEN + + keeplooping = .FALSE. + + kindv = 2 + lam = dpict + IF (l .lt. 0.d0) lam = -lam + sinal = 0.d0 + cosal2 = 1.d0 + test = 2.d0 + prev = test + sig = dpict - DABS(DATAN(sinu1/cosu1) + DATAN(sinu2/cosu2)) + sinsig = DSIN(sig) + cossig = DCOS(sig) + c = ((-3.d0*cosal2+4.d0)/rf+4.d0)*cosal2/rf/16.d0 + IF (DABS(sinal-prev) .lt. tol) THEN + converged = .TRUE. + EXIT + END IF + !***** costm = -2.d0*sinu1*sinu2/(cosal2+tol)+cossig + IF (DABS(cosal2) .lt. eps) THEN + !*** avoid div 0 + costm= -2.d0*(sinu1*sinu2/dsign(eps,cosal2))+cossig + ELSE + costm= -2.d0*(sinu1*sinu2/cosal2)+cossig + ENDIF + costm2=costm*costm + keeplooping = .TRUE. + ENDIF + + IF ( .NOT. keeplooping) THEN + + IF ( ((lam-test)*(test-prev)) .lt. 0.d0 .AND. it .gt. 5) lam = (2.d0*lam+3.d0*test + & + prev) / 6.d0 + + prev = test + test = lam + + EXIT + + END IF + + ELSE + + sinal= (lam-l)/d + !*** refined converge. + IF ( ((sinal-test)*(test-prev)) .lt. 0.d0 .AND. it.gt.5) sinal = (2.d0*sinal+3.d0*test + & + prev)/6.d0 + + prev = test + test = sinal + + cosal2 = -sinal**2+1.d0 + sinlam = sinal*sinsig/(cosu1*cosu2) + coslam = -DSQRT(DABS(-sinlam**2 +1.d0)) + lam = DATAN2(sinlam,coslam) + temp = cosu1*sinu2 - sinu1*cosu2*coslam + sinsig = DSQRT((cosu2*sinlam)**2+temp**2) + cossig = sinu1*sinu2 + cosu1*cosu2*coslam + sig = DATAN2(sinsig,cossig) + c = ((-3.d0*cosal2+4.d0)/rf+4.d0)*cosal2/rf/16.d0 + IF (DABS(sinal-prev) .lt. tol) THEN + converged = .TRUE. + EXIT + END IF + + !***** costm = -2.d0*sinu1*sinu2/(cosal2+tol)+cossig + IF (DABS(cosal2) .lt. eps) THEN + !*** avoid div 0 + costm = -2.d0*(sinu1*sinu2/dsign(eps,cosal2)) + cossig + ELSE + costm = -2.d0*(sinu1*sinu2/cosal2)+cossig + ENDIF + + costm2=costm*costm + + ENDIF + + END DO + + IF (converged) EXIT + + END DO + + !*** convergence + IF (kindv .eq. 2) THEN + !*** antipodal + faz= sinal/cosu1 + baz= DSQRT(-faz**2+1.d0) + IF (temp .lt. 0.d0) baz = -baz + faz = DATAN2(faz,baz) + tem1 = -sinal + tem2 = sinu1*sinsig - cosu1*cossig*baz + baz = DATAN2(tem1,tem2) + ELSE + !*** long-line + tem1 = cosu2*sinlam + tem2 = cosu1*sinu2 - sinu1*cosu2*coslam + faz = DATAN2(tem1,tem2) + tem1 = -cosu1*sinlam + tem2 = sinu1*cosu2 - cosu1*sinu2*coslam + baz = DATAN2(tem1,tem2) + ENDIF + + IF (faz .lt. 0.d0) faz = faz + dpict + dpict + IF (baz .lt. 0.d0) baz = baz + dpict + dpict + + !*** Helmert 1880 from Vincenty "Geodetic inverse solution between antipodal points" + ep2 = 1.d0/(boa*boa)-1.d0 + !*** 15 + bige = DSQRT(1.d0+ep2*cosal2) + !*** 16 + bigf = (bige-1.d0)/(bige+1.d0) + !*** 17 + biga = (1.d0+bigf*bigf/4.d0)/(1.d0-bigf) + !*** 18 + bigb = bigf*(1.d0-0.375d0*bigf*bigf) + + z =bigb/6.d0*costm*(-3.d0+4.d0*sinsig**2)*(-3.d0+4.d0*costm2) + + !*** 19 + dsig = bigb*sinsig*(costm+bigb/4.d0*(cossig*(-1.d0+2.d0*costm2)-z)) + !*** 20 + s = (boa*a)*biga*(sig-dsig) + + RETURN + + END SUBROUTINE inverse + + SUBROUTINE position_from_point(lon_pt, lat_pt, a, f, elips, dist, angle, lon_dist, lat_dist) + ! Subroutine to compute the coordinates at a given distance and angle from a given point + ! + ! Following software by NOAA + ! https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html + ! + ! Bessel, F. W., (1826) On the computation of geographical longitude and latitude from geodetic + ! measurements, Astronomische Nachrichten (Astronomical Notes), Band 4, No. 86, Spalten 241-254, + ! Altona 1826 + ! Vincenty, T., (1975) Direct and inverse solutions on the ellipsoid with application of nested + ! equations, Survey Review, Vol. 23, No. 176, pp. 88-93. + ! Vincenty, T., (1976) Correspondence: solutions of geodesics, survey Review, Vol. 23, No. 180, + ! p. 294. + + IMPLICIT NONE + + REAL, INTENT(in) :: lon_pt, lat_pt, dist, angle + REAL(kind=rr8), INTENT(in) :: a, f + CHARACTER(len=30), INTENT(in) :: elips + REAL, INTENT(out) :: lat_dist, lon_dist + + ! Local + REAL(kind=rr8) :: ddist + INTEGER(kind=ii4) :: lon_sign, lat_sign, angle_sign + REAL(kind=rr8) :: lon_deg, lat_deg, angle_deg, lon_min, & + lat_min, angle_min, lon_sec, lat_sec, angle_sec + REAL(kind=rr8) :: glat1, glon1, faz + REAL(kind=rr8) :: glat2, glon2, baz + INTEGER(kind=ii4) :: glon2_sign, glat2_sign + REAL(kind=rr8) :: dms_angle + INTEGER(kind=ii4) :: glon2_deg, glat2_deg, glon2_min, & + glat2_min + REAL(kind=rr8) :: glon2_sec, glat2_sec + REAL(kind=rr8) :: lat2, lon2 + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon_pt, lat_pt: coordinates of the point + ! a: semimajor axis equatorial (in meters) + ! f: flattening + ! elips: name of the ellipsis being used + ! dist: distance [m] + ! angle: angle [deg] + fname = 'position_from_point' + + glat1 = lat_pt / drad + glon1 = lon_pt / drad + faz = angle / drad + + ddist = dist*1.d0 + + ! Assuming we are well behaved... + ! ! Checking distance + ! IF ( ddist > dd_max ) THEN + ! PRINT (errormsg) + ! PRINT *, ' Invalid Distance: ', ddist, ' ... Please re-enter it ' + ! PRINT *, " maximum distance for projection '" // TRIM(elips) // "' is ", dd_max, ' m' + ! STOP -1 + ! END IF + + !PRINT *,' ' // TRIM(fname) // ": ", a, f, glat1, glon1, ddist + CALL dirct1(a, f, glat1, glon1, glat2, glon2, faz, baz, ddist) + + lat_dist = glat2*drad + lon_dist = glon2*drad + + RETURN + + END SUBROUTINE position_from_point + + SUBROUTINE DIRCT1(A,F,GLAT1,GLON1,GLAT2,GLON2,FAZ,BAZ,S) +! +! *** SOLUTION OF THE GEODETIC DIRECT PROBLEM AFTER T.VINCENTY +! *** MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS +! *** EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL +! +! *** A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID +! *** F IS THE FLATTENING OF THE REFERENCE ELLIPSOID +! *** LATITUDES AND LONGITUDES IN RADIANS POSITIVE NORTH AND EAST +! *** AZIMUTHS IN RADIANS CLOCKWISE FROM NORTH +! *** GEODESIC DISTANCE S ASSUMED IN UNITS OF SEMI-MAJOR AXIS A +! +! *** PROGRAMMED FOR CDC-6600 BY LCDR L.PFEIFER NGS ROCKVILLE MD 20FEB75 +! *** MODIFIED FOR SYSTEM 360 BY JOHN G GERGEN NGS ROCKVILLE MD 750608 +! +! Following software by NOAA +! https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html +! +! Bessel, F. W., (1826) On the computation of geographical longitude and latitude from geodetic +! measurements, Astronomische Nachrichten (Astronomical Notes), Band 4, No. 86, Spalten 241-254, +! Altona 1826 +! Vincenty, T., (1975) Direct and inverse solutions on the ellipsoid with application of nested +! equations, Survey Review, Vol. 23, No. 176, pp. 88-93. +! Vincenty, T., (1976) Correspondence: solutions of geodesics, survey Review, Vol. 23, No. 180, +! p. 294. + + IMPLICIT NONE + + REAL(kind=rr8), INTENT(in) :: GLAT1, GLON1, FAZ, A, F, S + REAL(kind=rr8), INTENT(out) :: GLAT2, GLON2, BAZ + + ! Local + INTEGER :: i + REAL(kind=rr8) :: C, D, EPS + REAL(kind=rr8) :: R, TU, SF, CF, CU, SU, SA, C2A, X, E + REAL(kind=rr8) :: Y, SY, CY, CZ + LOGICAL :: converged + + CHARACTER(len=S32) :: fname + + DATA EPS /0.5D-13/ + + fname = 'DIRCT1' + + R=1.-F + TU=R*DSIN(GLAT1)/DCOS(GLAT1) + SF=DSIN(FAZ) + CF=DCOS(FAZ) + BAZ=0. + IF (CF.NE.0.) BAZ=DATAN2(TU,CF)*2. + + CU=1./DSQRT(TU*TU+1.) + SU=TU*CU + SA=CU*SF + C2A=-SA*SA+1. + X=DSQRT((1./R/R-1.)*C2A+1.)+1. + X=(X-2.)/X + C=1.-X + C=(X*X/4.+1)/C + D=(0.375D0*X*X-1.)*X + TU=S/R/A/C + Y=TU + + converged = .FALSE. + DO i=1, 100 + SY=DSIN(Y) + CY=DCOS(Y) + CZ=DCOS(BAZ+Y) + E=CZ*CZ*2.-1. + C=Y + X=E*CY + Y=E+E-1. + Y=(((SY*SY*4.-3.)*Y*CZ*D/6.+X)*D/4.-CZ)*SY*D+TU + IF (DABS(Y-C).LE.EPS) THEN + converged = .TRUE. + EXIT + END IF + END DO + IF (.NOT.converged) THEN + PRINT *, TRIM(errormsg) + PRINT *, TRIM(fname) // ': computation of Y,C did not converged !!' + PRINT *, ' DABS(Y-C)> EPS:', DABS(Y-C), ' > ', EPS + STOP -1 + END IF + + BAZ=CU*CY*CF-SU*SY + C=R*DSQRT(SA*SA+BAZ*BAZ) + D=SU*CY+CU*SY*CF + GLAT2=DATAN2(D,C) + C=CU*CY-SU*SY*CF + X=DATAN2(SY*SF,C) + C=((-3.*C2A+4.)*F+4.)*C2A*F/16. + D=((E*CY*C+CZ)*SY*C+Y)*SA + GLON2=GLON1+X-(1.-C)*D*F + BAZ=DATAN2(SA,BAZ)+dpict + + RETURN + + END SUBROUTINE DIRCT1 + + SUBROUTINE radar_inside (ims, ime, jms, jme, lon, lat, bnds_lon, bnds_lat, rlon, rlat, rdist, & + missing_value, gatehdist, inradar, radnum, inside, idbg) + ! Subroutine to determine if a radar fells inside the tile + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ims, ime, jms, jme + REAL, DIMENSION(ims:ime,jms:jme), INTENT(in) :: lon, lat + REAL, DIMENSION(ims:ime,jms:jme,4), INTENT(in) :: bnds_lon, bnds_lat + REAL, INTENT(in) :: rlon, rlat, rdist, missing_value + LOGICAL, OPTIONAL :: idbg + REAL, DIMENSION(ims:ime,jms:jme), INTENT(out) :: gatehdist + INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(out) :: inradar, radnum + LOGICAL, INTENT(out) :: inside + + ! Local + INTEGER :: i, j + INTEGER :: ii, ie, ji, je + INTEGER :: ibnd + LOGICAL, DIMENSION(ims:ime,jms:jme) :: mask + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon, lat: 2D array longitude, latitude matrices at the center + ! bnds_lon, bnds_lat: 2D array longitude, latitude matrices of the vertices + ! rlon, rlat: location of the radar + ! rdist: maximum distance to the radar + ! gatehdist: distance to the radar in the model grid + ! inradar: whether the model grid lays inside the measurement area of any radar + ! radnum: pertinence of the model grid to a given radar + fname = 'radar_inside' + + dbg = set_dbg(idbg) + + ! Is this radar inside the range of the tile? + inside = .FALSE. + + gatehdist = missing_value + inradar = 0 + radnum = fillValueI + + ! Do we need to be efficient? It is not working.... + ii = ims + ie = ime + ji = jms + je = jme + CALL is_inside(ii, ie, ji, je, lon(ii:ie,ji:je), lat(ii:ie,ji:je), rlon, rlat, rdist, inside, dbg) + + ! Going to the borders + DO ibnd=1, 4 + CALL is_inside(ii, ie, ji, je, bnds_lon(ii:ie,ji:je,ibnd), bnds_lat(ii:ie,ji:je,ibnd), rlon, & + rlat, rdist, inside, dbg) + END DO + + !! S border + IF (.NOT.inside) THEN + ii = ims + ie = ime + ji = jms + je = jms + CALL is_inside(ii, ie, ji, je, lon(ii:ie,ji:je), lat(ii:ie,ji:je), rlon, rlat, rdist, inside, dbg) + + ! Going to the borders + DO ibnd=1, 4 + CALL is_inside(ii, ie, ji, je, bnds_lon(ii:ie,ji:je,ibnd), bnds_lat(ii:ie,ji:je,ibnd), rlon, & + rlat, rdist, inside, dbg) + END DO + END IF + + !! W border + IF (.NOT.inside) THEN + ii = ims + ie = ims + ji = jms + je = jme + CALL is_inside(ii, ie, ji, je, lon(ii:ie,ji:je), lat(ii:ie,ji:je), rlon, rlat, rdist, inside, dbg) + + ! Going to the borders + DO ibnd=1, 4 + CALL is_inside(ii, ie, ji, je, bnds_lon(ii:ie,ji:je,ibnd), bnds_lat(ii:ie,ji:je,ibnd), rlon, & + rlat, rdist, inside, dbg) + END DO + END IF + + !! N border + IF (.NOT.inside) THEN + ii = ims + ie = ime + ji = jme + je = jme + CALL is_inside(ii, ie, ji, je, lon(ii:ie,ji:je), lat(ii:ie,ji:je), rlon, rlat, rdist, inside, dbg) + + ! Going to the borders + DO ibnd=1, 4 + CALL is_inside(ii, ie, ji, je, bnds_lon(ii:ie,ji:je,ibnd), bnds_lat(ii:ie,ji:je,ibnd), rlon, & + rlat, rdist, inside, dbg) + END DO + END IF + + !! E border + IF (.NOT.inside) THEN + ii = ime + ie = ime + ji = jms + je = jme + CALL is_inside(ii, ie, ji, je, lon(ii:ie,ji:je), lat(ii:ie,ji:je), rlon, rlat, rdist, inside, dbg) + + ! Going to the borders + DO ibnd=1, 4 + CALL is_inside(ii, ie, ji, je, bnds_lon(ii:ie,ji:je,ibnd), bnds_lat(ii:ie,ji:je,ibnd), rlon, & + rlat, rdist, inside, dbg) + END DO + END IF + + ! Filling other values with a single loop + IF (inside) THEN + CALL real_distances(ims, ime, jms, jme, lon, lat, rlon, rlat, gatehdist) + DO i=ims, ime + DO j=jms, jme + IF (gatehdist(i,j) <= rdist) THEN + inradar(i,j) = 1 + radnum(i,j) = 1 + ELSE + gatehdist(i,j) = missing_value + END IF + END DO + END DO + mask = radnum == 1 + IF (dbg) PRINT *,' ' // TRIM(fname) // ": found ", SUM(radnum, mask=mask), " inside" + END IF + + RETURN + + END SUBROUTINE radar_inside + + SUBROUTINE is_inside(ims, ime, jms, jme, lon, lat, rlon, rlat, rdist, inside, idbg) + ! Subroutine to check if a tile gets a radar + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ims, ime, jms, jme + REAL, DIMENSION(ims:ime,jms:jme), INTENT(in) :: lon, lat + REAL, INTENT(in) :: rlon, rlat, rdist + LOGICAL, OPTIONAL :: idbg + LOGICAL, INTENT(inout) :: inside + + ! Local + REAL :: mindist + REAL, DIMENSION(:,:), ALLOCATABLE :: distborder + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon, lat: 2D array longitude, latitude matrices at the center + ! rlonv, rlat: location of the radar + ! rdist: maximum distance to the radar + ! gatehdist: distance to the radar in the model grid + fname = 'is_inside' + + dbg = set_dbg(idbg) + + IF (ALLOCATED(distborder)) DEALLOCATE(distborder) + ALLOCATE(distborder(ims:ime,jms:jme)) + CALL real_distances(ims, ime, jms, jme, & + lon(ims:ime,jms:jme), lat(ims:ime,jms:jme), & + rlon, rlat, & + distborder) + mindist = MINVAL(distborder) + + IF (mindist < rdist) inside = .TRUE. + + IF (dbg) PRINT *, ' ' //TRIM(fname)//' for ims,ime;jms,jme ',ims,' , ',ime,' ; ',jms,' , ',jme, & + ' min,max lon ', MINVAL(lon), ',', MAXVAL(lon), ' min,max lat ', MINVAL(lat), ',', MAXVAL(lat), & + ' mindist: ', mindist, ' inside ', inside + + IF (ALLOCATED(distborder)) DEALLOCATE(distborder) + + RETURN + + END SUBROUTINE is_inside + + SUBROUTINE beam_height(hgt_pt, a, angle, azimuth, dist, hgt, lhgt) + ! Function to compute the height of the beam + + ! After + ! + ! [PURPOSE:] Common radar procedures and variables. + ! + ! [HISTORY:] + ! 02/19/2014 Juan Ruiz created + ! + + IMPLICIT NONE + + REAL, INTENT(in) :: hgt_pt, angle, azimuth, dist + REAL(kind=rr8), INTENT(in) :: a + REAL, INTENT(out) :: hgt, lhgt + + ! Local + REAL :: tmp1, tmp2 + REAL(kind=rr8) :: height, local_elevation + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! hgt_pt: height of the radar station [m] + ! a: radar ellipsoide semimajor axis equatorial [m] + ! angle: angle of the beam [rad] + ! azimuth: azimuth of the beam [rad] + ! dist: distance respect the radar [m] + fname = 'beam_height' + + ! Initializing + height = 0.d0 + local_elevation = 0.d0 + + tmp1 = dist**2 + (ke*a)**2 + 2.*dist*ke*a*sin(azimuth) + height = hgt_pt + SQRT( tmp1 ) - ke*a + + !Compute the local elevation angle tacking into account the efective earth radius. + tmp1= dist * cos( azimuth ) + tmp2= dist * sin( azimuth ) + ke*a + local_elevation = azimuth + atan(tmp1/tmp2) + + IF (height < 0.) THEN + PRINT *,' ' // TRIM(ErrorMsg) + PRINT *,' dist: ', dist, ' azimuth: ', azimuth + PRINT *,' radar height: ', height + CALL ErrMsg('Negative radar height !', fname, -1) + END IF + + hgt = REAL(height) + lhgt = REAL(local_elevation) + + RETURN + + END SUBROUTINE beam_height + + SUBROUTINE real_distances_ellipsoide(ims, ime, jms, jme, a, f, elips & + ,lon, lat & + ,lonv, latv & + ,dist, idbg) +! ,dlonpt, dlatpt) + ! Subroutine to compute the diferential of coordinattes of a distance from a given point on an + ! ellipsoide + ! Following software by NOAA: + ! https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html + ! + ! Bessel, F. W., (1826) On the computation of geographical longitude and latitude from geodetic + ! measurements, Astronomische Nachrichten (Astronomical Notes), Band 4, No. 86, Spalten 241-254, + ! Altona 1826 + ! Vincenty, T., (1975) Direct and inverse solutions on the ellipsoid with application of nested + ! equations, Survey Review, Vol. 23, No. 176, pp. 88-93. + ! Vincenty, T., (1976) Correspondence: solutions of geodesics, survey Review, Vol. 23, No. 180, + ! p. 294. + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ims, ime, jms, jme + REAL(kind=rr8), INTENT(in) :: a, f + CHARACTER(len=30), INTENT(in) :: elips + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(in) :: lon, lat + REAL, INTENT(in) :: lonv, latv + LOGICAL, OPTIONAL :: idbg + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(out) :: dist +! REAL, INTENT(out) :: dlonpt, dlatpt + + ! Local + INTEGER :: i,j + REAL :: angleab + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon, lat: 2D array longitude, latitude matrices + ! lonv, latv: point to compute the distance + ! a: semimajor axis equatorial (in meters) + ! f: flattening + ! elips: name of the ellipsis being used + + fname = 'real_distance_ellipsoide' + + dbg = set_dbg(idbg) + + dist = 0. + + DO i=ims, ime + DO j=jms, jme + CALL distance_between_points(lonv, latv, lon(i,j), lat(i,j), a, f, elips, dist(i,j), angleab) + END DO + END DO + + RETURN + + END SUBROUTINE real_distances_ellipsoide + + SUBROUTINE real_distances(ims, ime, jms, jme & + ,lon, lat & + ,lonv, latv & + ,dist) +! ,dlonpt, dlatpt) + ! Subroutine to compute the diferential of coordinattes of a distance from a given point + ! Following: https://en.wikipedia.org/wiki/Geographical_distance + ! + ! D = 2*R * SQRT( [sin(dlat /2)*cos(dlon/2)]**2 + [cos(latm)*sin(dlon/2)]**2 ) + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ims, ime, jms, jme + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(in) :: lon, lat + REAL, INTENT(in) :: lonv, latv + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(out) :: dist +! REAL, INTENT(out) :: dlonpt, dlatpt + + ! Local +! REAL, PARAMETER :: R = 6371.009 ! km + REAL, PARAMETER :: R = 6371009. + REAL, DIMENSION(ims:ime,jms:jme) :: dlat, dlon, latm, sec1, sec2 + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! lon, lat: 2D array longitude, latitude matrices + ! lonv, latv: point to compute the distance + fname = 'real_distance' + + dist = 0. + + dlat = (latv - lat)/DegRad + dlon = (lonv - lon)/DegRad + latm = (lat + latv)/DegRad / 2. + + sec1 = sin(dlat /2.)*cos(dlon/2.) + sec2 = cos(latm)*sin(dlon/2.) + dist = 2.*R * SQRT( (sec1)**2 + (sec2)**2 ) + + RETURN + + END SUBROUTINE real_distances + +!!!! basic + + LOGICAL FUNCTION set_dbg(idbg) + ! Function to setup debug value + + IMPLICIT NONE + + LOGICAL, OPTIONAL :: idbg + + set_dbg = .FALSE. + + IF (.NOT.PRESENT(idbg)) THEN + set_dbg = .FALSE. + ELSE + set_dbg = idbg + END IF + + RETURN + + END FUNCTION set_dbg + + SUBROUTINE ErrMsgIAvail(msg, Navail, avail, funcn) + ! Subroutine to stop execution by non exsitent integer value from available ones and provide an + ! error message + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(in) :: msg, funcn + INTEGER, INTENT(in) :: Navail + INTEGER, DIMENSION(Navail), INTENT(in) :: avail + + ! Local + INTEGER :: i + CHARACTER (LEN=1000) :: diag_message + + !!!!!!! Variables + ! msg: message to print with the error + ! funcn: name of the funciton, section to localize the error + ! Navail: amount of available values + ! avail: available values + + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " // TRIM(funcn) // ": " // TRIM(msg) // & + ' not ready !! Available ones: ', (TRIM(ItoS(avail(i))) // ', ', i=1, Navail) + CALL wrf_error_fatal ( diag_message ) + + RETURN + + END SUBROUTINE ErrMsgIAvail + + SUBROUTINE ErrMsgSAvail(msg, Navail, avail, funcn) + ! Subroutine to stop execution by non exsitent string value from available ones and provide an error + ! message + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(in) :: msg, funcn + INTEGER, INTENT(in) :: Navail + CHARACTER(len=*), DIMENSION(Navail), INTENT(in) :: avail + + ! Local + INTEGER :: i + CHARACTER (LEN=1000) :: diag_message + + !!!!!!! Variables + ! msg: message to print with the error + ! funcn: name of the funciton, section to localize the error + ! Navail: amount of available values + ! avail: available values + + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " // TRIM(funcn) // ": " // TRIM(msg) // & + ' not ready !! Available ones: ', ("'" // TRIM(avail(i)) // "', ", i=1, Navail) + CALL wrf_error_fatal ( diag_message ) + + RETURN + + END SUBROUTINE ErrMsgSAvail + + SUBROUTINE ErrMsgAlloc(msg, arrayn, rank, shapev, funcn, errN) + ! Subroutine to stop execution and provide an error message when allocation problems + + IMPLICIT NONE + + INTEGER, INTENT(in) :: errN, rank + CHARACTER(LEN=*), INTENT(in) :: msg, funcn, arrayn + INTEGER, DIMENSION(rank), INTENT(in) :: shapev + + ! Local + INTEGER :: i + CHARACTER (LEN=1000) :: diag_message + + !!!!!!! Variables + ! msg: message to print with the error + ! rank: rank of the array to allocate + ! shapev: shape of the array to allocate + ! funcn: name of the funciton, section to localize the error + ! errN: number of the error provided for a given FORTRAN function + + IF (errN /= 0) THEN + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " // TRIM(funcn) // ": " // TRIM(msg) // & + ' Problem allocating : ' // "'" // TRIM(arrayn) // "' (", (TRIM(ItoS(shapev(i))) // ', ',i=1, & + rank), " )' !!" + CALL wrf_error_fatal ( diag_message ) + END IF + + + RETURN + + END SUBROUTINE ErrMsgAlloc + + SUBROUTINE ErrMsg(msg, funcn, errN) + ! Subroutine to stop execution and provide an error message + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(in) :: msg, funcn + INTEGER, INTENT(in) :: errN + + ! Local + CHARACTER (LEN=1000) :: diag_message + CHARACTER(len=S32) :: emsg + + !!!!!!! Variables + ! msg: message to print with the error + ! funcn: name of the funciton, section to localize the error + ! errN: number of the error provided for a given FORTRAN function + + IF (errN /= 0) THEN + PRINT *,TRiM(emsg) + PRINT *,' ' // TRIM(funcn) // ': ' // TRIM(msg) // ' !!' + PRINT *,' error number:', errN + STOP + + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " // TRIM(funcn) // ": " // TRIM(msg) // ' !!' + CALL wrf_error_fatal ( diag_message ) + END IF + + RETURN + + END SUBROUTINE ErrMsg + + CHARACTER(len=S32)FUNCTION DtoS(Dval) + ! DtoS: Function to transform a double real to String + + IMPLICIT NONE + + REAL(kind=rr8), INTENT(IN) :: Dval + + WRITE(DtoS,'(F21.10)')Dval + + END FUNCTION DtoS + + CHARACTER(len=S32)FUNCTION RtoS(Rval) + ! RtoS: Function to transform a real to String + + IMPLICIT NONE + + REAL, INTENT(IN) :: Rval + + WRITE(RtoS,'(F21.10)')Rval + + END FUNCTION RtoS + + CHARACTER(len=S32)FUNCTION ItoS(Ival) + ! ItoS: Function to transform an integer to String + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: Ival + + ! Local + CHARACTER(len=S32) :: itoS0 + + WRITE(ItoS0,'(I50)')Ival + CALL removeNONnum(ItoS0, ItoS) + + END FUNCTION ItoS + + INTEGER FUNCTION StoI(String) + ! Function to transform a String to an integer + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(IN) :: String + ! Local + INTEGER :: Lstring + + Lstring = LEN_TRIM(String) + + READ(String,'(I' // TRIM(ItoS(Lstring)) // ')')StoI + + END FUNCTION StoI + + REAL FUNCTION StoR(String) + ! Function to transform a String to a real + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(IN) :: String + + ! Local + INTEGER :: Lstring + + Lstring = LEN_TRIM(String) + + READ(String,'(F' // TRIM(ItoS(Lstring)) // '.0)')StoR + + END FUNCTION StoR + + LOGICAL FUNCTION StoL(String) + ! Function to transform a String to a boolean + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(IN) :: String + + IF (TRIM(String) == '.T.' .OR. TRIM(String) == '.true.' .OR. TRIM(String) == '.TRUE.' & + .OR. TRIM(String) == 'yes' .OR. TRIM(String) == 'YES' .OR. TRIM(String) == 'y' ) THEN + StoL = .TRUE. + ELSE + StoL = .FALSE. + END IF + + END FUNCTION StoL + + LOGICAL FUNCTION FileExist(filen) + ! Checks if a file exists + + IMPLICIT NONE + + CHARACTER(len=*), INTENT(in) :: filen + + FileExist = .FALSE. + INQUIRE(file=TRIM(filen), exist=FileExist) + + RETURN + + END FUNCTION FileExist + + INTEGER FUNCTION freeunit() + ! provides the number of a free unit in which open a file + + IMPLICIT NONE + + LOGICAL :: is_used + + is_used = .true. + DO freeunit=10,100 + INQUIRE(unit=freeunit, opened=is_used) + IF (.not. is_used) EXIT + END DO + + RETURN + + END FUNCTION freeunit + + SUBROUTINE removeNONnum(String, newString) + ! Subroutine to remove non numeric characters from a string + + IMPLICIT NONE + + CHARACTER(LEN=*), INTENT(IN) :: String + CHARACTER(LEN=*), INTENT(OUT) :: newString + + ! Local + INTEGER :: ic, inc, Lstring + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! String: string to remove non-numeric characters + ! newString: resultant string + fname = 'removeNONnum' + + Lstring = LEN_TRIM(String) + newString = '' + inc = 1 + DO ic=1, Lstring + IF (ICHAR(String(ic:ic)) >= ICHAR('0') .AND. ICHAR(String(ic:ic)) <= ICHAR('9')) THEN + newString(inc:inc) = String(ic:ic) + inc = inc + 1 + END IF + END DO + + END SUBROUTINE removeNONnum + +!!!! netCDF + + SUBROUTINE open_NCfile(filename, ncid, idbg) + ! Subroutine to open a basic netCDF file + + USE netcdf + + IMPLICIT NONE + + INCLUDE 'netcdf.inc' + + CHARACTER(len=S128), INTENT(in) :: filename + LOGICAL, OPTIONAL :: idbg + INTEGER, INTENT(out) :: ncid + + ! Local + INTEGER :: rcode + LOGICAL :: fexist + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! filename: name of the file to create + + fname = 'open_NCfile' + + dbg = set_dbg(idbg) + + ncid = -1 + + fexist = FileExist(filename) + IF (.NOT.fexist) THEN + CALL ErrMsg("netCDF file '" // TRIM(filename) // "' does not exist !!", fname, -1) + END IF + + rcode = NF90_OPEN(filename, mode=NF90_NOWRITE, ncid=ncid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + RETURN + + END SUBROUTINE open_NCfile + + SUBROUTINE get_Ndims_ncid(ncid, Ndims) + ! Subroutine to get the amount of dimensions from an open netCDF file + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid + INTEGER, INTENT(out) :: Ndims + + ! Local + INTEGER :: idd + INTEGER :: rcode + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! ncid: netCDF open object + fname = 'get_Ndims_ncid' + + rcode = nf90_inquire(ncid, nDimensions=Ndims) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + RETURN + + END SUBROUTINE get_Ndims_ncid + + SUBROUTINE get_Nvars_ncid(ncid, Nvars) + ! Subroutine to get the amount of variables from an open netCDF file + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid + INTEGER, INTENT(out) :: Nvars + + ! Local + INTEGER :: idd + INTEGER :: rcode + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! ncid: netCDF open object + fname = 'get_Nvars_ncid' + + rcode = nf90_inquire(ncid, nVariables=Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + RETURN + + END SUBROUTINE get_Nvars_ncid + + LOGICAL FUNCTION isin_ncunit(nid, varname) + ! Function to tell if a given variable is inside a netcdf file unit + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: nid + CHARACTER(LEN=*), INTENT(in) :: varname + + ! Local + INTEGER :: vid, Ndims, Nvars + INTEGER :: iv, rcode + CHARACTER(LEN=1000) :: varinfile + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! nid: number of the open netCDF + ! varname: name of the variable + + fname = 'isin_ncunit' + + ! Initializing + isin_ncunit = .FALSE. + + rcode = nf90_inquire(nid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + DO iv=1, Nvars + rcode = nf90_inquire_variable(nid, iv, name=varinfile) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (TRIM(varinfile) == TRIM(varname)) THEN + isin_ncunit = .TRUE. + EXIT + ELSE + isin_ncunit = .FALSE. + END IF + END DO + + RETURN + + END FUNCTION isin_ncunit + + LOGICAL FUNCTION isin_dim_ncid(nid, dimname) + ! Function to tell if a given dimension is inside a netcdf file unit + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: nid + CHARACTER(LEN=*), INTENT(in) :: dimname + + ! Local + INTEGER :: vid, Ndims, Nvars + INTEGER :: iv, id, rcode + CHARACTER(len=S16) :: diminfile + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! nid: number of the open netCDF + ! dimname: name of the dimension + + fname = 'isin_dim_ncid' + + ! Initializing + isin_dim_ncid = .FALSE. + + rcode = nf90_inquire(nid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + DO id=1, Ndims + rcode = nf90_inquire_dimension(nid, id, name=diminfile) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (TRIM(diminfile) == TRIM(dimname)) THEN + isin_dim_ncid = .TRUE. + EXIT + END IF + END DO + + RETURN + + END FUNCTION isin_dim_ncid + + SUBROUTINE get_dimlength_ncid(ncid, dimn, dimlength) + ! Subroutine to get the length of a dimension from an existing to an open netCDF + + USE netcdf + + IMPLICIT NONE + + INCLUDE 'netcdf.inc' + + INTEGER, INTENT(in) :: ncid + CHARACTER(len=*), INTENT(in) :: dimn + INTEGER, INTENT(out) :: dimlength + + ! Local + INTEGER :: idd + INTEGER :: Ndims + INTEGER, DIMENSION(:), ALLOCATABLE :: dimids + CHARACTER(len=S32), DIMENSION(:), ALLOCATABLE :: dimns + INTEGER :: rcode + LOGICAL :: infile + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! ncid: netCDF open object + ! dimn: name of dimension + fname = 'get_dimlength_ncid' + + dimlength = -1 + + infile = isin_dim_ncid(ncid, dimn) + + IF (.NOT.infile) THEN + PRINT *, Errormsg + CALL get_Ndims_ncid(ncid, Ndims) + + IF (ALLOCATED(dimns)) DEALLOCATE(dimns) + ALLOCATE(dimns(Ndims)) + CALL get_dims_ncid(ncid, Ndims, dimns) + + CALL ErrMsgSAvail("ncid does not contain dimension named '" // TRIM(dimn) // "' !!", Ndims, & + dimns, fname) + END IF + + rcode = nf90_inq_dimid(ncid, dimn, idd) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + rcode = nf90_inquire_dimension(ncid, idd, len=dimlength) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + RETURN + + END SUBROUTINE get_dimlength_ncid + + SUBROUTINE get_dims_ncid(ncid, Ndims, dimns, dimlengths) + ! Subroutine to get the dimensions from an open netCDF file + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid, Ndims + CHARACTER(len=S32), DIMENSION(Ndims), INTENT(out) :: dimns + INTEGER, DIMENSION(Ndims), OPTIONAL :: dimlengths + + ! Local + INTEGER :: idd + INTEGER :: rcode + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! ncid: netCDF open object + ! Ndims: amount of fimensions in the file + fname = 'get_dims_ncid' + + dimns = '' + IF (PRESENT(dimlengths)) dimlengths = -1 + + DO idd=1, Ndims + IF (PRESENT(dimlengths)) THEN + rcode = nf90_inquire_dimension(ncid, idd, name=dimns(idd), len=dimlengths(idd)) + ELSE + rcode = nf90_inquire_dimension(ncid, idd, name=dimns(idd)) + END IF + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (PRESENT(dimlengths)) THEN + IF (dimlengths(idd) == NF90_UNLIMITED) PRINT *,' ' // TRIM(fname) // ": dimension '" // & + TRIM(dimns(idd)) // "' is unlimited ", NF90_UNLIMITED + END IF + END DO + + RETURN + + END SUBROUTINE get_dims_ncid + + SUBROUTINE get_varS0D_ncid(ncid, StrL, vname, vals, idbg) + ! Subroutine to get a 0D String variable from a netCDF file unit passing starting and ending + ! First dimension accounts for length of the string values + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid, StrL + CHARACTER(LEN=S32), INTENT(in) :: vname + LOGICAL, OPTIONAL :: idbg + CHARACTER(len=1), DIMENSION(StrL), INTENT(out) :: vals + + ! Local + INTEGER :: i, Ndims, Nvars + INTEGER :: rcode, varid, vartype, rank + LOGICAL :: vfound + INTEGER, DIMENSION(:), ALLOCATABLE :: invarids + CHARACTER(len=S32), DIMENSION(:), ALLOCATABLE :: invarns + LOGICAL :: infile + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! ncid: netCDF file identifier + ! StrL: length of the String values + ! vals: values to get + ! vname: name of the variable to get + fname = 'get_varS0D_ncid' + + dbg = set_dbg(idbg) + + ! Checking existence of variable + infile = isin_ncunit(ncid, vname) + IF (.NOT.infile) THEN + rcode = nf90_inquire(ncid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + ALLOCATE(invarids(Nvars)) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + ALLOCATE(invarns(Nvars)) + + CALL get_vars_ncid(ncid, Nvars, invarids, invarns) + PRINT *,' Available variales in netcdf unit: ', ("'" // TRIM(invarns(i)) // "', ", i=1, Nvars) + CALL ErrMsg("Passed variable '" // TRIM(vname) // "' not found", fname, -1) + END IF + + rcode = nf90_inq_varid(ncid, vname, varid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + ! Checking variable + rcode = nf90_inquire_variable(ncid, varid, xtype=vartype, ndims=rank) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + ! First dimension accounts for length of the string values + rank = rank - 1 + IF (vartype /= NF90_CHAR) CALL ErrMSg("Wrong type of variable '" // TRIM(vname) // & + "'. I expected NF90_CHAR= " // TRIM(ItoS(NF90_DOUBLE)) // ' instead got ' // & + TRIM(ItoS(vartype)), fname, -1) + IF (rank /= 0) CALL ErrMSg("Wrong rank of variable '" //TRIM(vname)// "'. I expected rank= 0 " // & + ' instead got ' // TRIM(ItoS(rank)), fname, -1) + + rcode = nf90_get_var(ncid, varid, vals) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + + RETURN + + END SUBROUTINE get_varS0D_ncid + + SUBROUTINE get_varR0D_ncid(ncid, vname, vals, idbg) + ! Subroutine to get a 0D Real variable from a netCDF file unit passing starting and ending + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid + CHARACTER(LEN=S32), INTENT(in) :: vname + LOGICAL, OPTIONAL :: idbg + REAL, INTENT(out) :: vals + + ! Local + INTEGER :: i, Ndims, Nvars + INTEGER :: rcode, varid, vartype, rank + LOGICAL :: vfound + INTEGER, DIMENSION(:), ALLOCATABLE :: invarids + CHARACTER(len=S32), DIMENSION(:), ALLOCATABLE :: invarns + LOGICAL :: infile + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! ncid: netCDF file identifier + ! vals: values to get + ! vname: name of the variable to get + fname = 'get_varR0D_ncid' + + dbg = set_dbg(idbg) + + ! Checking existence of variable + infile = isin_ncunit(ncid, vname) + IF (.NOT.infile) THEN + rcode = nf90_inquire(ncid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + ALLOCATE(invarids(Nvars)) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + ALLOCATE(invarns(Nvars)) + + CALL get_vars_ncid(ncid, Nvars, invarids, invarns) + PRINT *,' Available variales in netcdf unit: ', ("'" // TRIM(invarns(i)) // "', ", i=1, Nvars) + CALL ErrMsg("Passed variable '" // TRIM(vname) // "' not found", fname, -1) + END IF + + rcode = nf90_inq_varid(ncid, vname, varid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + ! Checking variable + rcode = nf90_inquire_variable(ncid, varid, xtype=vartype, ndims=rank) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (vartype /= NF90_FLOAT) CALL ErrMSg("Wrong type of variable '" // TRIM(vname) // & + "'. I expected NF90_FLOAT= " // TRIM(ItoS(NF90_FLOAT)) // ' instead got ' // & + TRIM(ItoS(vartype)), fname, -1) + IF (rank /= 0) CALL ErrMSg("Wrong rank of variable '" //TRIM(vname)// "'. I expected rank= 0 " // & + ' instead got ' // TRIM(ItoS(rank)), fname, -1) + + rcode = nf90_get_var(ncid, varid, vals) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + + RETURN + + END SUBROUTINE get_varR0D_ncid + + SUBROUTINE get_varR1D_ncid(ncid, id1, ed1, vname, vals, idbg) + ! Subroutine to get a 1D Real variable from a netCDF file unit passing starting and ending + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid, id1, ed1 + CHARACTER(LEN=S32), INTENT(in) :: vname + LOGICAL, OPTIONAL :: idbg + REAL, DIMENSION(id1:ed1),INTENT(out) :: vals + + ! Local + INTEGER :: i, Ndims, Nvars + INTEGER :: rcode, varid, vartype, rank + INTEGER :: d1 + LOGICAL :: vfound + INTEGER, DIMENSION(:), ALLOCATABLE :: invarids + CHARACTER(len=S32), DIMENSION(:), ALLOCATABLE :: invarns + LOGICAL :: infile + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! ncid: netCDF file identifier + ! id1,ed1: shape of the matrix + ! vals: values to get + ! vname: name of the variable to get + fname = 'get_varR1D_ncid' + + dbg = set_dbg(idbg) + + ! Checking existence of variable + infile = isin_ncunit(ncid, vname) + IF (.NOT.infile) THEN + rcode = nf90_inquire(ncid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + ALLOCATE(invarids(Nvars)) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + ALLOCATE(invarns(Nvars)) + + CALL get_vars_ncid(ncid, Nvars, invarids, invarns) + PRINT *,' Available variales in netcdf unit: ', ("'" // TRIM(invarns(i)) // "', ", i=1, Nvars) + CALL ErrMsg("Passed variable '" // TRIM(vname) // "' not found", fname, -1) + END IF + + rcode = nf90_inq_varid(ncid, vname, varid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + ! Checking variable + rcode = nf90_inquire_variable(ncid, varid, xtype=vartype, ndims=rank) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (vartype /= NF90_FLOAT) CALL ErrMSg("Wrong type of variable '" // TRIM(vname) // & + "'. I expected NF90_FLOAT= " // TRIM(ItoS(NF90_FLOAT)) // ' instead got ' // & + TRIM(ItoS(vartype)), fname, -1) + IF (rank /= 1) CALL ErrMSg("Wrong rank of variable '" //TRIM(vname)// "'. I expected rank= 1 " // & + ' instead got ' // TRIM(ItoS(rank)), fname, -1) + + d1 = ed1-id1+1 + rcode = nf90_get_var(ncid, varid, vals, start=(/id1/), count=(/d1/)) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + + RETURN + + END SUBROUTINE get_varR1D_ncid + + SUBROUTINE get_varD3D_ncid(ncid, id1, id2, id3, ed1, ed2, ed3, vname, vals, idbg) + ! Subroutine to get a 3D double Real variable from a netCDF file unit passing starting and ending + + USE netcdf + + IMPLICIT NONE + + INTEGER, INTENT(in) :: ncid, id1, ed1, id2, ed2, id3, ed3 + CHARACTER(LEN=S32), INTENT(in) :: vname + LOGICAL, OPTIONAL :: idbg + REAL(kind=rr8), DIMENSION(id1:ed1,id2:ed2,id3:ed3), & + INTENT(out) :: vals + + ! Local + INTEGER :: i, Ndims, Nvars + INTEGER :: rcode, varid, vartype, rank + INTEGER :: d1, d2, d3 + LOGICAL :: vfound + INTEGER, DIMENSION(:), ALLOCATABLE :: invarids + CHARACTER(len=S32), DIMENSION(:), ALLOCATABLE :: invarns + LOGICAL :: infile + LOGICAL :: dbg + CHARACTER(len=S32) :: fname + + !!!!!!! Variables + ! ncid: netCDF file identifier + ! id1,ed1,id2,ed2,id3,ed3: shape of the matrix + ! vals: values to get + ! vname: name of the variable to get + fname = 'get_varD3D_ncid' + + dbg = set_dbg(idbg) + + ! Checking existence of variable + infile = isin_ncunit(ncid, vname) + IF (.NOT.infile) THEN + rcode = nf90_inquire(ncid, Ndims, Nvars) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + ALLOCATE(invarids(Nvars)) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + ALLOCATE(invarns(Nvars)) + + CALL get_vars_ncid(ncid, Nvars, invarids, invarns) + PRINT *,' Available variales in netcdf unit: ', ("'" // TRIM(invarns(i)) // "', ", i=1, Nvars) + CALL ErrMsg("Passed variable '" // TRIM(vname) // "' not found", fname, -1) + END IF + + rcode = nf90_inq_varid(ncid, vname, varid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + ! Checking variable + rcode = nf90_inquire_variable(ncid, varid, xtype=vartype, ndims=rank) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + IF (vartype /= NF90_DOUBLE) CALL ErrMSg("Wrong type of variable '" // TRIM(vname) // & + "'. I expected NF90_DOUBLE= " // TRIM(ItoS(NF90_DOUBLE)) // ' instead got ' // & + TRIM(ItoS(vartype)), fname, -1) + IF (rank /= 3) CALL ErrMSg("Wrong rank of variable '" //TRIM(vname)// "'. I expected rank= 3 " // & + ' instead got ' // TRIM(ItoS(rank)), fname, -1) + + d1 = ed1-id1+1 + d2 = ed2-id2+1 + d3 = ed3-id3+1 + rcode = nf90_get_var(ncid, varid, vals, start=(/id1,id2,id3/), count=(/d1,d2,d3/)) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + IF (ALLOCATED(invarids)) DEALLOCATE(invarids) + IF (ALLOCATED(invarns)) DEALLOCATE(invarns) + + RETURN + + END SUBROUTINE get_varD3D_ncid + + ! From UNIDATA: https://www.unidata.ucar.edu/software/netcdf/docs/netcdf-f90.html + SUBROUTINE handle_errf(st, fn) + ! Subroutine to provide the error message when something with netCDF went wrong (including fname) + + USE netcdf + + INTEGER, INTENT(in) :: st + CHARACTER(len=*), INTENT(in) :: fn + CHARACTER (LEN=1000) :: diag_message + + !!!!!!! Variables + ! st: netCDF status number + ! fn: function name from which it is used + + IF (st /= nf90_noerr) THEN + WRITE (diag_message , * ) & + TRIM(errormsg) // ": " // TRIM(fn) // ": " // TRIM(nf90_strerror(st)) + CALL wrf_error_fatal ( diag_message ) + + END IF + + END SUBROUTINE handle_errf + + SUBROUTINE close_NCfile(filename, ncid, idbg) + ! Subroutine to close a netCDF file + + USE netcdf + + IMPLICIT NONE + + INCLUDE 'netcdf.inc' + + CHARACTER(len=S128), INTENT(in) :: filename + LOGICAL, OPTIONAL :: idbg + INTEGER, INTENT(out) :: ncid + + ! Local + INTEGER :: rcode + LOGICAL :: fexist + LOGICAL :: dbg + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! filename: name of the file to create + + fname = 'close_NCfile' + + dbg = set_dbg(idbg) + + rcode = NF90_CLOSE(ncid) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + + RETURN + + END SUBROUTINE close_NCfile + + SUBROUTINE get_vars_ncid(ncid, Nvars, invarids, invarns) + ! Subroutine to provide all the ids and names of the variables existing to an open netCDF + ! previously knowing how many variables the file has via Nvars + ! rcode = nf90_inq_nvars(ncid, Nvars) + + USE netcdf + + IMPLICIT NONE + + INCLUDE 'netcdf.inc' + + INTEGER, INTENT(in) :: ncid, Nvars + INTEGER, DIMENSION(Nvars), INTENT(out) :: invarids + CHARACTER(len=S32), DIMENSION(Nvars), INTENT(out) :: invarns + + ! Local + INTEGER :: iv + INTEGER :: rcode + CHARACTER(len=Sl) :: fname + + !!!!!!! Variables + ! ncid: netCDF open object + ! Nvars: amount of variables + ! invarids: Ids of the variables + ! invarns: name of the variables + fname = 'get_vars_ncid' + + invarids = -1 + invarns = '' + + DO iv=1, Nvars + invarids(iv) = iv + rcode = nf90_inquire_variable(ncid, invarids(iv), invarns(iv)) + IF (rcode /= NF90_NOERR) CALL handle_errf(rcode, fname) + END DO + + RETURN + + END SUBROUTINE get_vars_ncid + + +#endif +END MODULE module_simradar_tools diff --git a/share/output_wrf.F b/share/output_wrf.F index 69ebcf31fc..e2f9176f87 100644 --- a/share/output_wrf.F +++ b/share/output_wrf.F @@ -1000,6 +1000,30 @@ SUBROUTINE output_wrf ( fid , grid , config_flags, switch , ierr ) CALL wrf_put_dom_ti_integer ( fid , 'COMPRESSION' , config_flags%compression , 1 , ierr ) ENDIF +#ifdef SIMRADAR +! No character vectors as attribute +! CALL wrf_put_dom_ti_integer ( fid , 'SIMRADAR_RADAR_NAMES' , model_config_rec%radar_name , 1 , ierr ) +!rconfig character radar_name namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_name" "name of the radar" + + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADAR_CLON' , model_config_rec%radar_clon(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADAR_CLAT' , model_config_rec%radar_clat(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADAR_CHGT' , model_config_rec%radar_chgt(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) +! No character vectors as attribute +!rconfig character radar_band namelist,simradar max_nradar (/'-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-','-'/) rh{11} "radar_band" "band of the radar: 'c', 's', or 'x'" "" + CALL wrf_put_dom_ti_integer ( fid , 'SIMRADAR_RADAR_SCAN' , model_config_rec%radar_scan(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADAR_DISTMAX' , model_config_rec%radar_distmax(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADAR_HGTTMAX' , model_config_rec%radar_hgtmax(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_integer ( fid , 'SIMRADAR_EM_RAY' , model_config_rec%em_ray(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADARDS' , model_config_rec%radards(1:model_config_rec%nradar) , model_config_rec%nradar , ierr ) + CALL wrf_put_dom_ti_integer ( fid , 'SIMRADAR_RADARINTERP' , model_config_rec%radinterp , 1 , ierr ) + CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_RADII_DIST' , model_config_rec%radii_dist , 1 , ierr ) + CALL wrf_put_dom_ti_integer ( fid , 'SIMRADAR_ELOPTION' , model_config_rec%eloption , 1 , ierr ) +! CALL wrf_put_dom_ti_char ( fid , 'SIMRADAR_ELLIPSE_NAME' , TRIM(model_config_rec%elips_name) , ierr ) +! CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_ELLIPSE_A' , model_config_rec%elips_a , 1 , ierr ) +! CALL wrf_put_dom_ti_real ( fid , 'SIMRADAR_ELLIPSE_F' , model_config_rec%elips_f , 1 , ierr ) + CALL wrf_put_dom_ti_char ( fid , 'SIMRADAR_TMATRIX-LUT_FILEN' , TRIM(model_config_rec%lutfilen) , ierr ) +#endif + #if ( (EM_CORE == 1) && (DA_CORE != 1) ) save_topo_orig = grid%save_topo_from_real