17 use camp_util,
only: string_t
19 use camp_chem_spec_data
20 use camp_aero_rep_data
21 use camp_aero_rep_single_particle
31 #ifdef PMC_USE_MOSAIC_MULTI_OPT
59 character(len=AERO_NAME_LEN),
allocatable :: name(:)
62 integer,
allocatable :: mosaic_index(:)
64 real(kind=dp),
allocatable :: wavelengths(:)
66 real(kind=dp),
allocatable :: density(:)
68 integer,
allocatable :: num_ions(:)
70 real(kind=dp),
allocatable :: molec_weight(:)
72 real(kind=dp),
allocatable :: kappa(:)
74 real(kind=dp),
allocatable :: abifm_m(:)
76 real(kind=dp),
allocatable :: abifm_c(:)
78 character(len=AERO_SOURCE_NAME_LEN),
allocatable :: source_name(:)
80 character(len=AERO_SOURCE_NAME_LEN),
allocatable :: weight_class_name(:)
85 class(aero_rep_data_t),
pointer :: aero_rep_ptr
88 integer,
allocatable :: camp_particle_spec_id(:)
91 integer :: camp_particle_state_size = -1
110 real(kind=dp),
intent(in) :: v
125 real(kind=dp),
intent(in) :: v
140 real(kind=dp),
intent(in) :: r
155 real(kind=dp),
intent(in) :: d
173 real(kind=dp),
intent(in) :: v
176 aero_data%fractal, v)
192 real(kind=dp),
intent(in) :: v
194 real(kind=dp),
intent(in) :: temp
196 real(kind=dp),
intent(in) :: pressure
199 aero_data%fractal, v, temp, pressure)
208 mobility_rad, temp, pressure)
213 real(kind=dp),
intent(in) :: mobility_rad
215 real(kind=dp),
intent(in) :: temp
217 real(kind=dp),
intent(in) :: pressure
220 aero_data%fractal, mobility_rad, temp, pressure)
229 mobility_rad, temp, pressure)
234 real(kind=dp),
intent(in) :: mobility_rad
236 real(kind=dp),
intent(in) :: temp
238 real(kind=dp),
intent(in) :: pressure
242 mobility_rad, temp, pressure)
254 if (
allocated(aero_data%name))
then
270 if (
allocated(aero_data%source_name))
then
286 if (
allocated(aero_data%weight_class_name))
then
303 character(len=*),
intent(in) :: name
310 if (name == aero_data%name(i))
then
332 character(len=*),
intent(in) :: name
334 if (.not.
allocated(aero_data%source_name))
then
341 aero_data%source_name = [aero_data%source_name, &
356 character(len=*),
intent(in) :: name
358 if (.not.
allocated(aero_data%weight_class_name))
then
364 aero_data%weight_class_name, name)
366 aero_data%weight_class_name = [aero_data%weight_class_name, &
394 integer,
parameter :: n_mosaic_spec = 19
395 character(AERO_NAME_LEN),
parameter,
dimension(n_mosaic_spec) :: &
396 mosaic_spec_name = [ &
397 "SO4 ",
"NO3 ",
"Cl ",
"NH4 ",
"MSA ",
"ARO1 ", &
398 "ARO2 ",
"ALK1 ",
"OLE1 ",
"API1 ",
"API2 ",
"LIM1 ", &
399 "LIM2 ",
"CO3 ",
"Na ",
"Ca ",
"OIN ",
"OC ", &
402 integer :: i_spec, i_mosaic_spec, i
404 aero_data%mosaic_index = 0
407 aero_data%name(i_spec))
421 integer,
intent(in) :: i_part
423 integer,
intent(in) :: i_spec
426 call assert(106669451,
allocated(aero_data%camp_particle_spec_id))
427 call assert(278731889, aero_data%camp_particle_state_size .ge. 0)
428 camp_spec_id = (i_part - 1) * aero_data%camp_particle_state_size + &
429 aero_data%camp_particle_spec_id(i_spec)
439 subroutine spec_file_read_aero_data(file, aero_data)
446 integer :: n_species, species, i
447 character(len=SPEC_LINE_MAX_VAR_LEN),
allocatable :: species_name(:)
448 real(kind=dp),
allocatable :: species_data(:,:)
489 n_species =
size(species_data, 1)
490 if (.not. ((
size(species_data, 2) == 6) .or. (n_species == 0)))
then
491 call die_msg(428926381,
'each line in ' // trim(file%name) &
492 //
' should contain exactly 7 values')
507 aero_data%density(i) = species_data(i,1)
508 aero_data%num_ions(i) = nint(species_data(i,2))
509 aero_data%molec_weight(i) = species_data(i,3)
510 aero_data%kappa(i) = species_data(i,4)
511 aero_data%abifm_m(i) = species_data(i,5)
512 aero_data%abifm_c(i) = species_data(i,6)
514 (aero_data%num_ions(i) == 0) .or. (aero_data%kappa(i) == 0d0), &
515 "ions and kappa both non-zero for species " &
516 // trim(aero_data%name(i)) //
" in " // trim(file%name))
517 if (species_name(i) ==
"H2O")
then
518 aero_data%i_water = i
520 aero_data%density(i), const%water_density), &
521 "input H2O density not equal to const%water_density (" &
525 aero_data%molec_weight(i),const%water_molec_weight), &
526 "input H2O molec_weight not equal " &
527 //
"to const%water_molec_weight (" &
535 end subroutine spec_file_read_aero_data
545 character(len=*),
intent(in) :: name
549 integer,
allocatable :: species_list(:)
551 type(spec_line_t) :: line
557 do i = 1,
size(line%data)
561 'unknown species: ' // trim(line%data(i)))
563 species_list(i) = spec
599 character,
intent(inout) :: buffer(:)
601 integer,
intent(inout) :: position
606 integer :: prev_position
608 prev_position = position
634 character,
intent(inout) :: buffer(:)
636 integer,
intent(inout) :: position
641 integer :: prev_position
643 prev_position = position
674 integer,
intent(in) :: ncid
676 integer,
intent(out) :: dimid_aero_species
678 integer :: status, i_spec
679 integer :: varid_aero_species
680 integer :: aero_species_centers(aero_data_n_spec(aero_data))
681 character(len=(AERO_NAME_LEN * aero_data_n_spec(aero_data))) :: &
685 status = nf90_inq_dimid(ncid,
"aero_species", dimid_aero_species)
686 if (status == nf90_noerr)
return
693 aero_data_n_spec(aero_data), dimid_aero_species))
694 aero_species_names =
""
695 do i_spec = 1,aero_data_n_spec(aero_data)
696 aero_species_names((len_trim(aero_species_names) + 1):) &
697 = trim(aero_data%name(i_spec))
698 if (i_spec < aero_data_n_spec(aero_data))
then
699 aero_species_names((len_trim(aero_species_names) + 1):) =
","
702 call pmc_nc_check(nf90_def_var(ncid,
"aero_species", nf90_int, &
703 dimid_aero_species, varid_aero_species))
704 call pmc_nc_check(nf90_put_att(ncid, varid_aero_species,
"names", &
706 call pmc_nc_check(nf90_put_att(ncid, varid_aero_species,
"description", &
707 "dummy dimension variable (no useful value) - read species names " &
708 //
"as comma-separated values from the 'names' attribute"))
712 do i_spec = 1,aero_data_n_spec(aero_data)
713 aero_species_centers(i_spec) = i_spec
715 call pmc_nc_check(nf90_put_var(ncid, varid_aero_species, &
716 aero_species_centers))
731 integer,
intent(in) :: ncid
733 integer,
intent(out) :: dimid_aero_source
735 integer :: status, i_source
736 integer :: varid_aero_source
737 integer :: aero_source_centers(aero_data_n_source(aero_data))
738 character(len=(AERO_SOURCE_NAME_LEN * aero_data_n_source(aero_data))) &
742 status = nf90_inq_dimid(ncid,
"aero_source", dimid_aero_source)
743 if (status == nf90_noerr)
return
750 aero_data_n_source(aero_data), dimid_aero_source))
751 aero_source_names =
""
752 do i_source = 1,aero_data_n_source(aero_data)
753 aero_source_names((len_trim(aero_source_names) + 1):) &
754 = trim(aero_data%source_name(i_source))
755 if (i_source < aero_data_n_source(aero_data))
then
756 aero_source_names((len_trim(aero_source_names) + 1):) =
","
759 call pmc_nc_check(nf90_def_var(ncid,
"aero_source", nf90_int, &
760 dimid_aero_source, varid_aero_source))
761 call pmc_nc_check(nf90_put_att(ncid, varid_aero_source,
"names", &
763 call pmc_nc_check(nf90_put_att(ncid, varid_aero_source,
"description", &
764 "dummy dimension variable (no useful value) - read source names " &
765 //
"as comma-separated values from the 'names' attribute"))
769 do i_source = 1,aero_data_n_source(aero_data)
770 aero_source_centers(i_source) = i_source
772 call pmc_nc_check(nf90_put_var(ncid, varid_aero_source, &
773 aero_source_centers))
783 dimid_aero_weight_class)
788 integer,
intent(in) :: ncid
790 integer,
intent(out) :: dimid_aero_weight_class
792 integer :: status, i_source, i_class
793 integer :: varid_aero_weight_class
794 integer :: aero_weight_class_centers(aero_data_n_weight_class(aero_data))
795 character(len=(AERO_SOURCE_NAME_LEN &
* aero_data_n_weight_class(aero_data))) :: aero_weight_class_names
798 status = nf90_inq_dimid(ncid,
"aero_weight_class", dimid_aero_weight_class)
799 if (status == nf90_noerr)
return
805 call pmc_nc_check(nf90_def_dim(ncid,
"aero_weight_class", &
806 aero_data_n_weight_class(aero_data), dimid_aero_weight_class))
807 aero_weight_class_names =
""
808 do i_class = 1,aero_data_n_weight_class(aero_data)
809 aero_weight_class_names((len_trim(aero_weight_class_names) + 1):) &
810 = trim(aero_data%weight_class_name(i_class))
811 if (i_class < aero_data_n_weight_class(aero_data))
then
812 aero_weight_class_names((len_trim(aero_weight_class_names) + 1):) &
816 call pmc_nc_check(nf90_def_var(ncid,
"aero_data_weight_class", nf90_int, &
817 dimid_aero_weight_class, varid_aero_weight_class))
818 call pmc_nc_check(nf90_put_att(ncid, varid_aero_weight_class,
"names", &
819 aero_weight_class_names))
820 call pmc_nc_check(nf90_put_att(ncid, varid_aero_weight_class, &
822 "dummy dimension variable (no useful value) - read source names " &
823 //
"as comma-separated values from the 'names' attribute"))
827 do i_class = 1,aero_data_n_weight_class(aero_data)
828 aero_weight_class_centers(i_class) = i_class
830 call pmc_nc_check(nf90_put_var(ncid, varid_aero_weight_class, &
831 aero_weight_class_centers))
841 dimid_optical_wavelengths)
846 integer,
intent(in) :: ncid
848 integer,
intent(out) :: dimid_optical_wavelengths
853 status = nf90_inq_dimid(ncid,
"optical_wavelengths", &
854 dimid_optical_wavelengths)
855 if (status == nf90_noerr)
return
861 call pmc_nc_check(nf90_def_dim(ncid,
"optical_wavelengths", &
871 subroutine aero_data_output_netcdf(aero_data, ncid)
876 integer,
intent(in) :: ncid
878 integer :: dimid_aero_species, dimid_aero_source, &
879 dimid_aero_weight_classes, dimid_optical_wavelengths
916 dimid_aero_weight_classes)
918 dimid_optical_wavelengths)
921 "aero_mosaic_index", (/ dimid_aero_species /), &
922 long_name=
"MOSAIC indices of aerosol species")
924 "aero_optical_wavelengths", (/ dimid_optical_wavelengths /), &
925 long_name=
"wavelength of optical calculations", unit=
"m")
927 "aero_density", (/ dimid_aero_species /), unit=
"kg/m^3", &
928 long_name=
"densities of aerosol species")
930 "aero_num_ions", (/ dimid_aero_species /), &
931 long_name=
"number of ions after dissociation of aerosol species")
933 "aero_molec_weight", (/ dimid_aero_species /), unit=
"kg/mol", &
934 long_name=
"molecular weights of aerosol species")
936 "aero_kappa", (/ dimid_aero_species /), unit=
"1", &
937 long_name=
"hygroscopicity parameters (kappas) of aerosol species")
939 "aero_abifm_m", (/ dimid_aero_species /), unit=
"1", &
940 long_name=
"m parameter of ABIFM")
942 "aero_abifm_c", (/ dimid_aero_species /), unit=
"1", &
943 long_name=
"c parameter of ABIFM")
945 "aero_i_water", long_name=
"Index of aerosol water or " &
946 //
"0 if water does not exist.")
947 call fractal_output_netcdf(aero_data%fractal, ncid)
949 end subroutine aero_data_output_netcdf
959 integer,
intent(in) :: ncid
961 integer,
parameter :: MAX_SPECIES = 1000
962 integer,
parameter :: MAX_SOURCES = 1000
964 character(len=1000) :: name
965 integer :: dimid_aero_species, n_spec, varid_aero_species, i_spec, i
966 integer :: dimid_aero_source, n_source, varid_aero_source, i_source
967 integer :: dimid_aero_weight_class, n_class, i_class, &
968 varid_aero_weight_class
969 character(len=((AERO_NAME_LEN + 2) * MAX_SPECIES)) :: aero_species_names
970 character(len=:),
allocatable :: aero_source_names
971 character(len=:),
allocatable :: aero_weight_class_names
973 call pmc_nc_check(nf90_inq_dimid(ncid,
"aero_species", &
976 dimid_aero_species, name, n_spec))
977 call assert(141013948, n_spec < max_species)
982 dimid_aero_source, name, n_source))
983 call assert(739238793, n_source < max_sources)
988 "aero_optical_wavelengths")
996 call pmc_nc_check(nf90_inq_varid(ncid,
"aero_species", &
998 call pmc_nc_check(nf90_get_att(ncid, varid_aero_species,
"names", &
1004 do while ((aero_species_names(i:i) /=
" ") &
1005 .and. (aero_species_names(i:i) /=
","))
1008 call assert(852937292, i > 1)
1009 aero_data%name(i_spec) = aero_species_names(1:(i-1))
1010 aero_species_names = aero_species_names((i+1):)
1012 call assert(729138192, aero_species_names ==
"")
1014 call pmc_nc_check(nf90_inq_varid(ncid,
"aero_source", &
1016 allocate(
character(len=((AERO_SOURCE_NAME_LEN + 2) * MAX_SOURCES)) &
1017 :: aero_source_names)
1018 call pmc_nc_check(nf90_get_att(ncid, varid_aero_source,
"names", &
1024 do while ((aero_source_names(i:i) /=
" ") &
1025 .and. (aero_source_names(i:i) /=
","))
1028 call assert(840982478, i > 1)
1029 aero_data%source_name(i_source) = aero_source_names(1:(i-1))
1030 aero_source_names = aero_source_names((i+1):)
1032 call assert(377166446, aero_source_names ==
"")
1034 call pmc_nc_check(nf90_inq_varid(ncid,
"aero_data_weight_class", &
1035 varid_aero_weight_class))
1036 allocate(
character(len=((AERO_SOURCE_NAME_LEN + 2) * MAX_SOURCES)) &
1037 :: aero_weight_class_names)
1038 call pmc_nc_check(nf90_get_att(ncid, varid_aero_weight_class,
"names", &
1039 aero_weight_class_names))
1041 call pmc_nc_check(nf90_inq_dimid(ncid,
"aero_weight_class", &
1042 dimid_aero_weight_class))
1044 dimid_aero_weight_class, name, n_class))
1046 do i_class = 1,n_class
1048 do while ((aero_weight_class_names(i:i) /=
" ") &
1049 .and. (aero_weight_class_names(i:i) /=
","))
1052 call assert(840982472, i > 1)
1053 aero_data%weight_class_name(i_class) = aero_weight_class_names(1:(i-1))
1054 aero_weight_class_names = aero_weight_class_names((i+1):)
1056 call assert(377166448, aero_weight_class_names ==
"")
1069 subroutine aero_data_initialize(aero_data, camp_core)
1074 type(camp_core_t),
intent(in) :: camp_core
1076 character(len=:),
allocatable :: rep_name, prop_name, str_val
1077 type(string_t),
allocatable :: spec_names(:), tmp_spec_names(:)
1078 integer :: num_spec, i_spec, spec_type
1079 type(chem_spec_data_t),
pointer :: chem_spec_data
1080 type(property_t),
pointer :: property_set
1082 rep_name =
"PartMC single particle"
1083 if (.not. camp_core%get_aero_rep(rep_name, aero_data%aero_rep_ptr))
then
1084 call die_msg(418509983,
"Missing 'PartMC single particle' aerosol "// &
1088 call assert_msg(935419266, camp_core%get_chem_spec_data(chem_spec_data), &
1089 "No chemical species data in camp_core.")
1092 spec_names = aero_data%aero_rep_ptr%unique_names()
1093 allocate(tmp_spec_names(
size(spec_names)))
1095 do i_spec = 1,
size(spec_names)
1096 call assert(496388827, chem_spec_data%get_type( &
1097 aero_data%aero_rep_ptr%spec_name(spec_names(i_spec)%string), &
1099 if (spec_type /= chem_spec_variable .and. &
1100 spec_type /= chem_spec_constant .and. &
1101 spec_type /= chem_spec_pssa) cycle
1102 if (spec_names(i_spec)%string(1:3) /=
"P1.")
exit
1103 num_spec = num_spec + 1
1104 tmp_spec_names(num_spec)%string = spec_names(i_spec)%string(4:)
1106 deallocate(spec_names)
1107 allocate(spec_names(num_spec))
1108 spec_names(:) = tmp_spec_names(1:num_spec)
1109 deallocate(tmp_spec_names)
1111 allocate(aero_data%name(num_spec))
1112 allocate(aero_data%mosaic_index(num_spec))
1113 allocate(aero_data%wavelengths(
n_swbands))
1114 allocate(aero_data%density(num_spec))
1115 allocate(aero_data%num_ions(num_spec))
1116 allocate(aero_data%molec_weight(num_spec))
1117 allocate(aero_data%kappa(num_spec))
1118 allocate(aero_data%abifm_m(num_spec))
1119 allocate(aero_data%abifm_c(num_spec))
1120 allocate(aero_data%camp_particle_spec_id(num_spec))
1123 aero_data%i_water = 0
1125 do i_spec = 1, num_spec
1126 aero_data%name(i_spec) = spec_names(i_spec)%string
1127 if (.not. chem_spec_data%get_property_set( &
1128 aero_data%aero_rep_ptr%spec_name(
"P1." &
1129 // spec_names(i_spec)%string), property_set))
then
1130 call die_msg(934844845,
"Missing property set for aerosol species " &
1131 // spec_names(i_spec)%string)
1133 prop_name =
"density [kg m-3]"
1134 if (.not. property_set%get_real(prop_name, &
1135 aero_data%density(i_spec)))
then
1136 call die_msg(547508215,
"Missing density for aerosol species " &
1137 // spec_names(i_spec)%string)
1139 prop_name =
"num_ions"
1140 if (.not. property_set%get_int(prop_name, &
1141 aero_data%num_ions(i_spec)))
then
1142 call die_msg(324777059,
"Missing num_ions for aerosol species " &
1143 // spec_names(i_spec)%string)
1145 prop_name =
"molecular weight [kg mol-1]"
1146 if (.not. property_set%get_real(prop_name, &
1147 aero_data%molec_weight(i_spec)))
then
1148 call die_msg(549413749,
"Missing molec_weight for aerosol species " &
1149 // spec_names(i_spec)%string)
1152 if (.not. property_set%get_real(prop_name, &
1153 aero_data%kappa(i_spec)))
then
1154 call die_msg(944207343,
"Missing kappa for aerosol species " &
1155 // spec_names(i_spec)%string)
1157 prop_name =
"abifm_m"
1158 if (.not. property_set%get_real(prop_name, &
1159 aero_data%abifm_m(i_spec)))
then
1160 call die_msg(944207345,
"Missing abifm_m for aerosol species " &
1161 // spec_names(i_spec)%string)
1163 prop_name =
"abifm_c"
1164 if (.not. property_set%get_real(prop_name, &
1165 aero_data%abifm_c(i_spec)))
then
1166 call die_msg(944207346,
"Missing abifm_c for aerosol species " &
1167 // spec_names(i_spec)%string)
1169 prop_name =
"PartMC name"
1170 if (property_set%get_string(prop_name, str_val))
then
1171 if (str_val ==
"H2O")
then
1172 call assert_msg(227489086, aero_data%i_water == 0, &
1173 "Multiple aerosol water species")
1174 aero_data%i_water = i_spec
1177 aero_data%camp_particle_spec_id(i_spec) = &
1178 aero_data%aero_rep_ptr%spec_state_id(
"P1." &
1179 // spec_names(i_spec)%string)
1182 select type( aero_rep => aero_data%aero_rep_ptr)
1183 type is(aero_rep_single_particle_t)
1186 aero_data%camp_particle_state_size = aero_rep%per_particle_size()
1189 call die_msg(281737350,
"Wrong aerosol representation type")
1192 end subroutine aero_data_initialize