diff --git a/Registry/registry.var_chem b/Registry/registry.var_chem index 295a075a21..7415b18a1b 100644 --- a/Registry/registry.var_chem +++ b/Registry/registry.var_chem @@ -124,14 +124,96 @@ state real oc_a04 ijkftb chem 1 - i0{12}rhusdf=(bdy_ state real bc_a04 ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bc_a04" "Black carbon, aerosol bin 04" "ug/kg-dryair" #state real hysw_a04 ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "hysw_a04" "Hysteresis water, aerosol bin 04" "ug/kg-dryair" # + +## aerosol variables for RACM_SOA_VBS +state real so4aj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "so4aj" "Sulfate conc. Acc. mode" "ug/kg-dryair" +state real so4ai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "so4ai" "Sulfate conc. Aitken mode" "ug/kg-dryair" +state real nh4aj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nh4aj" "Ammonium conc. Acc. mode" "ug/kg-dryair" +state real nh4ai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nh4ai" "Ammonium conc. Aitken mode" "ug/kg-dryair" +state real no3aj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "no3aj" "Nitrate conc. Acc. mode" "ug/kg-dryair" +state real no3ai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "no3ai" "Nitrate conc. Aitken mode" "ug/kg-dryair" +state real naaj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "naaj" "Sodium conc. Acc. mode" "ug/kg-dryair" +state real naai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "naai" "Sodium conc. Aitken mode" "ug/kg-dryair" +state real claj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "claj" "Chloride conc. Acc. mode" "ug/kg-dryair" +state real clai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "clai" "Chloride conc. Aitken mode" "ug/kg-dryair" +state real orgpaj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgpaj" "Prim. anth. org. conc. from aromatics Acc. mode" "ug/kg-dryair" +state real orgpai ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgpai" "Prim. anth. org. conc. from aromatics Aitken mode" "ug/kg-dryair" +state real ecj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "ecj" "Elemental carbon Acc. mode" "ug/kg-dryair" +state real eci ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "eci" "Elemental carbon Aitken mode" "ug/kg-dryair" +state real p25j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "p25j" "Primary PM2.5 Acc. mode" "ug/kg-dryair" +state real p25i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "p25i" "Primary PM2.5 Aitken mode" "ug/kg-dryair" +state real antha ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "antha" "Coarse anthropogenic aerosols" "ug/kg-dryair" +state real seas ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "seas" "Coarse marine aerosols" "ug/kg-dryair" +state real soila ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "soila" "Coarse soil-derived aerosols" "ug/kg-dryair" +state real nu0 ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nu0" "Aitken mode number" "#/kg-dryair" +state real ac0 ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "ac0" "Accumulation mode number" "#/kg-dryair" +state real corn ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "corn" "Coarse mode number" "#/kg-dryair" +state real asoa1j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa1j" "SOA Anth. org. conc. with C* value=1 ug/m^3 Acc. mode" "ug/kg-dryair" +state real asoa1i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa1i" "SOA Anth. org. conc. with C* value=1 ug/m^3 Aitken mode" "ug/kg-dryair" +state real asoa2j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa2j" "SOA Anth. org. conc. with C* value=10 ug/m^3 Acc. mode" "ug/kg-dryair" +state real asoa2i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa2i" "SOA Anth. org. conc. with C* value=10 ug/m^3 Aitken mode" "ug/kg-dryair" +state real asoa3j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa3j" "SOA Anth. org. conc. with C* value=100 ug/m^3 Acc. mode" "ug/kg-dryair" +state real asoa3i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa3i" "SOA Anth. org. conc. with C* value=100 ug/m^3 Aitken mode" "ug/kg-dryair" +state real asoa4j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa4j" "SOA Anth. org. conc. with C* value=1000 ug/m^3 Acc. mode" "ug/kg-dryair" +state real asoa4i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "asoa4i" "SOA Anth. org. conc. with C* value=1000 ug/m^3 Aitken mode" "ug/kg-dryair" +state real bsoa1j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa1j" "SOA Biog. org. conc. with C* value=1 ug/m^3 Acc. mode" "ug/kg-dryair" +state real bsoa1i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa1i" "SOA Biog. org. conc. with C* value=1 ug/m^3 Aitken mode" "ug/kg-dryair" +state real bsoa2j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa2j" "SOA Biog. org. conc. with C* value=10 ug/m^3 Acc. mode" "ug/kg-dryair" +state real bsoa2i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa2i" "SOA Biog. org. conc. with C* value=10 ug/m^3 Aitken mode" "ug/kg-dryair" +state real bsoa3j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa3j" "SOA Biog. org. conc. with C* value=100 ug/m^3 Acc. mode" "ug/kg-dryair" +state real bsoa3i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa3i" "SOA Biog. org. conc. with C* value=100 ug/m^3 Aitken mode" "ug/kg-dryair" +state real bsoa4j ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa4j" "SOA Biog. org. conc. with C* value=1000 ug/m^3 Acc. mode" "ug/kg-dryair" +state real bsoa4i ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "bsoa4i" "SOA Biog. org. conc. with C* value=1000 ug/m^3 Aitken mode" "ug/kg-dryair" +# +##cloud-phase aerosol +#state real so4cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "so4cwj" "Sulfate conc. Acc. mode in cloud" "ug/kg-dryair" +#state real so4cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "so4cwi" "Sulfate conc. Aitken mode in cloud" "ug/kg-dryair" +#state real nh4cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nh4cwj" "Ammonium conc. Acc. mode in cloud" "ug/kg-dryair" +#state real nh4cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nh4cwi" "Ammonium conc. Aitken mode in cloud" "ug/kg-dryair" +#state real no3cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "no3cwj" "Nitrate conc. Acc. mode in cloud" "ug/kg-dryair" +#state real no3cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "no3cwi" "Nitrate conc. Aitken mode in cloud" "ug/kg-dryair" +#state real nacwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nacwj" "Sodium conc. Acc. mode in cloud" "ug/kg-dryair" +#state real nacwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nacwi" "Sodium conc. Aitken mode in cloud" "ug/kg-dryair" +#state real clcwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "clcwj" "Chloride conc. Acc. mode in cloud" "ug/kg-dryair" +#state real clcwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "clcwi" "Chloride conc. Aitken mode in cloud" "ug/kg-dryair" +#state real orgaro1cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgaro1cwj" "SOA Anth. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgaro1cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgaro1cwi" "SOA Anth. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgaro2cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgaro2cwj" "SOA Anth. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgaro2cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgaro2cwi" "SOA Anth. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgalk1cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgalk1cwj" "SOA Anth. org. conc. from alkanes and others except aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgalk1cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgalk1cwi" "SOA Anth. org. conc. from alkanes and others except aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgole1cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgole1cwj" "SOA Anth. org. conc. from alkenes and others except aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgole1cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgole1cwi" "SOA Anth. org. conc. from alkenes and others except aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgba1cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba1cwj" "SOA Biog. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgba1cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba1cwi" "SOA Biog. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgba2cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba2cwj" "SOA Biog. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgba2cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba2cwi" "SOA Biog. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgba3cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba3cwj" "SOA Biog. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgba3cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba3cwi" "SOA Biog. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgba4cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba4cwj" "SOA Biog. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgba4cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgba4cwi" "SOA Biog. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real orgpacwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgpacwj" "Prim. anth. org. conc. from aromatics Acc. mode in cloud" "ug/kg-dryair" +#state real orgpacwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "orgpacwi" "Prim. anth. org. conc. from aromatics Aitken mode in cloud" "ug/kg-dryair" +#state real eccwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "eccwj" "Elemental carbon Acc. mode in cloud" "ug/kg-dryair" +#state real eccwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "eccwi" "Elemental carbon Aitken mode in cloud" "ug/kg-dryair" +#state real p25cwj ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "p25cwj" "Primary PM2.5 Acc. mode in cloud" "ug/kg-dryair" +#state real p25cwi ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "p25cwi" "Primary PM2.5 Aitken mode in cloud" "ug/kg-dryair" +#state real anthcw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "anthcw" "Coarse anthropogenic aerosols in cloud" "ug/kg-dryair" +#state real seascw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "seascw" "Coarse marine aerosols in cloud" "ug/kg-dryair" +#state real soilcw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "soilcw" "Coarse soil-derived aerosols in cloud" "ug/kg-dryair" +#state real nu0cw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "nu0cw" "Aitken mode number in cloud" "#/kg-dryair" +#state real ac0cw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "ac0cw" "Accumulation mode number in cloud" "#/kg-dryair" +#state real corncw ijkftb chem 1 - i0{12}rhusdf=(bdy_interp:dt) "corncw" "Coarse mode number in cloud" "#/kg-dryair" +# + rconfig integer io_style_emissions namelist,chem 1 1 r "io_style_emissions" "" "" # Namelist variables specific to DA-Chem -rconfig integer chem_cv_options namelist,wrfvarchem 1 0 - "chem_cv_options" "0:off, 10:gocart, 20:mosaic " "" +rconfig integer chem_cv_options namelist,wrfvarchem 1 0 - "chem_cv_options" "0:off, 10:gocart, 20:mosaic 108:racm_soa_vbs " "" rconfig integer chemicda_opt namelist,wrfvarchem 1 1 - "chemicda_opt" " " " " -rconfig real var_scaling12 namelist,wrfvarchem 36*max_outer_iterations 1.0 - "var_scaling12" "" "" -rconfig real len_scaling12 namelist,wrfvarchem 36*max_outer_iterations 1.0 - "len_scaling12" "" "" -rconfig real max_vert_var12 namelist,wrfvarchem 36*1 99.0 - "max_vert_var12" "" "" +rconfig real var_scaling12 namelist,wrfvarchem 45*max_outer_iterations 1.0 - "var_scaling12" "" "" +rconfig real len_scaling12 namelist,wrfvarchem 45*max_outer_iterations 1.0 - "len_scaling12" "" "" +rconfig real max_vert_var12 namelist,wrfvarchem 45*1 99.0 - "max_vert_var12" "" "" rconfig logical test_transformsch namelist,wrfvarchem 1 .false. - "test_transformsch" "" "" rconfig logical use_chemic_surfobs namelist,wrfvarchem 1 .true. - "use_chemic_surfobs" "" "" rconfig real max_error_chemic_surf namelist,wrfvarchem 1 500.0 - "max_error_chemic_surf" "" "" @@ -189,6 +271,44 @@ typedef xchem_type real chem_ic_oin_a01 ijkf chem_ic 1 - typedef xchem_type real chem_ic_oin_a02 ijkf chem_ic 1 - - typedef xchem_type real chem_ic_oin_a03 ijkf chem_ic 1 - - typedef xchem_type real chem_ic_oin_a04 ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_so4aj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_so4ai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_nh4aj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_nh4ai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_no3aj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_no3ai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_naaj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_naai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_claj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_clai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_orgpaj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_orgpai ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_ecj ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_eci ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_p25j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_p25i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_antha ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_seas ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_soila ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_nu0 ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_ac0 ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_corn ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa1j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa1i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa2j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa2i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa3j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa3i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa4j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_asoa4i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa1j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa1i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa2j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa2i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa3j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa3i ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa4j ijkf chem_ic 1 - - +typedef xchem_type real chem_ic_bsoa4i ijkf chem_ic 1 - - ##END XCHEM_TYPE DEFINITION @@ -201,9 +321,10 @@ state xchem_type vchem2 - - #CHEMISTRY PACKAGE DEFINITIONS package gocart_test chem_cv_options==10 - chem:sulf,p10,p25,so2,no2,o3,co,bc1,bc2,oc1,oc2,dust_1,dust_2,dust_3,dust_4,seas_1,seas_2,seas_3,seas_4;chem_ic:chem_ic_sulf,chem_ic_p10,chem_ic_p25,chem_ic_so2,chem_ic_no2,chem_ic_o3,chem_ic_co,chem_ic_bc1,chem_ic_bc2,chem_ic_oc1,chem_ic_oc2,chem_ic_dust_1,chem_ic_dust_2,chem_ic_dust_3,chem_ic_dust_4,chem_ic_seas_1,chem_ic_seas_2,chem_ic_seas_3,chem_ic_seas_4 package mosaic_test chem_cv_options==20 - chem:bc_a01,bc_a02,bc_a03,bc_a04,oc_a01,oc_a02,oc_a03,oc_a04,so4_a01,so4_a02,so4_a03,so4_a04,no3_a01,no3_a02,no3_a03,no3_a04,nh4_a01,nh4_a02,nh4_a03,nh4_a04,cl_a01,cl_a02,cl_a03,cl_a04,na_a01,na_a02,na_a03,na_a04,oin_a01,oin_a02,oin_a03,oin_a04,so2,no2,o3,co;chem_ic:chem_ic_bc_a01,chem_ic_bc_a02,chem_ic_bc_a03,chem_ic_bc_a04,chem_ic_oc_a01,chem_ic_oc_a02,chem_ic_oc_a03,chem_ic_oc_a04,chem_ic_so4_a01,chem_ic_so4_a02,chem_ic_so4_a03,chem_ic_so4_a04,chem_ic_no3_a01,chem_ic_no3_a02,chem_ic_no3_a03,chem_ic_no3_a04,chem_ic_nh4_a01,chem_ic_nh4_a02,chem_ic_nh4_a03,chem_ic_nh4_a04,chem_ic_cl_a01,chem_ic_cl_a02,chem_ic_cl_a03,chem_ic_cl_a04,chem_ic_na_a01,chem_ic_na_a02,chem_ic_na_a03,chem_ic_na_a04,chem_ic_oin_a01,chem_ic_oin_a02,chem_ic_oin_a03,chem_ic_oin_a04,chem_ic_so2,chem_ic_no2,chem_ic_o3,chem_ic_co +package racm_soa_vbs_da chem_cv_options==108 - chem:so2,no2,o3,co,so4aj,so4ai,nh4aj,nh4ai,no3aj,no3ai,naaj,naai,claj,clai,asoa1j,asoa1i,asoa2j,asoa2i,asoa3j,asoa3i,asoa4j,asoa4i,bsoa1j,bsoa1i,bsoa2j,bsoa2i,bsoa3j,bsoa3i,bsoa4j,bsoa4i,orgpaj,orgpai,ecj,eci,p25j,p25i,antha,seas,soila;chem_ic:chem_ic_so2,chem_ic_no2,chem_ic_o3,chem_ic_co,chem_ic_so4aj,chem_ic_so4ai,chem_ic_nh4aj,chem_ic_nh4ai,chem_ic_no3aj,chem_ic_no3ai,chem_ic_naaj,chem_ic_naai,chem_ic_claj,chem_ic_clai,chem_ic_asoa1j,chem_ic_asoa1i,chem_ic_asoa2j,chem_ic_asoa2i,chem_ic_asoa3j,chem_ic_asoa3i,chem_ic_asoa4j,chem_ic_asoa4i,chem_ic_bsoa1j,chem_ic_bsoa1i,chem_ic_bsoa2j,chem_ic_bsoa2i,chem_ic_bsoa3j,chem_ic_bsoa3i,chem_ic_bsoa4j,chem_ic_bsoa4i,chem_ic_orgpaj,chem_ic_orgpai,chem_ic_ecj,chem_ic_eci,chem_ic_p25j,chem_ic_p25i,chem_ic_antha,chem_ic_seas,chem_ic_soila ### surface chem obs ### -state real - f chemic_surf - - - - "Chem IC Surface Obs." "" +state real - f chemic_surf - - - - "Chem IC Surface Obs." "" state real chemsi_pm25 f chemic_surf 1 - h{19} "chems_pm25" "pm25 observation" "ug/m^3" state real chemsi_pm10 f chemic_surf 1 - h{19} "chems_pm10" "pm10 observation" "ug/m^3" state real chemsi_so2 f chemic_surf 1 - h{19} "chems_so2" "so2 observation" "ug/m^3" diff --git a/var/da/da_chem_sfc/da_ao_stats_chem_sfc.inc b/var/da/da_chem_sfc/da_ao_stats_chem_sfc.inc index 40156f450b..d5976eac13 100644 --- a/var/da/da_chem_sfc/da_ao_stats_chem_sfc.inc +++ b/var/da/da_chem_sfc/da_ao_stats_chem_sfc.inc @@ -35,10 +35,6 @@ subroutine da_ao_stats_chem_sfc (stats_unit, iv, re, ob) do ichem=PARAM_FIRST_SCALAR, num_chemic_surf nchem=0 -! stats%maximum%chem(ichem) = maxmin_type (missing_r, 0, 0) -! stats%minimum%chem(ichem) = maxmin_type (missing_r, 0, 0) -! stats%average = residual_chem_sfc_type(0.0, 0.0, 0.0, 0.0, 0.0) -! stats%rms_err = stats%average do n=1, iv%info(chemic_surf)%nlocal if (iv%info(chemic_surf)%proc_domain(1,n)) then @@ -67,7 +63,7 @@ subroutine da_ao_stats_chem_sfc (stats_unit, iv, re, ob) write(unit=stats_unit, fmt='(/a/)') ' Diagnostics of AO for chem' if (chemicda_opt == 1 .or. chemicda_opt == 2) then call da_print_stats_chem_sfc(stats_unit, nchem, stats) - else if (chemicda_opt == 3 .or. chemicda_opt == 4 .or. chemicda_opt == 5) then + else if (chemicda_opt >= 3) then call da_print_stats_chem_sfc(stats_unit, nchem, stats, ichem) end if end if diff --git a/var/da/da_chem_sfc/da_calculate_grady_chem_sfc.inc b/var/da/da_chem_sfc/da_calculate_grady_chem_sfc.inc index d5ab3b7340..6d21681ba9 100644 --- a/var/da/da_chem_sfc/da_calculate_grady_chem_sfc.inc +++ b/var/da/da_chem_sfc/da_calculate_grady_chem_sfc.inc @@ -18,9 +18,8 @@ subroutine da_calculate_grady_chem_sfc(iv, re, jo_grad_y) do n=1, iv%info(chemic_surf)%nlocal if (iv%chemic_surf(n)%chem(ichem)%qc < obs_qc_pointer) re%chemic_surf(n)%chem(ichem) = 0.0 - if (iv%chemic_surf(n)%chem(ichem)%qc < obs_qc_pointer) re%chemic_surf(n)%chem(ichem) = 0.0 - - jo_grad_y%chemic_surf(n)%chem(ichem) = -re%chemic_surf(n)%chem(ichem) / (iv%chemic_surf(n)%chem(ichem)%error * iv%chemic_surf(n)%chem(ichem)%error) + jo_grad_y%chemic_surf(n)%chem(ichem) = -re%chemic_surf(n)%chem(ichem) / & + (iv%chemic_surf(n)%chem(ichem)%error * iv%chemic_surf(n)%chem(ichem)%error) end do end do diff --git a/var/da/da_chem_sfc/da_chem_sfc.f90 b/var/da/da_chem_sfc/da_chem_sfc.f90 index e60842aa77..dc47adc232 100644 --- a/var/da/da_chem_sfc/da_chem_sfc.f90 +++ b/var/da/da_chem_sfc/da_chem_sfc.f90 @@ -43,7 +43,16 @@ module da_chem_sfc p_chem_ic_so4_a01, p_chem_ic_so4_a02, p_chem_ic_so4_a03, p_chem_ic_so4_a04, p_chem_ic_no3_a01, p_chem_ic_no3_a02, p_chem_ic_no3_a03, p_chem_ic_no3_a04,& p_chem_ic_nh4_a01, p_chem_ic_nh4_a02, p_chem_ic_nh4_a03, p_chem_ic_nh4_a04, p_chem_ic_cl_a01, p_chem_ic_cl_a02, p_chem_ic_cl_a03, p_chem_ic_cl_a04, & p_chem_ic_na_a01, p_chem_ic_na_a02, p_chem_ic_na_a03, p_chem_ic_na_a04, p_chem_ic_oin_a01, p_chem_ic_oin_a02, p_chem_ic_oin_a03, p_chem_ic_oin_a04, & - p_chem_ic_so2, p_chem_ic_no2, p_chem_ic_o3, p_chem_ic_co + p_chem_ic_so2, p_chem_ic_no2, p_chem_ic_o3, p_chem_ic_co, & + p_chem_ic_so4aj, p_chem_ic_so4ai, p_chem_ic_nh4aj, p_chem_ic_nh4ai, & ! aerosols in racm_soa_vbs_da + p_chem_ic_no3aj, p_chem_ic_no3ai, p_chem_ic_naaj, p_chem_ic_naai, & + p_chem_ic_asoa1j, p_chem_ic_asoa1i, p_chem_ic_asoa2j, p_chem_ic_asoa2i, & + p_chem_ic_asoa3j, p_chem_ic_asoa3i, p_chem_ic_asoa4j, p_chem_ic_asoa4i, & + p_chem_ic_bsoa1j, p_chem_ic_bsoa1i, p_chem_ic_bsoa2j, p_chem_ic_bsoa2i, & + p_chem_ic_bsoa3j, p_chem_ic_bsoa3i, p_chem_ic_bsoa4j, p_chem_ic_bsoa4i, & + p_chem_ic_orgpaj, p_chem_ic_orgpai, p_chem_ic_ecj, p_chem_ic_eci, & + p_chem_ic_p25j, p_chem_ic_p25i, p_chem_ic_antha, p_chem_ic_seas, & + p_chem_ic_claj, p_chem_ic_clai, p_chem_ic_soila #endif use da_define_structures, only : iv_type, y_type, jo_type, maxmin_type, & diff --git a/var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc b/var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc index bb478504e9..fd06e55a14 100644 --- a/var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc +++ b/var/da/da_chem_sfc/da_get_innov_vector_chem_sfc.inc @@ -11,12 +11,18 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv) type(iv_type), intent(inout) :: iv ! O-B structure. integer, intent(inout) :: num_qcstat_conv(:,:,:,:) - integer :: n,ichem ! Loop counters + integer :: n,ichem,ipm ! Loop counters - real, allocatable :: model_chemic(:,:) real, allocatable :: model_chemic_surf(:,:) + real, allocatable :: model_chemic(:,:) real, allocatable :: model_rho(:,:) + real, parameter :: scaleL = 3000. ! scaling factor for obs error [meters] + + ! max thresholds of abs(o-b) in pm25, pm10, so2, no2, o3, co + real,dimension(6):: maxomb = (/100., 100., 0.20, 0.20, 0.20, 20./) + real :: err_o, err_r + if (trace_use) call da_trace_entry("da_get_innov_vector_chem_sfc") if ( it > 1 ) then @@ -29,60 +35,54 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv) ! [0.0] Get chem surf conc. - allocate (model_chemic(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) allocate (model_chemic_surf(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chemic_surf)) - allocate (model_rho(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_chemic (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_rho (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + ! [1.0] Get horizontal interpolation weights: do ichem = PARAM_FIRST_SCALAR ,num_chem - ! [1.0] Get horizontal interpolation weights: call da_interp_lin_2d (grid%xbchem%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) - call da_interp_lin_2d (grid%xb%rho(:,:,1), iv%info(chemic_surf), 1, model_rho(:,ichem)) + call da_interp_lin_2d (grid%xb %rho (:,:,1), iv%info(chemic_surf), 1, model_rho (:,ichem)) end do + ! [1.1] Compute prior observations at obs sites if (chem_cv_options == 10) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_p25)*(model_chemic(:,p_chem_ic_p25) + & - 1.375*(96.06/28.964*1000)*model_chemic(:,p_chem_ic_sulf) + & - model_chemic(:,p_chem_ic_bc1)+model_chemic(:,p_chem_ic_bc2)+1.8*(model_chemic(:,p_chem_ic_oc1)+model_chemic(:,p_chem_ic_oc2)) + & - model_chemic(:,p_chem_ic_dust_1)+0.286*model_chemic(:,p_chem_ic_dust_2)+model_chemic(:,p_chem_ic_seas_1) + & + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_p25)*(model_chemic(:,p_chem_ic_p25) + & + 1.375*(96.06/28.964*1000)*model_chemic(:,p_chem_ic_sulf) + & + model_chemic(:,p_chem_ic_bc1)+model_chemic(:,p_chem_ic_bc2)+1.8*(model_chemic(:,p_chem_ic_oc1)+model_chemic(:,p_chem_ic_oc2)) + & + model_chemic(:,p_chem_ic_dust_1)+0.286*model_chemic(:,p_chem_ic_dust_2)+model_chemic(:,p_chem_ic_seas_1) + & 0.942*model_chemic(:,p_chem_ic_seas_2)) -!!! for pm10 !!! -! model_chemic_surf(:,:,p_chemsi_pm10)=model_rho(:,:,p_chem_ic_p10)*(model_chemic(:,:,p_chem_ic_p25)+model_chemic(:,:,p_chem_ic_p10) + & -! 1.375*(96.06/28.964*1000)*model_chemic(:,:,p_chem_ic_sulf)+model_chemic(:,:,p_chem_ic_bc1)+model_chemic(:,:,p_chem_ic_bc2) + & -! 1.8*(model_chemic(:,:,p_chem_ic_oc1)+model_chemic(:,:,p_chem_ic_oc2))+model_chemic(:,:,p_chem_ic_dust_1)+model_chemic(:,:,p_chem_ic_dust_2) + & -! model_chemic(:,:,p_chem_ic_dust_3)+0.87*model_chemic(:,:,p_chem_ic_dust_4)+model_chemic(:,:,p_chem_ic_seas_1) + & -! model_chemic(:,:,p_chem_ic_seas_2)+model_chemic(:,:,p_chem_ic_seas_3)) - else if (chem_cv_options == 20) then if (chemicda_opt == 1 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01)+ & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02)+ & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) else if (chemicda_opt == 2 ) then - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02) + & - model_chemic(:,p_chem_ic_oc_a03)+model_chemic(:,p_chem_ic_oc_a04)+model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02) + & - model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a01)+model_chemic(:,p_chem_ic_no3_a02) + & - model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02) + & - model_chemic(:,p_chem_ic_cl_a03)+model_chemic(:,p_chem_ic_cl_a04)+model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02) + & - model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a01)+model_chemic(:,p_chem_ic_oin_a02) + & + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02) + & + model_chemic(:,p_chem_ic_oc_a03)+model_chemic(:,p_chem_ic_oc_a04)+model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02) + & + model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a01)+model_chemic(:,p_chem_ic_no3_a02)+ & + model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02)+ & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02) + & + model_chemic(:,p_chem_ic_cl_a03)+model_chemic(:,p_chem_ic_cl_a04)+model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02) + & + model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a01)+model_chemic(:,p_chem_ic_oin_a02) + & model_chemic(:,p_chem_ic_oin_a03)+model_chemic(:,p_chem_ic_oin_a04)) else if (chemicda_opt == 3 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & - model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & - model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01)+ & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02)+ & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & + model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) else if (chemicda_opt == 4 ) then model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) @@ -91,15 +91,15 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv) model_chemic_surf(:,p_chemsi_co)=model_rho(:,p_chem_ic_co)*(28.01/28.964*1000)*(model_chemic(:,p_chem_ic_co)) else if (chemicda_opt == 5 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01)+ & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02)+ & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & - model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & + model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) @@ -109,67 +109,100 @@ subroutine da_get_innov_vector_chem_sfc( it, num_qcstat_conv, grid, ob, iv) end if - end if + else if (chem_cv_options == 108) then ! racm_soa_vbs_da + + if (chemicda_opt == 1 .or. chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm2.5 + + model_chemic_surf(:,p_chemsi_pm25) = 0.0 + do ichem = p_chem_ic_so4aj, p_chem_ic_p25i + model_chemic_surf(:,p_chemsi_pm25) = model_chemic_surf(:,p_chemsi_pm25) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + + end if + + if (chemicda_opt == 2) then ! pm10 + + model_chemic_surf(:,p_chemsi_pm10) = 0.0 + do ichem = p_chem_ic_so4aj, p_chem_ic_soila + model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + + end if !(chemicda_opt == 2) then + + if (chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm10 after pm2.5 + + ! pm10 - pm2.5 residuals + model_chemic_surf(:,p_chemsi_pm10) = 0.0 + do ichem = p_chem_ic_antha, p_chem_ic_soila + model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + end if !(chemicda_opt == 3 .or. chemicda_opt == 5) + + if (chemicda_opt >= 4) then + + model_chemic_surf(:,p_chemsi_so2)=model_chemic(:,p_chem_ic_so2) + model_chemic_surf(:,p_chemsi_no2)=model_chemic(:,p_chem_ic_no2) + model_chemic_surf(:,p_chemsi_o3) =model_chemic(:,p_chem_ic_o3 ) + model_chemic_surf(:,p_chemsi_co) =model_chemic(:,p_chem_ic_co ) + + end if !(chemicda_opt >= 4 ) then + + end if ! (chem_cv_options == 108) + + ! Loop over chemic_surf (e.g., surfce concentration) do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + ipm = ichem - PARAM_FIRST_SCALAR + 1 do n=iv%info(chemic_surf)%n1,iv%info(chemic_surf)%n2 ! [2.0] Initialise components of innovation vector: + iv % chemic_surf(n) % chem(ichem) % inv = missing_r + iv % chemic_surf(n) % chem(ichem) % error = missing_r - iv % chemic_surf(n) % chem(ichem) % inv = 0.0 - - ! [3.0] To compute the 'inv': + ! [3.0] To compute the innovation (inv = o-f): if (ob % chemic_surf(n) % chem(ichem) > missing_r .and. & - iv % chemic_surf(n) % chem(ichem) % qc >= obs_qc_pointer) then - iv % chemic_surf(n) % chem(ichem) % inv = ob % chemic_surf(n) % chem(ichem) - model_chemic_surf(n,ichem) - -! iv % chemic_surf(n) % chem(ichem) % error = sqrt((1.0+0.0075*ob % chemic_surf(n) % chem(ichem))**2+ & -! 0.5*(1.0+0.0075*ob % chemic_surf(n) % chem(ichem))*sqrt(grid%dx/2000.0)) - - if (chemicda_opt == 3 .or. chemicda_opt == 5 ) then - if (ichem == PARAM_FIRST_SCALAR+1 ) then - iv % chemic_surf(n) % chem(ichem) % error = 1.0*sqrt((1.0+0.0075*ob % chemic_surf(n) % chem(ichem))**2+ & - 0.5*(1.0+0.0075*ob % chemic_surf(n) % chem(ichem))*sqrt(grid%dx/2000.0)) - else - iv % chemic_surf(n) % chem(ichem) % error = 1.5*sqrt((1.0+0.0075*ob % chemic_surf(n) % chem(ichem))**2+ & - 0.5*(1.0+0.0075*ob % chemic_surf(n) % chem(ichem))*sqrt(grid%dx/2000.0)) - end if - else - iv % chemic_surf(n) % chem(ichem) % error = 1.5*sqrt((1.0+0.0075*ob % chemic_surf(n) % chem(ichem))**2+ & - 0.5*(1.0+0.0075*ob % chemic_surf(n) % chem(ichem))*sqrt(grid%dx/2000.0)) - end if + iv % chemic_surf(n) % chem(ichem) % qc >= obs_qc_pointer) then + + iv % chemic_surf(n) % chem(ichem) % inv = & + ob % chemic_surf(n) % chem(ichem) - model_chemic_surf(n,ichem) + + if ( ichem <= PARAM_FIRST_SCALAR+1 ) then ! PM2.5 and PM10 + + err_o = 1.5 + 0.0075 * ob % chemic_surf(n) % chem(ichem) + err_r = 0.5 * err_o * sqrt(grid % dx / scaleL) + iv % chemic_surf(n) % chem(ichem) % error = sqrt(err_o**2 + err_r**2) - if (chemicda_opt == 4 .or. chemicda_opt == 5 ) then - if (ichem == num_chemic_surf ) then - if (abs(iv % chemic_surf(n) % chem(ichem) % inv) .ge. 120*1000) then - iv % chemic_surf(n) % chem(ichem) % inv = missing_r - iv % chemic_surf(n) % chem(ichem) % qc = missing - iv % chemic_surf(n) % chem(ichem) % error = missing_r - end if - else - if (abs(iv % chemic_surf(n) % chem(ichem) % inv) .ge. 120) then - iv % chemic_surf(n) % chem(ichem) % inv = missing_r - iv % chemic_surf(n) % chem(ichem) % qc = missing - iv % chemic_surf(n) % chem(ichem) % error = missing_r - end if - end if - else - if (abs(iv % chemic_surf(n) % chem(ichem) % inv) .ge. 120) then - iv % chemic_surf(n) % chem(ichem) % inv = missing_r - iv % chemic_surf(n) % chem(ichem) % qc = missing - iv % chemic_surf(n) % chem(ichem) % error = missing_r - end if + else ! gas species (so2, no2, o3, co) + + err_o = 0.10 * ob % chemic_surf(n) % chem(ichem) ! obs error as 10% of the observation value + iv % chemic_surf(n) % chem(ichem) % error = err_o + + !err_o = 1.0 + 0.0075 * ob % chemic_surf(n) % chem(ichem) ! Wei Sun + !err_r = 0.5 * err_o * sqrt(grid % dx / 2000.) + !iv % chemic_surf(n) % chem(ichem) % error = 1.5*sqrt(err_o**2 + err_r**2) + + end if !( ichem <= PARAM_FIRST_SCALAR+1 ) + + ! Quality Control + ! ---------------- + if (abs(iv % chemic_surf(n) % chem(ichem) % inv) .ge. maxomb(ipm)) then + iv % chemic_surf(n) % chem(ichem) % inv = missing_r + iv % chemic_surf(n) % chem(ichem) % qc = missing + iv % chemic_surf(n) % chem(ichem) % error = missing_r end if - else - iv % chemic_surf(n) % chem(ichem) % inv = missing_r - end if + + end if ! good obs + end do - end do + end do ! do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + deallocate(model_rho) deallocate(model_chemic) deallocate(model_chemic_surf) - deallocate(model_rho) ! ----------------------------------------------------------------------- ! [4.0] Perform optional maximum error check: diff --git a/var/da/da_chem_sfc/da_residual_chem_sfc.inc b/var/da/da_chem_sfc/da_residual_chem_sfc.inc index bb3dc88258..3f986dfe40 100644 --- a/var/da/da_chem_sfc/da_residual_chem_sfc.inc +++ b/var/da/da_chem_sfc/da_residual_chem_sfc.inc @@ -29,7 +29,8 @@ subroutine da_residual_chem_sfc(iv, y, re, np_missing, np_bad_data, np_obs_used, allocate ( re%chemic_surf(n)%chem(num_chemic_surf)) do ichem = PARAM_FIRST_SCALAR, num_chemic_surf np_available = np_available + 1 - re%chemic_surf(n)%chem(ichem) = da_residual(n, 0, y%chemic_surf(n)%chem(ichem), iv%chemic_surf(n)%chem(ichem), n_obs_bad % chemic_surf) + re%chemic_surf(n)%chem(ichem) = da_residual(n, 0, y%chemic_surf(n)%chem(ichem), & + iv%chemic_surf(n)%chem(ichem), n_obs_bad % chemic_surf) end do end do diff --git a/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc.inc b/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc.inc index 9259243e73..f751db53c8 100644 --- a/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc.inc +++ b/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc.inc @@ -14,109 +14,140 @@ subroutine da_transform_xtoy_chem_sfc (grid, iv, y) integer :: n,ichem ! Loop counter. - real, allocatable :: model_chemic(:,:) - real, allocatable :: model_chemic_surf(:,:) - real, allocatable :: model_rho(:,:) + real, allocatable :: model_rho(:,:) ! model rho for each species at obs sites [kg/m^3] + real, allocatable :: model_chemic(:,:) ! model increments in each aerosol/gas species at obs sites + real, allocatable :: model_chemic_surf(:,:) ! model increments in each observation variable at obs sites if (trace_use_dull) call da_trace_entry("da_transform_xtoy_chem_sfc") - if (sfc_assi_options == sfc_assi_options_1) then - allocate (model_chemic(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) - allocate (model_chemic_surf(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chemic_surf)) - allocate (model_rho(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_rho (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_chemic (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_chemic_surf(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chemic_surf)) - model_chemic = 0.0 - model_chemic_surf = 0.0 - model_rho = 0.0 + model_rho = 0.0 + model_chemic = 0.0 + model_chemic_surf = 0.0 + ! [1.0] Get horizontal interpolation weights: do ichem = PARAM_FIRST_SCALAR ,num_chem - ! [1.0] Get horizontal interpolation weights: - call da_interp_lin_2d (grid%xachem%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) - call da_interp_lin_2d (grid%xb%rho(:,:,1), iv%info(chemic_surf), 1, model_rho(:,ichem)) + call da_interp_lin_2d (grid%xachem%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) + call da_interp_lin_2d (grid%xb%rho(:,:,1), iv%info(chemic_surf), 1, model_rho (:,ichem)) end do + ! [1.1] Compute prior observations at obs sites if (chem_cv_options == 10) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_p25)*(model_chemic(:,p_chem_ic_p25) + & - 1.375*(96.06/28.964*1000)*model_chemic(:,p_chem_ic_sulf) + & - model_chemic(:,p_chem_ic_bc1)+model_chemic(:,p_chem_ic_bc2)+1.8*(model_chemic(:,p_chem_ic_oc1)+model_chemic(:,p_chem_ic_oc2)) + & - model_chemic(:,p_chem_ic_dust_1)+0.286*model_chemic(:,p_chem_ic_dust_2)+model_chemic(:,p_chem_ic_seas_1) + & - 0.942*model_chemic(:,p_chem_ic_seas_2)) - -! model_chemic_surf(:,:,p_chemsi_pm10)=model_rho(:,:,p_chem_ic_p10)*(model_chemic(:,:,p_chem_ic_p25)+model_chemic(:,:,p_chem_ic_p10) + & -! 1.375*(96.06/28.964*1000)*model_chemic(:,:,p_chem_ic_sulf)+model_chemic(:,:,p_chem_ic_bc1)+model_chemic(:,:,p_chem_ic_bc2) + & -! 1.8*(model_chemic(:,:,p_chem_ic_oc1)+model_chemic(:,:,p_chem_ic_oc2))+model_chemic(:,:,p_chem_ic_dust_1)+model_chemic(:,:,p_chem_ic_dust_2) + & -! model_chemic(:,:,p_chem_ic_dust_3)+0.87*model_chemic(:,:,p_chem_ic_dust_4)+model_chemic(:,:,p_chem_ic_seas_1) + & -! model_chemic(:,:,p_chem_ic_seas_2)+model_chemic(:,:,p_chem_ic_seas_3)) + + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_p25)*(model_chemic(:,p_chem_ic_p25) + & + 1.375*(96.06/28.964*1000)*model_chemic(:,p_chem_ic_sulf) + & + model_chemic(:,p_chem_ic_bc1)+model_chemic(:,p_chem_ic_bc2) + & + 1.8*(model_chemic(:,p_chem_ic_oc1)+model_chemic(:,p_chem_ic_oc2)) + & + model_chemic(:,p_chem_ic_dust_1) + 0.286*model_chemic(:,p_chem_ic_dust_2) + & + model_chemic(:,p_chem_ic_seas_1) + 0.942*model_chemic(:,p_chem_ic_seas_2)) else if (chem_cv_options == 20) then - if (chemicda_opt == 1 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & - model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) - else if (chemicda_opt == 2 ) then - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02) + & - model_chemic(:,p_chem_ic_oc_a03)+model_chemic(:,p_chem_ic_oc_a04)+model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02) + & - model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a01)+model_chemic(:,p_chem_ic_no3_a02) + & - model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02) + & - model_chemic(:,p_chem_ic_cl_a03)+model_chemic(:,p_chem_ic_cl_a04)+model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02) + & - model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a01)+model_chemic(:,p_chem_ic_oin_a02) + & - model_chemic(:,p_chem_ic_oin_a03)+model_chemic(:,p_chem_ic_oin_a04)) - else if (chemicda_opt == 3 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & - model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & - model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & - model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) - else if (chemicda_opt == 4 ) then - model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) - model_chemic_surf(:,p_chemsi_no2)=model_rho(:,p_chem_ic_no2)*(46.01/28.964*1000)*(model_chemic(:,p_chem_ic_no2)) - model_chemic_surf(:,p_chemsi_o3)=model_rho(:,p_chem_ic_o3)*(48.00/28.964*1000)*(model_chemic(:,p_chem_ic_o3)) - model_chemic_surf(:,p_chemsi_co)=model_rho(:,p_chem_ic_co)*(28.01/28.964*1000)*(model_chemic(:,p_chem_ic_co)) - - else if (chemicda_opt == 5 ) then - model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & - model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & - model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & - model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & - model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & - model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & - model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) - model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & - model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & - model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) - - model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) - model_chemic_surf(:,p_chemsi_no2)=model_rho(:,p_chem_ic_no2)*(46.01/28.964*1000)*(model_chemic(:,p_chem_ic_no2)) - model_chemic_surf(:,p_chemsi_o3)=model_rho(:,p_chem_ic_o3)*(48.00/28.964*1000)*(model_chemic(:,p_chem_ic_o3)) - model_chemic_surf(:,p_chemsi_co)=model_rho(:,p_chem_ic_co)*(28.01/28.964*1000)*(model_chemic(:,p_chem_ic_co)) - - end if - - end if - - do ichem = PARAM_FIRST_SCALAR, num_chemic_surf - ! [1.2] Interpolate horizontally: - do n=iv%info(chemic_surf)%n1,iv%info(chemic_surf)%n2 - y%chemic_surf(n)%chem(ichem) = model_chemic_surf(n,ichem) - end do + + if (chemicda_opt == 1 ) then + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) + else if (chemicda_opt == 2 ) then + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02) + & + model_chemic(:,p_chem_ic_oc_a03)+model_chemic(:,p_chem_ic_oc_a04)+model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02) + & + model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a01)+model_chemic(:,p_chem_ic_no3_a02) + & + model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02) + & + model_chemic(:,p_chem_ic_cl_a03)+model_chemic(:,p_chem_ic_cl_a04)+model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02) + & + model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a01)+model_chemic(:,p_chem_ic_oin_a02) + & + model_chemic(:,p_chem_ic_oin_a03)+model_chemic(:,p_chem_ic_oin_a04)) + else if (chemicda_opt == 3 ) then + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & + model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & + model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) + else if (chemicda_opt == 4 ) then + model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) + model_chemic_surf(:,p_chemsi_no2)=model_rho(:,p_chem_ic_no2)*(46.01/28.964*1000)*(model_chemic(:,p_chem_ic_no2)) + model_chemic_surf(:,p_chemsi_o3)=model_rho(:,p_chem_ic_o3)*(48.00/28.964*1000)*(model_chemic(:,p_chem_ic_o3)) + model_chemic_surf(:,p_chemsi_co)=model_rho(:,p_chem_ic_co)*(28.01/28.964*1000)*(model_chemic(:,p_chem_ic_co)) + else if (chemicda_opt == 5 ) then + model_chemic_surf(:,p_chemsi_pm25)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic(:,p_chem_ic_bc_a01)+model_chemic(:,p_chem_ic_bc_a02) + & + model_chemic(:,p_chem_ic_bc_a03)+model_chemic(:,p_chem_ic_oc_a01)+model_chemic(:,p_chem_ic_oc_a02)+model_chemic(:,p_chem_ic_oc_a03) + & + model_chemic(:,p_chem_ic_so4_a01)+model_chemic(:,p_chem_ic_so4_a02)+model_chemic(:,p_chem_ic_so4_a03)+model_chemic(:,p_chem_ic_no3_a01) + & + model_chemic(:,p_chem_ic_no3_a02)+model_chemic(:,p_chem_ic_no3_a03)+model_chemic(:,p_chem_ic_nh4_a01)+model_chemic(:,p_chem_ic_nh4_a02) + & + model_chemic(:,p_chem_ic_nh4_a03)+model_chemic(:,p_chem_ic_cl_a01)+model_chemic(:,p_chem_ic_cl_a02)+model_chemic(:,p_chem_ic_cl_a03) + & + model_chemic(:,p_chem_ic_na_a01)+model_chemic(:,p_chem_ic_na_a02)+model_chemic(:,p_chem_ic_na_a03)+model_chemic(:,p_chem_ic_oin_a01) + & + model_chemic(:,p_chem_ic_oin_a02)+model_chemic(:,p_chem_ic_oin_a03)) + model_chemic_surf(:,p_chemsi_pm10)=model_rho(:,p_chem_ic_bc_a04)*(model_chemic(:,p_chem_ic_bc_a04)+model_chemic(:,p_chem_ic_oc_a04) + & + model_chemic(:,p_chem_ic_so4_a04)+model_chemic(:,p_chem_ic_no3_a04)+model_chemic(:,p_chem_ic_nh4_a04)+model_chemic(:,p_chem_ic_cl_a04) + & + model_chemic(:,p_chem_ic_na_a04)+model_chemic(:,p_chem_ic_oin_a04)) + model_chemic_surf(:,p_chemsi_so2)=model_rho(:,p_chem_ic_so2)*(64.06/28.964*1000)*(model_chemic(:,p_chem_ic_so2)) + model_chemic_surf(:,p_chemsi_no2)=model_rho(:,p_chem_ic_no2)*(46.01/28.964*1000)*(model_chemic(:,p_chem_ic_no2)) + model_chemic_surf(:,p_chemsi_o3)=model_rho(:,p_chem_ic_o3)*(48.00/28.964*1000)*(model_chemic(:,p_chem_ic_o3)) + model_chemic_surf(:,p_chemsi_co)=model_rho(:,p_chem_ic_co)*(28.01/28.964*1000)*(model_chemic(:,p_chem_ic_co)) + end if + + else if (chem_cv_options == 108) then ! racm_soa_vbs_da + + if (chemicda_opt == 1 .or. chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm2.5 + model_chemic_surf(:,p_chemsi_pm25) = 0.0 + do ichem = p_chem_ic_so4aj, p_chem_ic_p25i + model_chemic_surf(:,p_chemsi_pm25) = model_chemic_surf(:,p_chemsi_pm25) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + end if + + if (chemicda_opt == 2) then ! pm10 only + model_chemic_surf(:,p_chemsi_pm10) = 0.0 + do ichem = p_chem_ic_so4aj, p_chem_ic_soila + model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + end if + + if (chemicda_opt == 3 .or. chemicda_opt == 5) then ! pm10 after pm2.5 + ! pm10 - pm2.5 residual + model_chemic_surf(:,p_chemsi_pm10) = 0.0 + do ichem = p_chem_ic_antha, p_chem_ic_soila + model_chemic_surf(:,p_chemsi_pm10) = model_chemic_surf(:,p_chemsi_pm10) + & + model_rho(:,ichem) * model_chemic(:,ichem) + end do + end if + + if ( chemicda_opt >= 4 ) then + model_chemic_surf(:,p_chemsi_so2)=model_chemic(:,p_chem_ic_so2) + model_chemic_surf(:,p_chemsi_no2)=model_chemic(:,p_chem_ic_no2) + model_chemic_surf(:,p_chemsi_o3) =model_chemic(:,p_chem_ic_o3 ) + model_chemic_surf(:,p_chemsi_co) =model_chemic(:,p_chem_ic_co ) + end if + + end if ! if (chem_cv_options == 108) then + + ! [1.2] Interpolate horizontally: + do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + do n=iv%info(chemic_surf)%n1,iv%info(chemic_surf)%n2 + if(iv%chemic_surf(n)%chem(ichem)%qc >= obs_qc_pointer) then + y%chemic_surf(n)%chem(ichem) = model_chemic_surf(n,ichem) + else + y%chemic_surf(n)%chem(ichem) = 0.0 + end if end do + end do - deallocate (model_chemic) - end if + deallocate (model_rho) + deallocate (model_chemic) + deallocate (model_chemic_surf) if (trace_use_dull) call da_trace_exit("da_transform_xtoy_chem_sfc") end subroutine da_transform_xtoy_chem_sfc - - diff --git a/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc_adj.inc b/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc_adj.inc index 83334a6118..0ade89a53e 100644 --- a/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc_adj.inc +++ b/var/da/da_chem_sfc/da_transform_xtoy_chem_sfc_adj.inc @@ -21,45 +21,45 @@ subroutine da_transform_xtoy_chem_sfc_adj(grid, iv, jo_grad_y, jo_grad_x) if (trace_use_dull) call da_trace_entry("da_transform_xtoy_chem_sfc_adj") - if (sfc_assi_options == sfc_assi_options_1) then - - allocate (model_chemic(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) - allocate (model_chemic_surf(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chemic_surf)) - allocate (model_rho(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) - - jo_grad_x%chem_ic = 0.0 - model_chemic = 0.0 - model_chemic_surf = 0.0 - model_rho = 0.0 - - do ichem = PARAM_FIRST_SCALAR, num_chemic_surf - ! [1.2] Interpolate horizontally: - do n=iv%info(chemic_surf)%n1,iv%info(chemic_surf)%n2 + allocate (model_rho (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_chemic (iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chem)) + allocate (model_chemic_surf(iv%info(chemic_surf)%n1:iv%info(chemic_surf)%n2,num_chemic_surf)) + + model_rho = 0.0 + model_chemic = 0.0 + model_chemic_surf = 0.0 + jo_grad_x%chem_ic = 0.0 + + do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + do n=iv%info(chemic_surf)%n1,iv%info(chemic_surf)%n2 + if(iv%chemic_surf(n)%chem(ichem)%qc >= obs_qc_pointer) then model_chemic_surf(n,ichem) = jo_grad_y%chemic_surf(n)%chem(ichem) - end do + end if end do - - do ichem = PARAM_FIRST_SCALAR ,num_chem - call da_interp_lin_2d (grid%xb%rho(:,:,1), iv%info(chemic_surf), 1, model_rho(:,ichem)) - end do - - if (chem_cv_options == 10) then - model_chemic(:,p_chem_ic_p25)=model_rho(:,p_chem_ic_p25)*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_sulf)=model_rho(:,p_chem_ic_sulf)*1.375*(model_chemic_surf(:,p_chemsi_pm25))*(96.06/28.964*1000) - model_chemic(:,p_chem_ic_bc1)=model_rho(:,p_chem_ic_bc1)*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_bc2)=model_rho(:,p_chem_ic_bc2)*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_oc1)=model_rho(:,p_chem_ic_oc1)*1.8*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_oc2)=model_rho(:,p_chem_ic_oc2)*1.8*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_dust_1)=model_rho(:,p_chem_ic_dust_1)*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_dust_2)=model_rho(:,p_chem_ic_dust_2)*0.286*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_seas_1)=model_rho(:,p_chem_ic_seas_1)*(model_chemic_surf(:,p_chemsi_pm25)) - model_chemic(:,p_chem_ic_seas_2)=model_rho(:,p_chem_ic_seas_2)*0.942*(model_chemic_surf(:,p_chemsi_pm25)) - - do ichem = PARAM_FIRST_SCALAR, num_chem + end do + + do ichem = PARAM_FIRST_SCALAR ,num_chem + call da_interp_lin_2d (grid%xb%rho(:,:,1), iv%info(chemic_surf), 1, model_rho(:,ichem)) + end do + + if (chem_cv_options == 10) then + model_chemic(:,p_chem_ic_p25)=model_rho(:,p_chem_ic_p25)*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_sulf)=model_rho(:,p_chem_ic_sulf)*1.375*(model_chemic_surf(:,p_chemsi_pm25))*(96.06/28.964*1000) + model_chemic(:,p_chem_ic_bc1)=model_rho(:,p_chem_ic_bc1)*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_bc2)=model_rho(:,p_chem_ic_bc2)*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_oc1)=model_rho(:,p_chem_ic_oc1)*1.8*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_oc2)=model_rho(:,p_chem_ic_oc2)*1.8*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_dust_1)=model_rho(:,p_chem_ic_dust_1)*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_dust_2)=model_rho(:,p_chem_ic_dust_2)*0.286*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_seas_1)=model_rho(:,p_chem_ic_seas_1)*(model_chemic_surf(:,p_chemsi_pm25)) + model_chemic(:,p_chem_ic_seas_2)=model_rho(:,p_chem_ic_seas_2)*0.942*(model_chemic_surf(:,p_chemsi_pm25)) + + do ichem = PARAM_FIRST_SCALAR, num_chem call da_interp_lin_2d_adj (jo_grad_x%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) - end do + end do + + else if (chem_cv_options == 20) then - else if (chem_cv_options == 20) then if (chemicda_opt == 1 ) then model_chemic(:,p_chem_ic_bc_a01)=model_rho(:,p_chem_ic_bc_a01)*(model_chemic_surf(:,p_chemsi_pm25)) model_chemic(:,p_chem_ic_bc_a02)=model_rho(:,p_chem_ic_bc_a02)*(model_chemic_surf(:,p_chemsi_pm25)) @@ -263,17 +263,46 @@ subroutine da_transform_xtoy_chem_sfc_adj(grid, iv, jo_grad_y, jo_grad_x) call da_interp_lin_2d_adj (jo_grad_x%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) end do + end if ! if (chemicda_opt == 1 ) then + + else if (chem_cv_options == 108) then ! racm_soa_vbs_da + + if ( chemicda_opt == 1 .or. chemicda_opt == 3 .or. chemicda_opt == 5 ) then ! pm2.5 + do ichem = p_chem_ic_so4aj, p_chem_ic_p25i + model_chemic(:,ichem) = model_chemic(:,ichem) + model_rho(:,ichem) * model_chemic_surf(:,p_chemsi_pm25) + end do + end if !(chemicda_opt == 1 .or. chemicda_opt == 3 .or. chemicda_opt == 5) + + if ( chemicda_opt == 2 ) then ! pm10 only + do ichem = p_chem_ic_so4aj, p_chem_ic_soila + model_chemic(:,ichem) = model_chemic(:,ichem) + model_rho(:,ichem) * model_chemic_surf(:,p_chemsi_pm10) + end do + end if !( chemicda_opt == 2 ) + + if ( chemicda_opt == 3 .or. chemicda_opt == 5 ) then ! pm10 after pm2.5 + ! pm10 - pm2.5 residual + do ichem = p_chem_ic_antha, p_chem_ic_soila + model_chemic(:,ichem) = model_chemic(:,ichem) + model_rho(:,ichem) * model_chemic_surf(:,p_chemsi_pm10) + end do + end if !( chemicda_opt == 3 .or. chemicda_opt == 5 ) + + if ( chemicda_opt >= 4 ) then + model_chemic(:,p_chem_ic_so2)=model_chemic_surf(:,p_chemsi_so2) + model_chemic(:,p_chem_ic_no2)=model_chemic_surf(:,p_chemsi_no2) + model_chemic(:,p_chem_ic_o3 )=model_chemic_surf(:,p_chemsi_o3 ) + model_chemic(:,p_chem_ic_co )=model_chemic_surf(:,p_chemsi_co ) end if - end if + do ichem = PARAM_FIRST_SCALAR, num_chem + call da_interp_lin_2d_adj (jo_grad_x%chem_ic(:,:,1,ichem), iv%info(chemic_surf), 1, model_chemic(:,ichem)) + end do - deallocate(model_chemic) - deallocate(model_chemic_surf) - deallocate(model_rho) - end if + end if ! if (chem_cv_options == 10) then + + deallocate(model_rho) + deallocate(model_chemic) + deallocate(model_chemic_surf) if (trace_use_dull) call da_trace_exit("da_transform_xtoy_chem_sfc_adj") end subroutine da_transform_xtoy_chem_sfc_adj - - diff --git a/var/da/da_define_structures/da_allocate_observations_chem_sfc.inc b/var/da/da_define_structures/da_allocate_observations_chem_sfc.inc index e6af9325b4..fb3d42cd23 100644 --- a/var/da/da_define_structures/da_allocate_observations_chem_sfc.inc +++ b/var/da/da_define_structures/da_allocate_observations_chem_sfc.inc @@ -12,43 +12,9 @@ subroutine da_allocate_observations_chem_sfc (iv) if (trace_use) call da_trace_entry("da_allocate_observations_chem_sfc") - if (iv%info(chemic_surf)%nlocal > 0) allocate(iv%chemic_surf (1:iv%info(chemic_surf)%nlocal)) - - i = chemic_surf - if (iv%info(i)%nlocal > 0) then - allocate (iv%info(i)%name(iv%info(i)%nlocal)) - allocate (iv%info(i)%platform(iv%info(i)%nlocal)) - allocate (iv%info(i)%id(iv%info(i)%nlocal)) - allocate (iv%info(i)%date_char(iv%info(i)%nlocal)) - allocate (iv%info(i)%levels(iv%info(i)%nlocal)) - allocate (iv%info(i)%lat(iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%lon(iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%elv(iv%info(i)%nlocal)) - allocate (iv%info(i)%pstar(iv%info(i)%nlocal)) - - allocate (iv%info(i)%slp(iv%info(i)%nlocal)) - allocate (iv%info(i)%pw(iv%info(i)%nlocal)) - - allocate (iv%info(i)%x (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%y (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%i (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%j (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%dx (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%dxm(kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%dy (kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%dym(kms:kme,iv%info(i)%nlocal)) - allocate (iv%info(i)%k (iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%dz (iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%dzm(iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%zk (iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%proc_domain(iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%thinned(iv%info(i)%max_lev,iv%info(i)%nlocal)) - allocate (iv%info(i)%obs_global_index(iv%info(i)%nlocal)) - - iv%info(i)%proc_domain(:,:) = .false. - iv%info(i)%thinned(:,:) = .false. - iv%info(i)%zk(:,:) = missing_r - end if + if (iv%info(chemic_surf)%nlocal > 0) allocate(iv%chemic_surf(1:iv%info(chemic_surf)%nlocal)) + + call da_allocate_obs_info(iv,chemic_surf) if (trace_use) call da_trace_exit("da_allocate_observations_chem_sfc") diff --git a/var/da/da_define_structures/da_deallocate_background_errors.inc b/var/da/da_define_structures/da_deallocate_background_errors.inc index 1a2b348f31..418b27cd5e 100644 --- a/var/da/da_define_structures/da_deallocate_background_errors.inc +++ b/var/da/da_define_structures/da_deallocate_background_errors.inc @@ -31,6 +31,8 @@ subroutine da_deallocate_background_errors (be) end if end do deallocate ( be % v12 ) + deallocate ( be % cv_mz_chemic ) + if ( num_chem >= PARAM_FIRST_SCALAR ) deallocate ( be % cv % sizechemic ) end if #endif ! Deallocate gridpoint errors: diff --git a/var/da/da_define_structures/da_define_structures.f90 b/var/da/da_define_structures/da_define_structures.f90 index 47afac264b..a703ff0aff 100644 --- a/var/da/da_define_structures/da_define_structures.f90 +++ b/var/da/da_define_structures/da_define_structures.f90 @@ -1111,7 +1111,7 @@ module da_define_structures type be_subtype integer :: mz ! Vertical truncation of errors. integer :: max_wave ! Global only - horizontal spectral truncation. - character*5 :: name ! Variable name. + character*10 :: name ! Variable name. real*8, pointer :: rf_alpha(:) ! RF scale length. real*8, pointer :: val(:,:) ! Local Standard dev./sqrt(eigenvalue). real*8, pointer :: evec(:,:,:) ! Local Vertical eigenvectors. diff --git a/var/da/da_main/da_update_firstguess.inc b/var/da/da_main/da_update_firstguess.inc index 4c8f8180cc..37802946b3 100755 --- a/var/da/da_main/da_update_firstguess.inc +++ b/var/da/da_main/da_update_firstguess.inc @@ -1198,9 +1198,7 @@ if ( use_chemic_surfobs ) then end_index1(2) = p%ed2 end_index1(3) = p%ed3 -!! do ic = 1, num_chem-1 do ic = PARAM_FIRST_SCALAR, num_chem -!!! do ic = 2,2 #ifdef DM_PARALLEL if (rootproc) then allocate( globbuf((ide-1-ids+3)*(jde-1-jds+3)*(kde-1-kds+3))) diff --git a/var/da/da_obs_io/da_final_write_obs_gas_sfc.inc b/var/da/da_obs_io/da_final_write_obs_gas_sfc.inc index 34045ee934..d1521e3aae 100644 --- a/var/da/da_obs_io/da_final_write_obs_gas_sfc.inc +++ b/var/da/da_obs_io/da_final_write_obs_gas_sfc.inc @@ -15,6 +15,10 @@ subroutine da_final_write_obs_gas_sfc(it,iv) character(len=filename_len), allocatable :: filename(:) character(len=filename_len) :: file + !local + character(len=4) :: typestr + typestr = 'chem' + if (chem_cv_options >= 108) typestr = 'gas' if (trace_use) call da_trace_entry("da_final_write_obs_gas_sfc") @@ -44,10 +48,10 @@ subroutine da_final_write_obs_gas_sfc(it,iv) end if call da_proc_sum_int(num_obs) if (num_obs > 0 .and. rootproc) then - write(omb_unit,'(a20,i8)')'chem', num_obs + write(omb_unit,'(a20,i8)') typestr, num_obs num_obs = 0 do k = 0,num_procs-1 - call da_read_omb_tmp(filename(k),iunit,num_obs,'chem',4,if_wind_sd) + call da_read_omb_tmp(filename(k),iunit,num_obs,typestr,3,if_wind_sd) end do end if diff --git a/var/da/da_obs_io/da_read_obs_chem_sfc.inc b/var/da/da_obs_io/da_read_obs_chem_sfc.inc index ef77d4b2ee..171acd70d9 100644 --- a/var/da/da_obs_io/da_read_obs_chem_sfc.inc +++ b/var/da/da_obs_io/da_read_obs_chem_sfc.inc @@ -4,7 +4,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) ! Purpose: Read surface chem IC observation files ! ! History: 03/2019 Creation (Wei Sun) - ! + ! 06/2021 Add maximum thresholds for each species (Soyoung Ha) !------------------------------------------------------------------------- implicit none @@ -13,12 +13,18 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) character(len=*), intent(in) :: filename type(domain), intent(in) :: grid ! first guess state. + ! Maximum thresholds for observation values - treat as missing above these values + ! FIXME: We may want to put these maximum thresholds in namelist later. + real, dimension(6) :: max_pm = (/500.,999.,2.0, 2.0, 2.0, 50./) + real :: pm25, pm10, so2, no2, o3, co + + character (len = 6) :: SiteID character (len = 10) :: fmt_name character (len = 160) :: info_string integer :: i, j, iost, nlevels, old_nlevels, fm,iunit - integer :: ichem + integer :: ichem, ipm real :: dattmp,lattmp,lontmp type (singl_level_type) :: platform @@ -33,7 +39,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) integer :: ndup, n, report, obs_index real*8 :: obs_time, analysis_time - integer :: iyear, imonth, iday, ihour, imin + integer :: iyear, imonth, iday, ihour, imin, isite real :: tdiff, dlat_earth,dlon_earth,crit integer :: itt,itx,iout real, allocatable :: in(:), out(:) @@ -43,6 +49,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) real :: v_p(kms:kme),v_h(kms:kme) if (trace_use) call da_trace_entry("da_read_obs_chem_sfc") + ! Initialize counts ntotal(:) = iv%info(:)%ptotal(iv%time-1) nlocal(:) = iv%info(:)%plocal(iv%time-1) @@ -54,6 +61,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) call cleangrids_conv(n) end do end if + ! open file ! ========= @@ -67,7 +75,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) if (iost /= 0) then write(unit=message(1),fmt='(A,I5,A)') & - "Error",iost," opening chemic obs file "//trim(filename) + "Error",iost," opening chem obs file "//trim(filename) call da_warning(__FILE__,__LINE__,message(1:1)) call da_free_unit(iunit) if (trace_use) call da_trace_exit("da_read_obs_chem_sfc") @@ -77,6 +85,10 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) allocate ( platform%chem (num_chemic_surf) ) obs_index = chemic_surf + write(*,'(A,6f7.1)') "da_read_obs_chem_sfc: pm25, pm10, so2, no2, o3, co obs are rejected if >= ", max_pm + + imin=0 ! minutes + ! loop over records ! ----------------- @@ -87,19 +99,16 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) report = report+1 - ! read station general info + ! Initialization ! ============================= + platform%chem(:)%inv = missing_r + platform%chem(:)%qc = missing + platform%chem(:)%error = missing_r -!!! platform%each (1) = chemic_surf_type(missing_r, missing, -1.0, & ! height -!!! field_type(missing_r, missing, missing_r, missing, missing_r)) - - platform%chem(:)%inv = missing_r - platform%chem(:)%qc = missing - platform%chem(:)%error = missing_r - - if (chem_cv_options == 10) then -!!!! ium data !!!! + ! read station general info + ! ============================= + if (chem_cv_options == 10) then !!!! ium data !!!! read (iunit, *, iostat = iost) & iyear, imonth, iday, ihour, imin, & platform%info%lat, & @@ -108,27 +117,11 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour platform%info%id = "333" -!!!! korus data !!!! -! read (iunit, *, iostat = iost) & -! platform%info%id, & -! platform%info%lat, & -! platform%info%lon, & -! platform%info%date_char, & -! platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! pm10 !!! -! platform%chem(PARAM_FIRST_SCALAR)%inv !!! pm25 !!! -!!! platform%chem(PARAM_FIRST_SCALAR+1)%inv -! read (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour -!!!!!!!!!!!!!!!!!!! - - if (platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR)%inv>=500) then -!!! if (platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR)%inv>=100) then + if (platform%chem(PARAM_FIRST_SCALAR)%inv==-999.or.platform%chem(PARAM_FIRST_SCALAR)%inv>=500) then platform%chem(PARAM_FIRST_SCALAR)%inv=missing_r else platform%chem(PARAM_FIRST_SCALAR)%qc=333 end if -!!! if (platform%chem(PARAM_FIRST_SCALAR+1)%inv==-999) platform%chem(PARAM_FIRST_SCALAR)%inv=missing_r - - imin=0 else if (chem_cv_options == 20) then if (chemicda_opt == 1) then @@ -142,7 +135,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) iyear, imonth, iday, ihour, imin, & platform%info%lat, & platform%info%lon, & - platform%chem(PARAM_FIRST_SCALAR)%inv, & + platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! pm25 read !!! platform%chem(PARAM_FIRST_SCALAR)%inv !!! pm10 only !!! else if (chemicda_opt == 3) then read (iunit, *, iostat = iost) & @@ -156,8 +149,8 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) iyear, imonth, iday, ihour, imin, & platform%info%lat, & platform%info%lon, & - platform%chem(PARAM_FIRST_SCALAR)%inv, & - platform%chem(PARAM_FIRST_SCALAR)%inv, & + platform%chem(PARAM_FIRST_SCALAR)%inv, & + platform%chem(PARAM_FIRST_SCALAR)%inv, & platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! so2 !!! platform%chem(PARAM_FIRST_SCALAR+1)%inv, & !!! no2 !!! platform%chem(PARAM_FIRST_SCALAR+2)%inv, & !!! o3 !!! @@ -239,7 +232,68 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) end do end if - end if + + else if (chem_cv_options == 108) then + + read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost) & + SiteID, & + platform%info%lat, & + platform%info%lon, & + iyear, imonth, iday, ihour, & + pm25, pm10, so2, no2, o3, co + + platform%info%id = trim(SiteID) + read(SiteID,'(I6)') isite + write (platform%info%date_char,'(i4,3(i2.2))') iyear, imonth, iday, ihour + + !!!!!!!!!!!!!!!!!!! + ! Caution: Observations are assumed to have the same unit of ppmv. + ! Unit conversion (to ppmv) must have been done in the preprocessing. + ! (ex. The original chinese data had ug/m3 for gas species.) + !!!!!!!!!!!!!!!!!!! + + select case ( chemicda_opt ) + case ( 1 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%info%name = "pm25" + case ( 2 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm10 + platform%info%name = "pm10" + case ( 3 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10 + platform%info%name = "all_pm" + case ( 4 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = so2 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = no2 + platform%chem(PARAM_FIRST_SCALAR+2)%inv = o3 + platform%chem(PARAM_FIRST_SCALAR+3)%inv = co + platform%info%name = "gas" + case ( 5 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10 + platform%chem(PARAM_FIRST_SCALAR+2)%inv = so2 + platform%chem(PARAM_FIRST_SCALAR+3)%inv = no2 + platform%chem(PARAM_FIRST_SCALAR+4)%inv = o3 + platform%chem(PARAM_FIRST_SCALAR+5)%inv = co + platform%info%name = "all" + case default + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%info%name = "" + end select + + ! Sanity check + !!!!!!!!!!!!!!!!!!! + do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + ipm = ichem - PARAM_FIRST_SCALAR + 1 + if (platform%chem(ichem)%inv<0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then + platform%chem(ichem)%inv=missing_r + else + platform%chem(ichem)%qc = 0 + endif + end do + + end if ! (chem_cv_options == 10) then if (iost /= 0) then ! FIX? This is expected, but its unclear how we handle failure @@ -257,18 +311,14 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) ! read model location ! ========================= -!!! read (iunit, fmt = fmt_srfc) & -!!! platform%loc%slp%inv, platform%loc%slp%qc, platform%loc%slp%error, & -!!! platform%loc%pw%inv, platform%loc%pw%qc, platform%loc%pw%error - ! Check if outside of the time range: call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time) if ( obs_time < time_slots(0) .or. & obs_time >= time_slots(num_fgat_time) ) then if (print_detail_obs) then - write(unit=stdout, fmt='(a)') '*** Outside of the time range:' - write(unit=stdout, fmt=fmt_info) & + write(unit=stdout, fmt='(a)') '*** Outside of the time range:' + write(unit=stdout, fmt=fmt_info) & platform%info%platform, & platform%info%date_char, & platform%info%name, & @@ -300,9 +350,12 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) if (print_detail_obs) then ! Simplistic approach, could be improved to get it all done on PE 0 if (.NOT. outside) then - write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') & - "Report",report," at",platform%info%lon,"E",platform%info%lat, & - "N on processor", myproc,", position", platform%loc%x,platform%loc%y + write(unit=stdout,fmt='(A,I5,A,A,3X,2F8.2,2I4,2F7.1,4F9.3)') & + "Report",report," at ",trim(platform%info%id),platform%info%lat,platform%info%lon, & + platform%loc%j, platform%loc%i, pm25, pm10, so2, no2, o3, co + !write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') & + ! "Report",report," at",platform%info%lon,"E",platform%info%lat, & + ! "N on processor", myproc,", position", platform%loc%x,platform%loc%y end if end if @@ -337,7 +390,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) platform%info%elv, & platform%info%id - write(unit=stdout, fmt = '(a,2i5,4e20.10)') & + write(unit=stdout, fmt = '(a,2i5,4f10.3)') & ' duplicating obs since loc% i,j,dx,dxm,dy & dym ', & platform%loc%i, platform%loc%j, & platform%loc%dx, platform%loc%dxm, & @@ -348,7 +401,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) dup_loop: do n = 1, ndup is_surface=.true. if (use_chemic_surfobs) then - if (.not. use_chemic_surfobs .or. ntotal(chemic_surf) == max_synop_input ) cycle reports + if ( ntotal(chemic_surf) == max_synop_input ) cycle reports if (n==1) ntotal(chemic_surf) = ntotal(chemic_surf)+1 if (outside) cycle reports if ( thin_conv_ascii ) then @@ -360,11 +413,12 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) ilocal(chemic_surf) = ilocal(chemic_surf) + 1 end if - allocate ( iv%chemic_surf(ilocal(chemic_surf))%chem(num_chemic_surf)) + allocate ( iv%chemic_surf(ilocal(chemic_surf))%chem(num_chemic_surf)) do ichem = PARAM_FIRST_SCALAR, num_chemic_surf iv%chemic_surf(ilocal(chemic_surf))%chem(ichem) = platform%chem(ichem) end do + else write(unit=message(1), fmt='(a)') 'unsaved obs found:' @@ -426,6 +480,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) ! thinning check if ( thin_conv_ascii ) then + print*, 'da_read_obs_chem_sfc: thin_conv_ascii ',thin_conv_ascii do n = 1, num_ob_indexes if (n==radar) cycle allocate ( in(thinning_grid_conv(n)%itxmax) ) @@ -466,6 +521,7 @@ subroutine da_read_obs_chem_sfc (iv, filename, grid) deallocate( out ) end do ! loop over num_ob_indexes end if ! thin_conv_ascii + if (trace_use) call da_trace_exit("da_read_obs_chem_sfc") end subroutine da_read_obs_chem_sfc diff --git a/var/da/da_obs_io/da_read_omb_tmp.inc b/var/da/da_obs_io/da_read_omb_tmp.inc index 915e7fc48a..e8be5874ef 100644 --- a/var/da/da_obs_io/da_read_omb_tmp.inc +++ b/var/da/da_obs_io/da_read_omb_tmp.inc @@ -18,7 +18,16 @@ subroutine da_read_omb_tmp(filename,unit_in,num,obs_type_in,nc,if_wind_sd) character*20 :: iv_type logical :: if_write, if_wind_sd - character*5 :: stn_id +! data format from da_write_obs_chem_sfc.inc + character(len=120) :: fmt_chem = '(i8,3x,a6,2f9.2,5(2f17.7,i8,2f17.7))' + character(len=120) :: fmt_chem1 = '(i8,2x,a6,2f11.6,2f11.2,i8,2f11.2)' + character(len=120) :: fmt_chem2 = '(i8,2x,a6,2f11.6,2(2f11.2,i8,2f11.2))' + character(len=120) :: fmt_chem4 = '(i8,2x,a6,2f11.6,4(2f12.3,i8,2f12.3))' + character(len=120) :: wrt_chem1 = '(a6,2f11.6,2f11.2,i8,2f11.2)' + character(len=120) :: wrt_chem2 = '(a6,2f11.6,2(2f11.2,i8,2f11.2))' + character(len=120) :: wrt_chem4 = '(a6,2f11.6,4(2f12.3,i8,2f12.3))' + + character*6 :: stn_id integer :: n, k, kk, l, levels, dummy_i real :: lat, lon, press, height, dummy real :: tpw_obs, tpw_inv, tpw_err, tpw_inc @@ -441,6 +450,7 @@ subroutine da_read_omb_tmp(filename,unit_in,num,obs_type_in,nc,if_wind_sd) #if (WRF_CHEM == 1) case ('chem' ) if (num_obs > 0) then + if (chem_cv_options > 1 .and. chem_cv_options < 100) then do n = 1, num_obs read(unit_in,'(2i8)')levels, ifgat if (if_write) then @@ -470,7 +480,7 @@ subroutine da_read_omb_tmp(filename,unit_in,num,obs_type_in,nc,if_wind_sd) lat, lon, dummy, & ! Lat/lon, dummy chem_obs, chem_inv, chem_qc, chem_err, chem_inc, & chem_obs2, chem_inv2, chem_qc2, chem_err2, chem_inc2 - else if (chemicda_opt == 4) then + else if (chemicda_opt == 4) then read(unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)& kk,l, stn_id, & ! Station lat, lon, dummy, & ! Lat/lon, dummy @@ -510,7 +520,64 @@ subroutine da_read_omb_tmp(filename,unit_in,num,obs_type_in,nc,if_wind_sd) end if end do end do - end if + + else if (chem_cv_options >= 108) then + + do n = 1, num_obs + if (chemicda_opt <= 2) then + read(unit_in, fmt = fmt_chem1, err= 1000)& + kk, stn_id, & ! Station + lat, lon, & ! Lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc + if (if_write) & + write(omb_unit, fmt = wrt_chem1, err= 1000)& + stn_id, lat, lon, & ! station, lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc + else if (chemicda_opt == 3 .or. chemicda_opt == 5) then + read(unit_in, fmt = fmt_chem2, err= 1000)& + kk, stn_id, & ! Station + lat, lon, & ! Lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc, & + chem_obs2, chem_inv2, chem_qc2, chem_err2, chem_inc2 + if (if_write) & + write(omb_unit, fmt = wrt_chem2, err= 1000)& + stn_id, lat, lon, & ! station, lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc, & + chem_obs2, chem_inv2, chem_qc2, chem_err2, chem_inc2 + end if + + end do ! n = 1, num_obs + + endif !(chem_cv_options >= 100) then + end if ! if (num_obs > 0) then + + if (if_write) exit reports + cycle reports + + case ('gas' ) + if (num_obs > 0) then + if (chem_cv_options >= 108) then + do n = 1, num_obs + if (chemicda_opt == 4 .or. chemicda_opt == 5) then + read(unit_in, fmt = fmt_chem4, err= 1000)& + kk, stn_id, & ! Station + lat, lon, & ! Lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc, & + chem_obs2, chem_inv2, chem_qc2, chem_err2, chem_inc2, & + chem_obs3, chem_inv3, chem_qc3, chem_err3, chem_inc3, & + chem_obs4, chem_inv4, chem_qc4, chem_err4, chem_inc4 + if (if_write) & + write(omb_unit, fmt = wrt_chem4, err= 1000)& + stn_id, lat, lon, & ! station, lat/lon + chem_obs, chem_inv, chem_qc, chem_err, chem_inc, & + chem_obs2, chem_inv2, chem_qc2, chem_err2, chem_inc2, & + chem_obs3, chem_inv3, chem_qc3, chem_err3, chem_inc3, & + chem_obs4, chem_inv4, chem_qc4, chem_err4, chem_inc4 + end if + end do ! do n = 1, num_obs + endif ! if (chem_cv_options >= 108) then + end if ! if (num_obs > 0) then + if (if_write) exit reports cycle reports #endif diff --git a/var/da/da_obs_io/da_scan_obs_chem_sfc.inc b/var/da/da_obs_io/da_scan_obs_chem_sfc.inc index 714cd0f961..a84f36d032 100644 --- a/var/da/da_obs_io/da_scan_obs_chem_sfc.inc +++ b/var/da/da_obs_io/da_scan_obs_chem_sfc.inc @@ -5,6 +5,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) ! ! History: 03/2019 Creation (Wei Sun) ! + ! 06/2021 Add maximum thresholds for each species (Soyoung Ha) !------------------------------------------------------------------------- implicit none @@ -13,12 +14,18 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) character(len=*), intent(in) :: filename type(domain), intent(in) :: grid ! first guess state. + ! Maximum thresholds for observation values - treat as missing above these values + ! FIXME: We may want to put these maximum thresholds in namelist later. + real, dimension(6) :: max_pm = (/500.,999.,2.0, 2.0, 2.0, 50./) + real :: pm25, pm10, so2, no2, o3, co + + character (len = 6) :: SiteID character (len = 10) :: fmt_name character (len = 160) :: info_string integer :: i, j, iost, nlevels, old_nlevels, fm,iunit - integer :: ichem + integer :: ichem, ipm real :: dattmp,lattmp,lontmp type (singl_level_type) :: platform @@ -33,7 +40,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) integer :: ndup, n, report, obs_index real*8 :: obs_time, analysis_time - integer :: iyear, imonth, iday, ihour, imin + integer :: iyear, imonth, iday, ihour, imin, isite real :: tdiff, dlat_earth,dlon_earth,crit integer :: itt,itx,iout real, allocatable :: in(:), out(:) @@ -43,6 +50,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) real :: v_p(kms:kme),v_h(kms:kme) if (trace_use) call da_trace_entry("da_scan_obs_chem_sfc") + ! Initialize counts if ( thin_conv_ascii ) then @@ -51,6 +59,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) call cleangrids_conv(n) end do end if + ! open file ! ========= @@ -64,7 +73,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) if (iost /= 0) then write(unit=message(1),fmt='(A,I5,A)') & - "Error",iost," opening chemic obs file "//trim(filename) + "Error",iost," opening chemic surface obs file "//trim(filename) call da_warning(__FILE__,__LINE__,message(1:1)) call da_free_unit(iunit) if (trace_use) call da_trace_exit("da_scan_obs_chem_sfc") @@ -74,6 +83,8 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) allocate ( platform%chem (num_chemic_surf) ) obs_index = chemic_surf + imin=0 ! minutes + ! loop over records ! ----------------- @@ -84,19 +95,16 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) report = report+1 - ! read station general info + ! Initialization ! ============================= + platform%chem(:)%inv = missing_r + platform%chem(:)%qc = missing + platform%chem(:)%error = missing_r -!!! platform%each (1) = chemic_surf_type(missing_r, missing, -1.0, & ! height -!!! field_type(missing_r, missing, missing_r, missing, missing_r)) - - platform%chem(:)%inv = missing_r - platform%chem(:)%qc = missing - platform%chem(:)%error = missing_r - - if (chem_cv_options == 10) then + ! read station general info + ! ============================= -!!!! ium data !!!! + if (chem_cv_options == 10) then !!!! ium data !!!! write(*,*) 'read ium data' read (iunit, *, iostat = iost) & @@ -107,20 +115,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour platform%info%id = "333" -!!!! korus data !!!! -! read (iunit, *, iostat = iost) & -! platform%info%id, & -! platform%info%lat, & -! platform%info%lon, & -! platform%info%date_char, & -! platform%chem(PARAM_FIRST_SCALAR)%inv -!!! platform%chem(PARAM_FIRST_SCALAR+1)%inv -! read (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour -!!!!!!!!!!!!!!!!!!! - - imin=0 - - else if (chem_cv_options == 20) then + else if (chem_cv_options == 20) then if (chemicda_opt == 1) then read (iunit, *, iostat = iost) & @@ -133,15 +128,15 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) iyear, imonth, iday, ihour, imin, & platform%info%lat, & platform%info%lon, & - platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! pm25 read !!! - platform%chem(PARAM_FIRST_SCALAR)%inv !!! pm10 only !!! + platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! pm25 read !!! + platform%chem(PARAM_FIRST_SCALAR)%inv !!! pm10 only !!! else if (chemicda_opt == 3) then read (iunit, *, iostat = iost) & iyear, imonth, iday, ihour, imin, & platform%info%lat, & platform%info%lon, & - platform%chem(PARAM_FIRST_SCALAR)%inv, & - platform%chem(PARAM_FIRST_SCALAR+1)%inv !!! pm10 & pm25 !!! + platform%chem(PARAM_FIRST_SCALAR)%inv, & + platform%chem(PARAM_FIRST_SCALAR+1)%inv !!! pm10 & pm25 !!! else if (chemicda_opt == 4) then read (iunit, *, iostat = iost) & iyear, imonth, iday, ihour, imin, & @@ -151,11 +146,8 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) platform%chem(PARAM_FIRST_SCALAR)%inv, & platform%chem(PARAM_FIRST_SCALAR)%inv, & !!! so2 !!! platform%chem(PARAM_FIRST_SCALAR+1)%inv, & !!! no2 !!! - platform%chem(PARAM_FIRST_SCALAR+2)%inv, & !!! o3 !!! - platform%chem(PARAM_FIRST_SCALAR+3)%inv !!! co !!! - - write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour - + platform%chem(PARAM_FIRST_SCALAR+2)%inv !!! o3 !!! + write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour else if (chemicda_opt == 5) then read (iunit, *, iostat = iost) & iyear, imonth, iday, ihour, imin, & @@ -167,13 +159,66 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) platform%chem(PARAM_FIRST_SCALAR+3)%inv, & !!! no2 !!! platform%chem(PARAM_FIRST_SCALAR+4)%inv, & !!! o3 !!! platform%chem(PARAM_FIRST_SCALAR+5)%inv !!! co !!! - - write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour - - end if - platform%info%id = "333" + write (platform%info%date_char,'(i4,3(i2))') iyear, imonth, iday, ihour + platform%info%id = "333" end if + else if (chem_cv_options == 108) then + + read (iunit, fmt='(A6,F10.6,1X,F10.6,1X,I4,3I2,2F9.1,4F9.3)', iostat=iost) & + SiteID, & + platform%info%lat, & + platform%info%lon, & + iyear, imonth, iday, ihour, & + pm25, pm10, so2, no2, o3, co + + platform%info%id = trim(SiteID) + read(SiteID,'(I6)') isite + write(platform%info%date_char,'(i4,3(i2.2))') iyear, imonth, iday, ihour + + select case ( chemicda_opt ) + case ( 1 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%info%name = "pm25" + case ( 2 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm10 + platform%info%name = "pm10" + case ( 3 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10 + platform%info%name = "all_pm" + case ( 4 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = so2 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = no2 + platform%chem(PARAM_FIRST_SCALAR+2)%inv = o3 + platform%chem(PARAM_FIRST_SCALAR+3)%inv = co + platform%info%name = "gas" + case ( 5 ) + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%chem(PARAM_FIRST_SCALAR+1)%inv = pm10 + platform%chem(PARAM_FIRST_SCALAR+2)%inv = so2 + platform%chem(PARAM_FIRST_SCALAR+3)%inv = no2 + platform%chem(PARAM_FIRST_SCALAR+4)%inv = o3 + platform%chem(PARAM_FIRST_SCALAR+5)%inv = co + platform%info%name = "all" + case default + platform%chem(PARAM_FIRST_SCALAR )%inv = pm25 + platform%info%name = "" + end select + + ! Sanity check + !!!!!!!!!!!!!!!!!!! + do ichem = PARAM_FIRST_SCALAR, num_chemic_surf + ipm = ichem - PARAM_FIRST_SCALAR + 1 + if (platform%chem(ichem)%inv<=0.or.abs(platform%chem(ichem)%inv)>=max_pm(ipm)) then + platform%chem(ichem)%inv=missing_r + else + platform%chem(ichem)%qc = 0 + endif + end do + + end if !(chem_cv_options == 108) + if (iost /= 0) then ! FIX? This is expected, but its unclear how we handle failure ! here without assuming the fortran2003 convention on @@ -190,10 +235,6 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) ! read model location ! ========================= -!!! read (iunit, fmt = fmt_srfc) & -!!! platform%loc%slp%inv, platform%loc%slp%qc, platform%loc%slp%error, & -!!! platform%loc%pw%inv, platform%loc%pw%qc, platform%loc%pw%error - ! Check if outside of the time range: call da_get_julian_time (iyear,imonth,iday,ihour,imin,obs_time) @@ -235,7 +276,7 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) if (.NOT. outside) then write(unit=stdout,fmt='(A,I5,A,F8.2,A,F8.2,A,I3,A,2F8.2)') & "Report",report," at",platform%info%lon,"E",platform%info%lat, & - "N on processor", myproc,", position", platform%loc%x,platform%loc%y + "N on proc", myproc," at x/y:", platform%loc%x,platform%loc%y end if end if @@ -305,7 +346,6 @@ subroutine da_scan_obs_chem_sfc (iv, filename, grid) end if - end do dup_loop end do reports diff --git a/var/da/da_obs_io/da_write_obs_chem_sfc.inc b/var/da/da_obs_io/da_write_obs_chem_sfc.inc index 0daf7171e2..2667c431c8 100644 --- a/var/da/da_obs_io/da_write_obs_chem_sfc.inc +++ b/var/da/da_obs_io/da_write_obs_chem_sfc.inc @@ -14,7 +14,15 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) integer :: n, k, num_obs, ios, ios2 integer :: ounit, ounit2 ! Output unit character(len=filename_len) :: filename, filename2 + + character(len=120) :: fmt_chem1 = '(i8,2x,a6,2f11.6,2f11.2,i8,2f11.2)' ! pm2.5 or pm10 + character(len=120) :: fmt_chem2 = '(i8,2x,a6,2f11.6,2(2f11.2,i8,2f11.2))' ! pm2.5 and pm10 + character(len=120) :: fmt_chem4 = '(i8,2x,a6,2f11.6,4(2f12.3,i8,2f12.3))' ! 4 gas species integer :: itime, ifgat + integer :: ipcc, jpcc, ic + integer, parameter :: ngas = 4 + integer, dimension(ngas) :: ipc + ipc = (/p_chemsi_so2, p_chemsi_no2, p_chemsi_o3, p_chemsi_co/) if (trace_use) call da_trace_entry("da_write_obs_chem_sfc") @@ -53,7 +61,7 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) iostat=ios) if (ios /= 0) then call da_error(__FILE__,__LINE__, & - (/"Cannot open conventional observation omb and oma file"//filename/)) + (/"Cannot open chemical observation omb and oma file"//filename/)) end if if (chemicda_opt == 5) then @@ -61,7 +69,7 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) iostat=ios2) if (ios2 /= 0) then call da_error(__FILE__,__LINE__, & - (/"Cannot open conventional observation omb and oma file"//filename2/)) + (/"Cannot open gas observation omb and oma file"//filename2/)) end if end if @@ -70,9 +78,13 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) if (iv%info(chemic_surf)%proc_domain(1,n)) num_obs = num_obs + 1 end do if (num_obs > 0) then - write(ounit,'(a20,i8)')'chem', num_obs + if (chemicda_opt == 4) then + write(ounit,'(a20,i8)') 'gas', num_obs + else + write(ounit,'(a20,i8)') 'chem', num_obs + endif if (chemicda_opt == 5) then - write(ounit2,'(a20,i8)')'chem', num_obs + write(ounit2,'(a20,i8)') 'gas', num_obs end if num_obs = 0 do n = 1, iv%info(chemic_surf)%nlocal @@ -86,95 +98,67 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) if (iv%info(chemic_surf)%proc_domain(1,n)) then num_obs = num_obs + 1 - write(ounit,'(2i8)') 1, ifgat - if (chemicda_opt == 5) then - write(ounit2,'(2i8)') 1, ifgat - end if - if (chemicda_opt == 1) then - write(ounit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_pm25), & - iv%chemic_surf(n)%chem(p_chemsi_pm25)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm25)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm25)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm25) - else if (chemicda_opt == 2) then - write(ounit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_pm10), & - iv%chemic_surf(n)%chem(p_chemsi_pm10)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm10)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm10)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm10) - - else if (chemicda_opt == 3) then - write(ounit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_pm25), & - iv%chemic_surf(n)%chem(p_chemsi_pm25)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm25)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm25)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm25), & - ob%chemic_surf(n)%chem(p_chemsi_pm10), & - iv%chemic_surf(n)%chem(p_chemsi_pm10)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm10)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm10)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm10) - + if (chem_cv_options > 1 .and. chem_cv_options < 108) then + write(ounit,'(2i8)') 1, ifgat + if (chemicda_opt == 5) then + write(ounit2,'(2i8)') 1, ifgat + end if + end if ! if (chem_cv_options > 1 .and. chem_cv_options < 108) + + ipcc = p_chemsi_pm25 + jpcc = p_chemsi_pm10 + if (chemicda_opt <= 2) then + if(chemicda_opt == 2) ipcc = p_chemsi_pm10 + write(ounit, fmt = fmt_chem1) & + num_obs , iv%info(chemic_surf)%id(n), & ! Station + iv%info(chemic_surf)%lat(1,n), & ! Latitude + iv%info(chemic_surf)%lon(1,n), & ! Longitude + ob%chemic_surf(n)%chem(ipcc), & ! observation + iv%chemic_surf(n)%chem(ipcc)%inv, & ! o-b + iv%chemic_surf(n)%chem(ipcc)%qc, & ! qc + iv%chemic_surf(n)%chem(ipcc)%error, & ! obs error + re%chemic_surf(n)%chem(ipcc) ! o-a else if (chemicda_opt == 4) then - write(ounit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_so2), & - iv%chemic_surf(n)%chem(p_chemsi_so2)%inv, iv%chemic_surf(n)%chem(p_chemsi_so2)%qc, iv%chemic_surf(n)%chem(p_chemsi_so2)%error, & - re%chemic_surf(n)%chem(p_chemsi_so2), & - ob%chemic_surf(n)%chem(p_chemsi_no2), & - iv%chemic_surf(n)%chem(p_chemsi_no2)%inv, iv%chemic_surf(n)%chem(p_chemsi_no2)%qc, iv%chemic_surf(n)%chem(p_chemsi_no2)%error, & - re%chemic_surf(n)%chem(p_chemsi_no2), & - ob%chemic_surf(n)%chem(p_chemsi_o3), & - iv%chemic_surf(n)%chem(p_chemsi_o3)%inv, iv%chemic_surf(n)%chem(p_chemsi_o3)%qc, iv%chemic_surf(n)%chem(p_chemsi_o3)%error, & - re%chemic_surf(n)%chem(p_chemsi_o3), & - ob%chemic_surf(n)%chem(p_chemsi_co), & - iv%chemic_surf(n)%chem(p_chemsi_co)%inv, iv%chemic_surf(n)%chem(p_chemsi_co)%qc, iv%chemic_surf(n)%chem(p_chemsi_co)%error, & - re%chemic_surf(n)%chem(p_chemsi_co) - - else if (chemicda_opt == 5) then - write(ounit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_pm25), & - iv%chemic_surf(n)%chem(p_chemsi_pm25)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm25)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm25)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm25), & - ob%chemic_surf(n)%chem(p_chemsi_pm10), & - iv%chemic_surf(n)%chem(p_chemsi_pm10)%inv, iv%chemic_surf(n)%chem(p_chemsi_pm10)%qc, iv%chemic_surf(n)%chem(p_chemsi_pm10)%error, & - re%chemic_surf(n)%chem(p_chemsi_pm10) - - write(ounit2,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')& - num_obs , 1, iv%info(chemic_surf)%id(n), & ! Station - iv%info(chemic_surf)%lat(1,n), & ! Latitude - iv%info(chemic_surf)%lon(1,n), & ! Longitude - iv%info(chemic_surf)%lon(1,n), & - ob%chemic_surf(n)%chem(p_chemsi_so2), & - iv%chemic_surf(n)%chem(p_chemsi_so2)%inv, iv%chemic_surf(n)%chem(p_chemsi_so2)%qc, iv%chemic_surf(n)%chem(p_chemsi_so2)%error, & - re%chemic_surf(n)%chem(p_chemsi_so2), & - ob%chemic_surf(n)%chem(p_chemsi_no2), & - iv%chemic_surf(n)%chem(p_chemsi_no2)%inv, iv%chemic_surf(n)%chem(p_chemsi_no2)%qc, iv%chemic_surf(n)%chem(p_chemsi_no2)%error, & - re%chemic_surf(n)%chem(p_chemsi_no2), & - ob%chemic_surf(n)%chem(p_chemsi_o3), & - iv%chemic_surf(n)%chem(p_chemsi_o3)%inv, iv%chemic_surf(n)%chem(p_chemsi_o3)%qc, iv%chemic_surf(n)%chem(p_chemsi_o3)%error, & - re%chemic_surf(n)%chem(p_chemsi_o3), & - ob%chemic_surf(n)%chem(p_chemsi_co), & - iv%chemic_surf(n)%chem(p_chemsi_co)%inv, iv%chemic_surf(n)%chem(p_chemsi_co)%qc, iv%chemic_surf(n)%chem(p_chemsi_co)%error, & - re%chemic_surf(n)%chem(p_chemsi_co) - end if - end if - end do - end if + write(ounit, fmt = fmt_chem4) & + num_obs , iv%info(chemic_surf)%id(n), & ! Station + iv%info(chemic_surf)%lat(1,n), & ! Latitude + iv%info(chemic_surf)%lon(1,n), & ! Longitude + (ob%chemic_surf(n)%chem(ipc(ic)), & ! observation + iv%chemic_surf(n)%chem(ipc(ic))%inv, & ! o-b + iv%chemic_surf(n)%chem(ipc(ic))%qc, & ! qc + iv%chemic_surf(n)%chem(ipc(ic))%error,& ! obs error + re%chemic_surf(n)%chem(ipc(ic)), ic=1,ngas) ! o-a + else ! chemicda_opt == 3 or 5 + write(ounit, fmt = fmt_chem2) & + num_obs , iv%info(chemic_surf)%id(n), & ! Station + iv%info(chemic_surf)%lat(1,n), & ! Latitude + iv%info(chemic_surf)%lon(1,n), & ! Longitude + ob%chemic_surf(n)%chem(ipcc), & ! observation + iv%chemic_surf(n)%chem(ipcc)%inv, & ! o-b + iv%chemic_surf(n)%chem(ipcc)%qc, & ! qc + iv%chemic_surf(n)%chem(ipcc)%error, & ! obs error + re%chemic_surf(n)%chem(ipcc), & ! o-a + ob%chemic_surf(n)%chem(jpcc), & ! observation + iv%chemic_surf(n)%chem(jpcc)%inv, & ! o-b + iv%chemic_surf(n)%chem(jpcc)%qc, & ! qc + iv%chemic_surf(n)%chem(jpcc)%error, & ! obs error + re%chemic_surf(n)%chem(jpcc) ! o-a + if(chemicda_opt == 5) then + write(ounit2, fmt = fmt_chem4) & + num_obs , iv%info(chemic_surf)%id(n), & ! Station + iv%info(chemic_surf)%lat(1,n), & ! Latitude + iv%info(chemic_surf)%lon(1,n), & ! Longitude + (ob%chemic_surf(n)%chem(ipc(ic)), & ! observation + iv%chemic_surf(n)%chem(ipc(ic))%inv, & ! o-b + iv%chemic_surf(n)%chem(ipc(ic))%qc, & ! qc + iv%chemic_surf(n)%chem(ipc(ic))%error,& ! obs error + re%chemic_surf(n)%chem(ipc(ic)), ic=1,ngas) ! o-a + endif ! (chemicda_opt == 5) + endif + + end if ! if (iv%info(chemic_surf)%proc_domain(1,n)) then + end do ! do n = 1, iv%info(chemic_surf)%nlocal + end if ! if (num_obs > 0) then close (ounit) call da_free_unit(ounit) @@ -187,5 +171,3 @@ subroutine da_write_obs_chem_sfc(it,ob, iv, re) if (trace_use) call da_trace_exit("da_write_obs_chem_sfc") end subroutine da_write_obs_chem_sfc - - diff --git a/var/da/da_setup_structures/da_get_vertical_truncation.inc b/var/da/da_setup_structures/da_get_vertical_truncation.inc index 10bc73fff7..7cf737e349 100644 --- a/var/da/da_setup_structures/da_get_vertical_truncation.inc +++ b/var/da/da_setup_structures/da_get_vertical_truncation.inc @@ -56,7 +56,8 @@ subroutine da_get_vertical_truncation( max_vert_var, eigenval, be_sub) do k = 1, kz cum_variance = cum_variance + eigenval(k) - if (eigenval(k) <= 0.0) then + if (eigenval(k) <= 1e-12) then ! Ha: Ignore trivial fields (ex. naaj, naai and seas) + !if (eigenval(k) <= 0.0) then be_sub % mz = k - 1 exit end if @@ -68,6 +69,7 @@ subroutine da_get_vertical_truncation( max_vert_var, eigenval, be_sub) end do if (max_vert_var == 0.0) be_sub % mz = 0 + if (be_sub % mz == 0.0) be_sub % mz = 1 ! HA - to avoid segmentation fault end if @@ -75,7 +77,7 @@ subroutine da_get_vertical_truncation( max_vert_var, eigenval, be_sub) 'Vertical truncation for ', name, & ' = ', be_sub % mz, '(', & max_vert_var, '%)' - write (unit=stdout,fmt='(A)') " " + !write (unit=stdout,fmt='(A)') " " if (trace_use_dull) call da_trace_exit("da_get_vertical_truncation") diff --git a/var/da/da_setup_structures/da_setup_background_errors.inc b/var/da/da_setup_structures/da_setup_background_errors.inc index 8cebb1e3ac..1913483890 100644 --- a/var/da/da_setup_structures/da_setup_background_errors.inc +++ b/var/da/da_setup_structures/da_setup_background_errors.inc @@ -77,22 +77,11 @@ subroutine da_setup_background_errors(grid, be) be % v11% mz = 0 #if (WRF_CHEM == 1) - num_cv_3d_basic = 0 - num_cv_3d_extra = 0 - num_cv_2d = 0 if (num_chem >= PARAM_FIRST_SCALAR) then -! allocate (be%v12 (num_chem) ) -! do ic = 1, num_chem -! be % v12(ic) % mz = 0 -! end do allocate (be%cv%sizechemic (num_chem) ) end if - be % ncv_mz_chemic = num_cv_3d_basic + num_cv_3d_extra + num_cv_2d - allocate ( be%cv_mz_chemic(1:be%ncv_mz_chemic) ) - be%cv_mz_chemic(:) = 0 ! initialize - #endif if (global) then diff --git a/var/da/da_setup_structures/da_setup_be_regional.inc b/var/da/da_setup_structures/da_setup_be_regional.inc index 64dcb2387a..62999fbddb 100644 --- a/var/da/da_setup_structures/da_setup_be_regional.inc +++ b/var/da/da_setup_structures/da_setup_be_regional.inc @@ -24,6 +24,18 @@ subroutine da_setup_be_regional(xb, be, grid) ,p_cl_a01, p_cl_a02, p_cl_a03, p_cl_a04 & ,p_na_a01, p_na_a02, p_na_a03, p_na_a04 & ,p_oin_a01, p_oin_a02, p_oin_a03, p_oin_a04 + + USE module_state_description, ONLY : p_so4aj, p_so4ai, p_nh4aj, p_nh4ai, & ! racm_soa_vbs_da + p_no3aj, p_no3ai, p_naaj, p_naai, & + p_asoa1j, p_asoa1i, p_asoa2j, p_asoa2i, & + p_asoa3j, p_asoa3i, p_asoa4j, p_asoa4i, & + p_bsoa1j, p_bsoa1i, p_bsoa2j, p_bsoa2i, & + p_bsoa3j, p_bsoa3i, p_bsoa4j, p_bsoa4i, & + p_orgpaj, p_orgpai, p_ecj, p_eci, & + p_p25j, p_p25i, p_antha, p_seas, & + p_claj, p_clai, p_soila +! p_so2, p_no2, p_o3, p_co ! gas species + #endif implicit none @@ -285,15 +297,20 @@ subroutine da_setup_be_regional(xb, be, grid) num_cv_3d_extra = 5 ! 5 hydrometeors, qc,qr,qi,qs,qg end if #if (WRF_CHEM == 1) - if ( chem_cv_options == 10 ) then - num_cv_3d_extra = num_cv_3d_extra + 19 + select case ( chem_cv_options ) + case ( 10 ) num_cv_3d_chem = 19 - end if - if ( chem_cv_options == 20 ) then - num_cv_3d_extra = num_cv_3d_extra + 36 + case ( 20 ) num_cv_3d_chem = 36 - end if - allocate (be % v12 (num_cv_3d_chem) ) + case ( 108 ) ! racm_soa_vbs_da + num_cv_3d_chem = 39 + case default ! mozcart_da + num_cv_3d_chem = 19 + end select + if ( chem_cv_options >= 1 ) then + allocate (be % v12 (num_cv_3d_chem) ) + num_cv_3d_extra = num_cv_3d_extra + num_cv_3d_chem + endif ! ( chem_cv_options >= 1 ) #endif if ( use_cv_w ) then num_cv_3d_extra = num_cv_3d_extra + 1 ! clouds+w @@ -888,8 +905,512 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('dust_3' ) - ispcc=p_dust_3-1 + case ('dust_3' ) + ispcc=p_dust_3-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('dust_4' ) + ispcc=p_dust_4-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('seas_1' ) + ispcc=p_seas_1-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('seas_2' ) + ispcc=p_seas_2-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('seas_3' ) + ispcc=p_seas_3-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('seas_4' ) + ispcc=p_seas_4-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc1' ) + ispcc=p_oc1-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc2' ) + ispcc=p_oc2-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('bc1' ) + ispcc=p_bc1-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('bc2' ) + ispcc=p_bc2-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('so2' ) + ispcc=p_so2-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('no2' ) + ispcc=p_no2-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('o3' ) + ispcc=p_o3-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('co' ) + ispcc=p_co-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('sulf' ) + ispcc=p_sulf-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('p25' ) + ispcc=p_p25-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('p10' ) + ispcc=p_p10-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('w' ) + !if ( use_cv_w ) then + ! if ( trim(adjustl(variable)) == 'w' ) then + be % v11 % name = trim(adjustl(variable)) + read (be_unit) be11_evec_glo + read (be_unit) be11_eval_glo + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b) + end do + ! end if + !end if + + case default; + message(1)=' Read problem in eigen vectors/values in BE file ' + write (unit=message(2),fmt='(A,A)') ' Trying to read Eigenvectors for variable: ',trim(adjustl(variable)) + write (unit=message(3),fmt='(A)') ' Make sure you are using the correct be.dat file for your cv_options setting!' + call da_error(__FILE__,__LINE__,message(1:3)) + end select + + end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1 + print*,"reading variables for chem_cv_options=10 is done" + end if ! chem_cv_options=10 + + + if ( chem_cv_options == 20 ) then + print*,"chem_cv_options == 20 variables from id ",num_cv_3d_basic+1 ,"to id", num_cv_3d_basic+num_cv_3d_extra + do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra + read (be_unit) variable + read (be_unit) nk, num_bins2d + print *, trim(adjustl(variable)), nk, num_bins2d + select case( trim(adjustl(variable)) ) + + case ('bc_1' ) + ispcc=p_bc_a01-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('bc_2' ) + ispcc=p_bc_a02-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('bc_3' ) + ispcc=p_bc_a03-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('bc_4' ) + ispcc=p_bc_a04-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc_1' ) + ispcc=p_oc_a01-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc_2' ) + ispcc=p_oc_a02-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc_3' ) + ispcc=p_oc_a03-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('oc_4' ) + ispcc=p_oc_a04-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('so4_1' ) + ispcc=p_so4_a01-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('so4_2' ) + ispcc=p_so4_a02-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('so4_3' ) + ispcc=p_so4_a03-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('so4_4' ) + ispcc=p_so4_a04-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('no3_1' ) + ispcc=p_no3_a01-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('no3_2' ) + ispcc=p_no3_a02-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('no3_3' ) + ispcc=p_no3_a03-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('no3_4' ) + ispcc=p_no3_a04-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('nh4_1' ) + ispcc=p_nh4_a01-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('nh4_2' ) + ispcc=p_nh4_a02-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('nh4_3' ) + ispcc=p_nh4_a03-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('nh4_4' ) + ispcc=p_nh4_a04-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -901,8 +1422,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('dust_4' ) - ispcc=p_dust_4-1 + case ('cl_1' ) + ispcc=p_cl_a01-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -914,8 +1435,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('seas_1' ) - ispcc=p_seas_1-1 + case ('cl_2' ) + ispcc=p_cl_a02-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -927,8 +1448,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('seas_2' ) - ispcc=p_seas_2-1 + case ('cl_3' ) + ispcc=p_cl_a03-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -940,8 +1461,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('seas_3' ) - ispcc=p_seas_3-1 + case ('cl_4' ) + ispcc=p_cl_a04-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -953,8 +1474,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('seas_4' ) - ispcc=p_seas_4-1 + case ('na_1' ) + ispcc=p_na_a01-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -966,8 +1487,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc1' ) - ispcc=p_oc1-1 + case ('na_2' ) + ispcc=p_na_a02-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -979,8 +1500,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc2' ) - ispcc=p_oc2-1 + case ('na_3' ) + ispcc=p_na_a03-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -992,8 +1513,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('bc1' ) - ispcc=p_bc1-1 + case ('na_4' ) + ispcc=p_na_a04-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1005,8 +1526,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('bc2' ) - ispcc=p_bc2-1 + case ('oin_1' ) + ispcc=p_oin_a01-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1018,8 +1539,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('so2' ) - ispcc=p_so2-1 + case ('oin_2' ) + ispcc=p_oin_a02-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1031,8 +1552,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('no2' ) - ispcc=p_no2-1 + case ('oin_3' ) + ispcc=p_oin_a03-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1044,8 +1565,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('o3' ) - ispcc=p_o3-1 + case ('oin_4' ) + ispcc=p_oin_a04-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1057,8 +1578,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('co' ) - ispcc=p_co-1 + case ('so2' ) + ispcc=p_so2-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1070,8 +1591,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('sulf' ) - ispcc=p_sulf-1 + case ('no2' ) + ispcc=p_no2-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1083,8 +1604,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('p25' ) - ispcc=p_p25-1 + case ('o3' ) + ispcc=p_o3-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1096,8 +1617,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('p10' ) - ispcc=p_p10-1 + case ('co' ) + ispcc=p_co-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1110,8 +1631,7 @@ subroutine da_setup_be_regional(xb, be, grid) end do case ('w' ) - !if ( use_cv_w ) then - ! if ( trim(adjustl(variable)) == 'w' ) then + if ( use_cv_w ) then be % v11 % name = trim(adjustl(variable)) read (be_unit) be11_evec_glo read (be_unit) be11_eval_glo @@ -1122,8 +1642,7 @@ subroutine da_setup_be_regional(xb, be, grid) be11_evec_loc(j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) be11_eval_loc(j,1:nk ) = eval_loc(1:nk,b) end do - ! end if - !end if + end if case default; message(1)=' Read problem in eigen vectors/values in BE file ' @@ -1133,20 +1652,21 @@ subroutine da_setup_be_regional(xb, be, grid) end select end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1 - print*,"reading variables for chem_cv_options=10 is done" - end if ! chem_cv_options=10 - + print*,"reading variables for chem_cv_options=20 is done" + end if ! chem_cv_options=20 - if ( chem_cv_options == 20 ) then - print*,"chem_cv_options == 20 variables from id ",num_cv_3d_basic+1 ,"to id", num_cv_3d_basic+num_cv_3d_extra + if ( chem_cv_options == 108 ) then ! racm_soa_vbs_da + print*,"chem_cv_options == 108 variables from id ",num_cv_3d_basic+1 ,"to id", & + num_cv_3d_basic+num_cv_3d_extra do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra read (be_unit) variable read (be_unit) nk, num_bins2d - print *, trim(adjustl(variable)), nk, num_bins2d + !print *, trim(adjustl(variable)), nk, num_bins2d + select case( trim(adjustl(variable)) ) - case ('bc_1' ) - ispcc=p_bc_a01-1 + case ('so4aj' ) + ispcc=p_so4aj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1154,12 +1674,12 @@ subroutine da_setup_be_regional(xb, be, grid) read (be_unit) eval_loc do j=1,nj b = bin2d(1,j) - be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('bc_2' ) - ispcc=p_bc_a02-1 + case ('so4ai' ) + ispcc=p_so4ai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1171,8 +1691,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('bc_3' ) - ispcc=p_bc_a03-1 + case ('nh4aj' ) + ispcc=p_nh4aj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1184,8 +1704,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('bc_4' ) - ispcc=p_bc_a04-1 + case ('nh4ai' ) + ispcc=p_nh4ai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1197,8 +1717,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc_1' ) - ispcc=p_oc_a01-1 + case ('no3aj' ) + ispcc=p_no3aj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1210,8 +1730,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc_2' ) - ispcc=p_oc_a02-1 + case ('no3ai' ) + ispcc=p_no3ai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1223,8 +1743,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc_3' ) - ispcc=p_oc_a03-1 + case ('naaj' ) + ispcc=p_naaj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1236,8 +1756,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oc_4' ) - ispcc=p_oc_a04-1 + case ('naai' ) + ispcc=p_naai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1249,8 +1769,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('so4_1' ) - ispcc=p_so4_a01-1 + case ('claj' ) + ispcc=p_claj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1262,8 +1782,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('so4_2' ) - ispcc=p_so4_a02-1 + case ('clai' ) + ispcc=p_clai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1275,8 +1795,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('so4_3' ) - ispcc=p_so4_a03-1 + case ('asoa1j' ) + ispcc=p_asoa1j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1288,8 +1808,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('so4_4' ) - ispcc=p_so4_a04-1 + case ('asoa1i' ) + ispcc=p_asoa1i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1301,8 +1821,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('no3_1' ) - ispcc=p_no3_a01-1 + case ('asoa2j' ) + ispcc=p_asoa2j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1314,8 +1834,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('no3_2' ) - ispcc=p_no3_a02-1 + case ('asoa2i' ) + ispcc=p_asoa2i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1327,8 +1847,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('no3_3' ) - ispcc=p_no3_a03-1 + case ('asoa3j' ) + ispcc=p_asoa3j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1340,8 +1860,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('no3_4' ) - ispcc=p_no3_a04-1 + case ('asoa3i' ) + ispcc=p_asoa3i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1353,8 +1873,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('nh4_1' ) - ispcc=p_nh4_a01-1 + case ('asoa4j' ) + ispcc=p_asoa4j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1366,8 +1886,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('nh4_2' ) - ispcc=p_nh4_a02-1 + case ('asoa4i' ) + ispcc=p_asoa4i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1379,8 +1899,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('nh4_3' ) - ispcc=p_nh4_a03-1 + case ('bsoa1j' ) + ispcc=p_bsoa1j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1392,8 +1912,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('nh4_4' ) - ispcc=p_nh4_a04-1 + case ('bsoa1i' ) + ispcc=p_bsoa1i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1405,8 +1925,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('cl_1' ) - ispcc=p_cl_a01-1 + case ('bsoa2j' ) + ispcc=p_bsoa2j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1418,8 +1938,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('cl_2' ) - ispcc=p_cl_a02-1 + case ('bsoa2i' ) + ispcc=p_bsoa2i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1431,8 +1951,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('cl_3' ) - ispcc=p_cl_a03-1 + case ('bsoa3j' ) + ispcc=p_bsoa3j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1444,8 +1964,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('cl_4' ) - ispcc=p_cl_a04-1 + case ('bsoa3i' ) + ispcc=p_bsoa3i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1457,8 +1977,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('na_1' ) - ispcc=p_na_a01-1 + case ('bsoa4j' ) + ispcc=p_bsoa4j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1470,8 +1990,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('na_2' ) - ispcc=p_na_a02-1 + case ('bsoa4i' ) + ispcc=p_bsoa4i-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1483,8 +2003,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('na_3' ) - ispcc=p_na_a03-1 + case ('orgpaj' ) + ispcc=p_orgpaj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1496,8 +2016,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('na_4' ) - ispcc=p_na_a04-1 + case ('orgpai' ) + ispcc=p_orgpai-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1509,8 +2029,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oin_1' ) - ispcc=p_oin_a01-1 + case ('ecj' ) + ispcc=p_ecj-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1522,8 +2042,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oin_2' ) - ispcc=p_oin_a02-1 + case ('eci' ) + ispcc=p_eci-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1535,8 +2055,8 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oin_3' ) - ispcc=p_oin_a03-1 + case ('p25j' ) + ispcc=p_p25j-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1548,8 +2068,47 @@ subroutine da_setup_be_regional(xb, be, grid) be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) end do - case ('oin_4' ) - ispcc=p_oin_a04-1 + case ('p25i' ) + ispcc=p_p25i-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('antha' ) + ispcc=p_antha-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('seas' ) + ispcc=p_seas-1 + be % v12(ispcc) % name = trim(adjustl(variable)) + read (be_unit) be12_evec_glo(ispcc,:,:) + read (be_unit) be12_eval_glo(ispcc,:) + read (be_unit) evec_loc + read (be_unit) eval_loc + do j=1,nj + b = bin2d(1,j) + be12_evec_loc(ispcc,j,1:nk,1:nk) = evec_loc(1:nk,1:nk,b) + be12_eval_loc(ispcc,j,1:nk ) = eval_loc(1:nk,b) + end do + + case ('soila' ) + ispcc=p_soila-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1562,7 +2121,7 @@ subroutine da_setup_be_regional(xb, be, grid) end do case ('so2' ) - ispcc=p_so2-1 + ispcc=p_so2-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1575,7 +2134,7 @@ subroutine da_setup_be_regional(xb, be, grid) end do case ('no2' ) - ispcc=p_no2-1 + ispcc=p_no2-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1588,7 +2147,7 @@ subroutine da_setup_be_regional(xb, be, grid) end do case ('o3' ) - ispcc=p_o3-1 + ispcc=p_o3-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1601,7 +2160,7 @@ subroutine da_setup_be_regional(xb, be, grid) end do case ('co' ) - ispcc=p_co-1 + ispcc=p_co-1 be % v12(ispcc) % name = trim(adjustl(variable)) read (be_unit) be12_evec_glo(ispcc,:,:) read (be_unit) be12_eval_glo(ispcc,:) @@ -1634,9 +2193,10 @@ subroutine da_setup_be_regional(xb, be, grid) call da_error(__FILE__,__LINE__,message(1:3)) end select + !print '(A6,3(A,I4),I5)', trim(adjustl(variable)), ' i =', i, ' , ispcc =', ispcc,' , nk =', nk, num_bins2d end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra-1 - print*,"reading variables for chem_cv_options=20 is done" - end if ! chem_cv_options=20 + print*,"reading variables for chem_cv_options=108 (racm_soa_vbs_da) is done" + end if ! ( chem_cv_options == 108 ) then ! racm_soa_vbs_da #endif deallocate (evec_loc) @@ -1684,7 +2244,7 @@ subroutine da_setup_be_regional(xb, be, grid) message(1:5) = '' write (unit=message(1),fmt='(3x,A,i2)') 'cloud_cv_options = ', cloud_cv_options #if (WRF_CHEM == 1) - write (unit=message(1),fmt='(3x,A,i2)') 'chem_cv_options = ', chem_cv_options + write (unit=message(1),fmt='(3x,A,i3)') 'chem_cv_options = ', chem_cv_options #endif write (unit=message(2),fmt='(3x,A,7A)') 'WRFDA dry control variables are: ', & trim(be % v1 % name), ', ', trim(be % v2 % name), ', ', & @@ -2006,7 +2566,141 @@ subroutine da_setup_be_regional(xb, be, grid) call da_error(__FILE__,__LINE__,message(1:2)) end select end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra - end if + end if ! if ( chem_cv_options == 20 ) then + + if ( chem_cv_options == 108 ) then ! racm_soa_vbs_da + do i = num_cv_3d_basic+1 , num_cv_3d_basic+num_cv_3d_extra + read (be_unit) variable + print *, trim(adjustl(variable)) + select case( trim(adjustl(variable)) ) + case ('so4aj') + ispcc=p_so4aj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('so4ai' ) + ispcc=p_so4ai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('nh4aj' ) + ispcc=p_nh4aj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('nh4ai' ) + ispcc=p_nh4ai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('no3aj' ) + ispcc=p_no3aj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('no3ai' ) + ispcc=p_no3ai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('naaj' ) + ispcc=p_naaj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('naai' ) + ispcc=p_naai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('claj' ) + ispcc=p_claj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('clai' ) + ispcc=p_clai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa1j' ) + ispcc=p_asoa1j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa1i' ) + ispcc=p_asoa1i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa2j' ) + ispcc=p_asoa2j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa2i' ) + ispcc=p_asoa2i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa3j' ) + ispcc=p_asoa3j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa3i' ) + ispcc=p_asoa3i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa4j' ) + ispcc=p_asoa4j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('asoa4i' ) + ispcc=p_asoa4i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa1j' ) + ispcc=p_bsoa1j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa1i' ) + ispcc=p_bsoa1i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa2j' ) + ispcc=p_bsoa2j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa2i' ) + ispcc=p_bsoa2i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa3j' ) + ispcc=p_bsoa3j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa3i' ) + ispcc=p_bsoa3i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa4j' ) + ispcc=p_bsoa4j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('bsoa4i' ) + ispcc=p_bsoa4i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('orgpaj' ) + ispcc=p_orgpaj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('orgpai' ) + ispcc=p_orgpai-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('ecj' ) + ispcc=p_ecj-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('eci' ) + ispcc=p_eci-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('p25j' ) + ispcc=p_p25j-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('p25i' ) + ispcc=p_p25i-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('antha' ) + ispcc=p_antha-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('seas' ) + ispcc=p_seas-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('soila' ) + ispcc=p_soila-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('so2' ) + ispcc=p_so2-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('no2' ) + ispcc=p_no2-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('o3' ) + ispcc=p_o3-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('co' ) + ispcc=p_co-1 + read(be_unit) be12_rf_lengthscale(ispcc,:) + case ('w') + if ( use_cv_w ) then + read(be_unit) be11_rf_lengthscale + end if + case default; + message(1)='Read problem in lengthscales in be.dat' + write(message(2),'("Trying to read lengthscales for variable ",I0,": ",A)')i,trim(adjustl(variable)) + call da_error(__FILE__,__LINE__,message(1:2)) + end select + end do ! num_cv_3d_basic+1 - num_cv_3d_basic+num_cv_3d_extra + endif !( chem_cv_options == 108 ) then ! racm_soa_vbs_da #endif ! Read in lengthscale of 2D Control variable ps_u @@ -2059,7 +2753,7 @@ subroutine da_setup_be_regional(xb, be, grid) !#if (WRF_CHEM == 1) ! if ( chem_cv_options >= 10 ) then ! do i=1,num_cv_3d_chem -! write(unit=message(i+6),fmt=write_fmt) 'be',(i+5),'= ',be12_rf_lengthscale(i,:) +! write(unit=message(i+6),fmt=write_fmt) 'be12 for ichem ',(i+5),'= ',be12_rf_lengthscale(i,:) ! end do ! call da_message(message(i+6:num_cv_3d_chem+6)) ! end if @@ -3023,7 +3717,12 @@ subroutine da_setup_be_regional(xb, be, grid) end if #endif if ( use_cv_w ) then - be%cv_mz(num_cv_3d_basic+num_cv_3d_extra+num_cv_2d) = be%v11%mz + if (be%ncv_mz ==8 .or. be%ncv_mz ==13 ) then + be%cv_mz(be%ncv_mz-2) = be%v11%mz + else + be%cv_mz(num_cv_3d_basic+num_cv_3d_extra+num_cv_2d) = be%v11%mz + end if + !be%cv_mz(num_cv_3d_basic+num_cv_3d_extra+num_cv_2d) = be%v11%mz ! This is a bug. (Soyoung on Aug-2021) end if be%cv_mz(be%ncv_mz-1) = be%alpha%mz be%cv_mz(be%ncv_mz) = be%ne diff --git a/var/da/da_setup_structures/da_setup_firstguess_wrf.inc b/var/da/da_setup_structures/da_setup_firstguess_wrf.inc index cff3baf9b7..28a0d0f2ba 100644 --- a/var/da/da_setup_structures/da_setup_firstguess_wrf.inc +++ b/var/da/da_setup_structures/da_setup_firstguess_wrf.inc @@ -4,14 +4,6 @@ subroutine da_setup_firstguess_wrf(xbx, grid, config_flags, ens) ! Purpose: Define/allocate components of WRF model state. !--------------------------------------------------------------------------- -#if (WRF_CHEM == 1) - USE module_state_description, ONLY : & - num_chem & - ,p_chem_ic_bc1, p_chem_ic_bc2, p_chem_ic_oc1, p_chem_ic_oc2 & - ,p_chem_ic_so2, p_chem_ic_no2, p_chem_ic_o3, p_chem_ic_co & - ,p_chem_ic_p10, p_chem_ic_p25, p_chem_ic_sulf -#endif - implicit none type (xbx_type), intent(out) :: xbx ! Header & non-gridded vars. diff --git a/var/da/da_setup_structures/da_setup_obs_structures.inc b/var/da/da_setup_structures/da_setup_obs_structures.inc index da28434f70..97a4909f69 100644 --- a/var/da/da_setup_structures/da_setup_obs_structures.inc +++ b/var/da/da_setup_structures/da_setup_obs_structures.inc @@ -323,13 +323,9 @@ subroutine da_setup_obs_structures( grid, ob, iv, j_cost) write(unit=stdout, fmt='(a)') 'Observation summary' do i=1,num_fgat_time - write(unit=stdout, fmt='(3x,a,i2)') 'ob time ', i + write(unit=stdout, fmt='(3x,2(a,i3))') 'ob time', i,' num_ob_indexes=',num_ob_indexes do j=1,num_ob_indexes - if(j.eq.27 & -#if (WRF_CHEM == 1) - .or. j.eq.29 .or. j.eq.30 & -#endif - )cycle + if(j.eq.27) cycle ! skip tamdar_sfc (HA) if (use_gpsztdobs .and.obs_names(j) == 'gpspw ' ) then ob_name = 'gpsztd ' else @@ -344,10 +340,10 @@ subroutine da_setup_obs_structures( grid, ob, iv, j_cost) iv%info(j)%thin_plocal(i) - iv%info(j)%thin_plocal(i-1), ' local (post-thinning)' else #if (WRF_CHEM == 1) - write(unit=stdout, fmt='(6x,a,i10,a,i8,a,i8,a,i8,a)') & - ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', & + write(unit=stdout, fmt='(i3,3x,a,i10,a,i8,a,i8,a,i8,a)') & + j,ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', & iv%info(j)%plocal(i) - iv%info(j)%plocal(i-1), ' local', & - iv%info(chemic_surf)%ntotal, 'chem total', iv%info(chemic_surf)%nlocal, 'chem local' + iv%info(chemic_surf)%ntotal, ' chem total', iv%info(chemic_surf)%nlocal, ' chem local' #else write(unit=stdout, fmt='(6x,a,i10,a,i8,a)') & ob_name, iv%info(j)%ptotal(i) - iv%info(j)%ptotal(i-1), ' global,', & diff --git a/var/da/da_setup_structures/da_setup_obs_structures_chem_sfc.inc b/var/da/da_setup_structures/da_setup_obs_structures_chem_sfc.inc index b09036df51..846f251745 100644 --- a/var/da/da_setup_structures/da_setup_obs_structures_chem_sfc.inc +++ b/var/da/da_setup_structures/da_setup_obs_structures_chem_sfc.inc @@ -25,17 +25,15 @@ subroutine da_setup_obs_structures_chem_sfc(ob, iv, grid) !-------------------------------------------------------------------------- ! [1.0] Scan data - only necessary if an alternative format is implemented !-------------------------------------------------------------------------- - call da_scan_obs_chem_sfc(iv, filename, grid) - + iv%time = 1 !!! Fixed by Soyoung Ha (June-2021) iv%info(chemic_surf)%plocal(iv%time) = iv%info(chemic_surf)%nlocal iv%info(chemic_surf)%ptotal(iv%time) = iv%info(chemic_surf)%ntotal !-------------------------------------------------------------------------- ! Allocate the ob based on input number of obs: !-------------------------------------------------------------------------- - call da_allocate_observations_chem_sfc (iv) !-------------------------------------------------------------------------- @@ -46,7 +44,6 @@ subroutine da_setup_obs_structures_chem_sfc(ob, iv, grid) !-------------------------------------------------------------------------- ! [3.0] Calculate innovation vector (O-B) and create (smaller) ob structure: !-------------------------------------------------------------------------- - call da_fill_obs_structures_chem_sfc(iv, ob) if (trace_use) call da_trace_exit("da_setup_obs_structures_chem_sfc") diff --git a/var/da/da_statistics/da_analysis_stats.inc b/var/da/da_statistics/da_analysis_stats.inc index b3ff3be384..672946b14e 100644 --- a/var/da/da_statistics/da_analysis_stats.inc +++ b/var/da/da_statistics/da_analysis_stats.inc @@ -120,6 +120,17 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR 'dust_2 i j ', & 'seas_1 i j ', & 'seas_2 i j' + case ( 108 ) + write(unit=stats_unit2, fmt='(9a/)') & + ' Lvl, so4aj i j, so4ai i j, nh4aj i j, nh4ai i j, no3aj i j,', & + 'no3ai i j, naaj i j, naai i j, claj i j, clai i j,', & + 'asoa1j i j, asoa1i i j, asoa2j i j, asoa2i i j, ', & + 'asoa3j i j, asoa3i i j, asoa4j i j, asoa4i i j,', & + 'bsoa1j i j, bsoa1i i j, bsoa2j i j, bsoa2i i j, ', & + 'bsoa3j i j, bsoa3i i j, bsoa4j i j, bsoa4i i j,', & + 'orgpaj i j, orgpai i j, ecj i j, eci i j,', & + 'p25j i j, p25i i j, antha i j, seas i j, soila i j,', & + 'so2 i j, no2 i j, o3 i j, co i j' end select #endif end if @@ -201,7 +212,7 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR min_chem(k, p_chem_ic_bc1), min_chem(k, p_chem_ic_bc2), min_chem(k, p_chem_ic_oc1), min_chem(k, p_chem_ic_oc2), & min_chem(k, p_chem_ic_dust_1), min_chem(k, p_chem_ic_dust_2), min_chem(k, p_chem_ic_seas_1), min_chem(k, p_chem_ic_seas_1) case ( 20 ) - write(unit=stats_unit2, fmt='(i4,24(f12.4,2i5)))') k, & + write(unit=stats_unit2, fmt='(i4,24(f12.4,2i5))') k, & min_chem(k, p_chem_ic_bc_a01), min_chem(k, p_chem_ic_bc_a02), min_chem(k, p_chem_ic_bc_a03), & min_chem(k, p_chem_ic_oc_a01), min_chem(k, p_chem_ic_oc_a02), min_chem(k, p_chem_ic_oc_a03), & min_chem(k, p_chem_ic_so4_a01), min_chem(k, p_chem_ic_so4_a02), min_chem(k, p_chem_ic_so4_a03), & @@ -210,6 +221,8 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR min_chem(k, p_chem_ic_cl_a01), min_chem(k, p_chem_ic_cl_a02), min_chem(k, p_chem_ic_cl_a03), & min_chem(k, p_chem_ic_na_a01), min_chem(k, p_chem_ic_na_a02), min_chem(k, p_chem_ic_na_a03), & min_chem(k, p_chem_ic_oin_a01), min_chem(k, p_chem_ic_oin_a02), min_chem(k, p_chem_ic_oin_a03) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,i4,39(f8.3,2i4))') 'k=',k, (min_chem(k,ic),ic=PARAM_FIRST_SCALAR,num_chem) end select #endif end if ! rootproc @@ -265,6 +278,8 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR chemm(p_chem_ic_cl_a01), chemm(p_chem_ic_cl_a02), chemm(p_chem_ic_cl_a03), & chemm(p_chem_ic_na_a01), chemm(p_chem_ic_na_a02), chemm(p_chem_ic_na_a03), & chemm(p_chem_ic_oin_a01), chemm(p_chem_ic_oin_a02), chemm(p_chem_ic_oin_a03) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,39f8.3)') 'Column_Min', (chemm(ic),ic=PARAM_FIRST_SCALAR,num_chem) end select #endif @@ -359,6 +374,8 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR max_chem(k, p_chem_ic_cl_a01), max_chem(k, p_chem_ic_cl_a02), max_chem(k, p_chem_ic_cl_a03), & max_chem(k, p_chem_ic_na_a01), max_chem(k, p_chem_ic_na_a02), max_chem(k, p_chem_ic_na_a03), & max_chem(k, p_chem_ic_oin_a01), max_chem(k, p_chem_ic_oin_a02), max_chem(k, p_chem_ic_oin_a03) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,i4,39(f8.3,2i4))') 'k=',k, (max_chem(k,ic),ic=PARAM_FIRST_SCALAR,num_chem) end select #endif end if @@ -414,6 +431,8 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR chemm(p_chem_ic_cl_a01), chemm(p_chem_ic_cl_a02), chemm(p_chem_ic_cl_a03), & chemm(p_chem_ic_na_a01), chemm(p_chem_ic_na_a02), chemm(p_chem_ic_na_a03), & chemm(p_chem_ic_oin_a01), chemm(p_chem_ic_oin_a02), chemm(p_chem_ic_oin_a03) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,39f8.3)') 'Column_Max', (chemm(ic),ic=PARAM_FIRST_SCALAR,num_chem) end select #endif @@ -533,7 +552,11 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR chemv(k, p_chem_ic_cl_a01)*rij_g, chemv(k, p_chem_ic_cl_a02)*rij_g, chemv(k, p_chem_ic_cl_a03)*rij_g, & chemv(k, p_chem_ic_na_a01)*rij_g, chemv(k, p_chem_ic_na_a02)*rij_g, chemv(k, p_chem_ic_na_a03)*rij_g, & chemv(k, p_chem_ic_oin_a01)*rij_g, chemv(k, p_chem_ic_oin_a02)*rij_g, chemv(k, p_chem_ic_oin_a03)*rij_g + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,i4,39f8.3)') 'k=',k, & + (chemv(k,ic)*rij_g, ic=PARAM_FIRST_SCALAR,num_chem) end select + do ic=PARAM_FIRST_SCALAR, num_chem chemm(ic)=chemm(ic)+chemv(k,ic) end do @@ -572,23 +595,25 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR end select #if (WRF_CHEM == 1) - select case ( chem_cv_options ) + select case ( chem_cv_options ) case ( 10 ) write(unit=stats_unit2, fmt='(a,10e14.4)') ' ALL', & - chemv(k, p_chem_ic_p25)*rijk_g, chemv(k, p_chem_ic_sulf)*rijk_g, & - chemv(k, p_chem_ic_bc1)*rijk_g, chemv(k, p_chem_ic_bc2)*rijk_g, chemv(k, p_chem_ic_oc1)*rijk_g, chemv(k, p_chem_ic_oc2)*rijk_g, & - chemv(k, p_chem_ic_dust_1)*rijk_g, chemv(k, p_chem_ic_dust_2)*rijk_g, chemv(k, p_chem_ic_seas_2)*rijk_g, chemv(k, p_chem_ic_seas_2)*rijk_g + chemm(p_chem_ic_p25)*rijk_g, chemm(p_chem_ic_sulf)*rijk_g, & + chemm(p_chem_ic_bc1)*rijk_g, chemm(p_chem_ic_bc2)*rijk_g, chemm(p_chem_ic_oc1)*rijk_g, chemm(p_chem_ic_oc2)*rijk_g, & + chemm(p_chem_ic_dust_1)*rijk_g, chemm(p_chem_ic_dust_2)*rijk_g, chemm(p_chem_ic_seas_2)*rijk_g, chemm(p_chem_ic_seas_2)*rijk_g case ( 20 ) write(unit=stats_unit2, fmt='(a,24e14.4)') ' ALL', & - chemv(k, p_chem_ic_bc_a01)*rijk_g, chemv(k, p_chem_ic_bc_a02)*rijk_g, chemv(k, p_chem_ic_bc_a03)*rijk_g, & - chemv(k, p_chem_ic_oc_a01)*rijk_g, chemv(k, p_chem_ic_oc_a02)*rijk_g, chemv(k, p_chem_ic_oc_a03)*rijk_g, & - chemv(k, p_chem_ic_so4_a01)*rijk_g, chemv(k, p_chem_ic_so4_a02)*rijk_g, chemv(k, p_chem_ic_so4_a03)*rijk_g, & - chemv(k, p_chem_ic_no3_a01)*rijk_g, chemv(k, p_chem_ic_no3_a02)*rijk_g, chemv(k, p_chem_ic_no3_a03)*rijk_g, & - chemv(k, p_chem_ic_nh4_a01)*rijk_g, chemv(k, p_chem_ic_nh4_a02)*rijk_g, chemv(k, p_chem_ic_nh4_a03)*rijk_g, & - chemv(k, p_chem_ic_cl_a01)*rijk_g, chemv(k, p_chem_ic_cl_a02)*rijk_g, chemv(k, p_chem_ic_cl_a03)*rijk_g, & - chemv(k, p_chem_ic_na_a01)*rijk_g, chemv(k, p_chem_ic_na_a02)*rijk_g, chemv(k, p_chem_ic_na_a03)*rijk_g, & - chemv(k, p_chem_ic_oin_a01)*rijk_g, chemv(k, p_chem_ic_oin_a02)*rijk_g, chemv(k, p_chem_ic_oin_a03)*rijk_g - end select + chemm(p_chem_ic_bc_a01)*rijk_g, chemm(p_chem_ic_bc_a02)*rijk_g, chemm(p_chem_ic_bc_a03)*rijk_g, & + chemm(p_chem_ic_oc_a01)*rijk_g, chemm(p_chem_ic_oc_a02)*rijk_g, chemm(p_chem_ic_oc_a03)*rijk_g, & + chemm(p_chem_ic_so4_a01)*rijk_g, chemm(p_chem_ic_so4_a02)*rijk_g, chemm(p_chem_ic_so4_a03)*rijk_g, & + chemm(p_chem_ic_no3_a01)*rijk_g, chemm(p_chem_ic_no3_a02)*rijk_g, chemm(p_chem_ic_no3_a03)*rijk_g, & + chemm(p_chem_ic_nh4_a01)*rijk_g, chemm(p_chem_ic_nh4_a02)*rijk_g, chemm(p_chem_ic_nh4_a03)*rijk_g, & + chemm(p_chem_ic_cl_a01)*rijk_g, chemm(p_chem_ic_cl_a02)*rijk_g, chemm(p_chem_ic_cl_a03)*rijk_g, & + chemm(p_chem_ic_na_a01)*rijk_g, chemm(p_chem_ic_na_a02)*rijk_g, chemm(p_chem_ic_na_a03)*rijk_g, & + chemm(p_chem_ic_oin_a01)*rijk_g, chemm(p_chem_ic_oin_a02)*rijk_g, chemm(p_chem_ic_oin_a03)*rijk_g + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,39f8.3)') 'ALL', (chemm(ic)*rijk_g,ic=PARAM_FIRST_SCALAR,num_chem) + end select #endif @@ -758,6 +783,9 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR sqrt(chemv(k, p_chem_ic_cl_a01)*rij_g), sqrt(chemv(k, p_chem_ic_cl_a02)*rij_g), sqrt(chemv(k, p_chem_ic_cl_a03)*rij_g), & sqrt(chemv(k, p_chem_ic_na_a01)*rij_g), sqrt(chemv(k, p_chem_ic_na_a02)*rij_g), sqrt(chemv(k, p_chem_ic_na_a03)*rij_g), & sqrt(chemv(k, p_chem_ic_oin_a01)*rij_g), sqrt(chemv(k, p_chem_ic_oin_a02)*rij_g), sqrt(chemv(k, p_chem_ic_oin_a03)*rij_g) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,i4,39f8.3)') 'k=',k, & + (sqrt(chemv(k,ic)*rij_g), ic=PARAM_FIRST_SCALAR,num_chem) end select do ic=PARAM_FIRST_SCALAR, num_chem chemm(ic)=chemm(ic)+chemv(k,ic) @@ -817,6 +845,8 @@ use module_state_description, only : num_chem, PARAM_FIRST_SCALAR sqrt(chemm(p_chem_ic_cl_a01)*rijk_g), sqrt(chemm(p_chem_ic_cl_a02)*rijk_g), sqrt(chemm(p_chem_ic_cl_a03)*rijk_g), & sqrt(chemm(p_chem_ic_na_a01)*rijk_g), sqrt(chemm(p_chem_ic_na_a02)*rijk_g), sqrt(chemm(p_chem_ic_na_a03)*rijk_g), & sqrt(chemm(p_chem_ic_oin_a01)*rijk_g), sqrt(chemm(p_chem_ic_oin_a02)*rijk_g), sqrt(chemm(p_chem_ic_oin_a03)*rijk_g) + case ( 108 ) ! racm_soa_vbs_da + write(unit=stats_unit2,fmt='(A,39f8.3)') 'ALL', (sqrt(chemm(ic)*rijk_g),ic=PARAM_FIRST_SCALAR,num_chem) end select #endif end if @@ -828,5 +858,3 @@ contains #include "da_maxmin_in_field.inc" end subroutine da_analysis_stats - - diff --git a/var/da/da_statistics/da_print_qcstat.inc b/var/da/da_statistics/da_print_qcstat.inc index c8f93c4603..4ef0fddf03 100644 --- a/var/da/da_statistics/da_print_qcstat.inc +++ b/var/da/da_statistics/da_print_qcstat.inc @@ -12,7 +12,7 @@ subroutine da_print_qcstat(it, iv, num_qcstat_conv) integer :: ii, i, j, k ,n, num, ounit, ios character(len=filename_len) :: file - + character*4 :: ob_field ! Added by Soyoung (Aug-2021) logical :: write_head character(1),dimension(115):: cline @@ -81,12 +81,17 @@ subroutine da_print_qcstat(it, iv, num_qcstat_conv) write_head = .true. end if do j = 1, num_ob_vars + ob_field = ob_vars(j) ! Added by Soyoung (Aug-2021) +#if (WRF_CHEM == 1) + if (i.eq.num_ob_indexes) ob_field = 'chem' +#endif if( num_qcstat_conv(1,i,j,npres_print+1) > 0 ) then - write(ounit,700) obs_names(i),ob_vars(j),typx(1),& + write(ounit,700) obs_names(i),ob_field,typx(1),& ! Soyoung (Aug-2021) +! write(ounit,700) obs_names(i),ob_vars(j),typx(1),& ((num_qcstat_conv(1,i,j,k) - num_qcstat_conv(2,i,j,k)),k=1,npres_print+1) write(ounit,701) typx(2), (num_qcstat_conv(2,i,j,k),k=1,npres_print+1) -700 format(1x,a10,a3,2x,a4,3x, 25(1x,i6) ) -701 format(16x,a4,3x,25(1x,i6) ) +700 format(1x,a12,a4,2x,a4,1x, 25(1x,i6) ) +701 format(19x,a4,1x,25(1x,i6) ) end if end do end do diff --git a/var/da/da_transfer_model/da_transfer_model.f90 b/var/da/da_transfer_model/da_transfer_model.f90 index eac892c340..94ed01cbdf 100644 --- a/var/da/da_transfer_model/da_transfer_model.f90 +++ b/var/da/da_transfer_model/da_transfer_model.f90 @@ -29,7 +29,25 @@ module da_transfer_model p_chem_ic_bc_a01, p_chem_ic_bc_a02, p_chem_ic_bc_a03, p_chem_ic_oc_a01, p_chem_ic_oc_a02, p_chem_ic_oc_a03, & p_chem_ic_so4_a01, p_chem_ic_so4_a02, p_chem_ic_so4_a03, p_chem_ic_no3_a01, p_chem_ic_no3_a02, p_chem_ic_no3_a03, & p_chem_ic_nh4_a01, p_chem_ic_nh4_a02, p_chem_ic_nh4_a03, p_chem_ic_cl_a01, p_chem_ic_cl_a02, p_chem_ic_cl_a03, & - p_chem_ic_na_a01, p_chem_ic_na_a02, p_chem_ic_na_a03, p_chem_ic_oin_a01, p_chem_ic_oin_a02, p_chem_ic_oin_a03 + p_chem_ic_na_a01, p_chem_ic_na_a02, p_chem_ic_na_a03, p_chem_ic_oin_a01, p_chem_ic_oin_a02, p_chem_ic_oin_a03, & + p_so4aj, p_so4ai, p_nh4aj, p_nh4ai, & ! aerosols in racm_soa_vbs_kpp + p_no3aj, p_no3ai, p_naaj, p_naai, & + p_asoa1j, p_asoa1i, p_asoa2j, p_asoa2i, & + p_asoa3j, p_asoa3i, p_asoa4j, p_asoa4i, & + p_bsoa1j, p_bsoa1i, p_bsoa2j, p_bsoa2i, & + p_bsoa3j, p_bsoa3i, p_bsoa4j, p_bsoa4i, & + p_orgpaj, p_orgpai, p_ecj, p_eci, & + p_p25j, p_p25i, p_antha, p_seas, & + p_claj, p_clai, p_soila, & + p_chem_ic_so4aj, p_chem_ic_so4ai, p_chem_ic_nh4aj, p_chem_ic_nh4ai, & ! aerosols in racm_soa_vbs_da + p_chem_ic_no3aj, p_chem_ic_no3ai, p_chem_ic_naaj, p_chem_ic_naai, & + p_chem_ic_asoa1j, p_chem_ic_asoa1i, p_chem_ic_asoa2j, p_chem_ic_asoa2i, & + p_chem_ic_asoa3j, p_chem_ic_asoa3i, p_chem_ic_asoa4j, p_chem_ic_asoa4i, & + p_chem_ic_bsoa1j, p_chem_ic_bsoa1i, p_chem_ic_bsoa2j, p_chem_ic_bsoa2i, & + p_chem_ic_bsoa3j, p_chem_ic_bsoa3i, p_chem_ic_bsoa4j, p_chem_ic_bsoa4i, & + p_chem_ic_orgpaj, p_chem_ic_orgpai, p_chem_ic_ecj, p_chem_ic_eci, & + p_chem_ic_p25j, p_chem_ic_p25i, p_chem_ic_antha, p_chem_ic_seas, & + p_chem_ic_claj, p_chem_ic_clai, p_chem_ic_soila #endif use module_dm, only : wrf_dm_sum_real, wrf_dm_sum_reals diff --git a/var/da/da_transfer_model/da_transfer_xatoanalysis.inc b/var/da/da_transfer_model/da_transfer_xatoanalysis.inc index 557ee789f7..f7fb1ec171 100644 --- a/var/da/da_transfer_model/da_transfer_xatoanalysis.inc +++ b/var/da/da_transfer_model/da_transfer_xatoanalysis.inc @@ -4,14 +4,6 @@ subroutine da_transfer_xatoanalysis(it, xbx, grid, config_flags) ! Purpose: Transfer xb and xa (increments) to analysis. !--------------------------------------------------------------------------- -#if (WRF_CHEM == 1) - USE module_state_description, ONLY : & - num_chem & !inline gocart - ,p_chem_ic_bc1, p_chem_ic_bc2, p_chem_ic_oc1, p_chem_ic_oc2 & !inline gocart - ,p_chem_ic_so2, p_chem_ic_no2, p_chem_ic_o3, p_chem_ic_co & - ,p_chem_ic_p10, p_chem_ic_p25, p_chem_ic_sulf -#endif - implicit none integer, intent(in) :: it ! outer-loop index @@ -22,7 +14,6 @@ subroutine da_transfer_xatoanalysis(it, xbx, grid, config_flags) character*4 filnam character(len=256) :: timestr - integer n,ichem if (trace_use) call da_trace_entry("da_transfer_xatoanalysis") !--------------------------------------------------------------------------- diff --git a/var/da/da_vtox_transforms/da_transform_vchemtox_adj.inc b/var/da/da_vtox_transforms/da_transform_vchemtox_adj.inc index a7aacb46cf..5fa22950f6 100644 --- a/var/da/da_vtox_transforms/da_transform_vchemtox_adj.inc +++ b/var/da/da_vtox_transforms/da_transform_vchemtox_adj.inc @@ -26,8 +26,7 @@ subroutine da_transform_vchemtox_adj(grid, vchem, be) !!! vchem%chem_ic(:,grid%j_start(ij):grid%j_end(ij),:,ic) = 0.D0 do m = kts,kte !LEVELS do n = 1, be % v12(ic-1) % mz !EIGENMODES -! do j = grid%j_start(ij), grid%j_end(ij) - do j = jts, jte + do j = grid%j_start(ij), grid%j_end(ij) ! do j = jts, jte temp = be % v12(ic-1) % evec (j,m,n) * be % v12(ic-1) % val (j,n) do i = its, ite vchem%chem_ic (i,j,n,ic) = vchem%chem_ic (i,j,n,ic) + &