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 character(len=AERO_SOURCE_NAME_LEN),
allocatable :: source_name(:)
76 character(len=AERO_SOURCE_NAME_LEN),
allocatable :: weight_class_name(:)
81 class(aero_rep_data_t),
pointer :: aero_rep_ptr
84 integer,
allocatable :: camp_particle_spec_id(:)
87 integer :: camp_particle_state_size = -1
106 real(kind=dp),
intent(in) :: v
121 real(kind=dp),
intent(in) :: v
136 real(kind=dp),
intent(in) :: r
151 real(kind=dp),
intent(in) :: d
169 real(kind=dp),
intent(in) :: v
172 aero_data%fractal, v)
188 real(kind=dp),
intent(in) :: v
190 real(kind=dp),
intent(in) :: temp
192 real(kind=dp),
intent(in) :: pressure
195 aero_data%fractal, v, temp, pressure)
204 mobility_rad, temp, pressure)
209 real(kind=dp),
intent(in) :: mobility_rad
211 real(kind=dp),
intent(in) :: temp
213 real(kind=dp),
intent(in) :: pressure
216 aero_data%fractal, mobility_rad, temp, pressure)
225 mobility_rad, temp, pressure)
230 real(kind=dp),
intent(in) :: mobility_rad
232 real(kind=dp),
intent(in) :: temp
234 real(kind=dp),
intent(in) :: pressure
238 mobility_rad, temp, pressure)
250 if (
allocated(aero_data%name))
then
266 if (
allocated(aero_data%source_name))
then
282 if (
allocated(aero_data%weight_class_name))
then
299 character(len=*),
intent(in) :: name
306 if (name == aero_data%name(i))
then
328 character(len=*),
intent(in) :: name
330 if (.not.
allocated(aero_data%source_name))
then
337 aero_data%source_name = [aero_data%source_name, &
352 character(len=*),
intent(in) :: name
354 if (.not.
allocated(aero_data%weight_class_name))
then
360 aero_data%weight_class_name, name)
362 aero_data%weight_class_name = [aero_data%weight_class_name, &
390 integer,
parameter :: n_mosaic_spec = 19
391 character(AERO_NAME_LEN),
parameter,
dimension(n_mosaic_spec) :: &
392 mosaic_spec_name = [ &
393 "SO4 ",
"NO3 ",
"Cl ",
"NH4 ",
"MSA ",
"ARO1 ", &
394 "ARO2 ",
"ALK1 ",
"OLE1 ",
"API1 ",
"API2 ",
"LIM1 ", &
395 "LIM2 ",
"CO3 ",
"Na ",
"Ca ",
"OIN ",
"OC ", &
398 integer :: i_spec, i_mosaic_spec, i
400 aero_data%mosaic_index = 0
403 aero_data%name(i_spec))
417 integer,
intent(in) :: i_part
419 integer,
intent(in) :: i_spec
422 call assert(106669451,
allocated(aero_data%camp_particle_spec_id))
423 call assert(278731889, aero_data%camp_particle_state_size .ge. 0)
424 camp_spec_id = (i_part - 1) * aero_data%camp_particle_state_size + &
425 aero_data%camp_particle_spec_id(i_spec)
435 subroutine spec_file_read_aero_data(file, aero_data)
442 integer :: n_species, species, i
443 character(len=SPEC_LINE_MAX_VAR_LEN),
allocatable :: species_name(:)
444 real(kind=dp),
allocatable :: species_data(:,:)
483 n_species =
size(species_data, 1)
484 if (.not. ((
size(species_data, 2) == 4) .or. (n_species == 0)))
then
485 call die_msg(428926381,
'each line in ' // trim(file%name) &
486 //
' should contain exactly 5 values')
499 aero_data%density(i) = species_data(i,1)
500 aero_data%num_ions(i) = nint(species_data(i,2))
501 aero_data%molec_weight(i) = species_data(i,3)
502 aero_data%kappa(i) = species_data(i,4)
504 (aero_data%num_ions(i) == 0) .or. (aero_data%kappa(i) == 0d0), &
505 "ions and kappa both non-zero for species " &
506 // trim(aero_data%name(i)) //
" in " // trim(file%name))
507 if (species_name(i) ==
"H2O")
then
508 aero_data%i_water = i
510 aero_data%density(i), const%water_density), &
511 "input H2O density not equal to const%water_density (" &
515 aero_data%molec_weight(i),const%water_molec_weight), &
516 "input H2O molec_weight not equal " &
517 //
"to const%water_molec_weight (" &
525 end subroutine spec_file_read_aero_data
535 character(len=*),
intent(in) :: name
539 integer,
allocatable :: species_list(:)
541 type(spec_line_t) :: line
547 do i = 1,
size(line%data)
551 'unknown species: ' // trim(line%data(i)))
553 species_list(i) = spec
587 character,
intent(inout) :: buffer(:)
589 integer,
intent(inout) :: position
594 integer :: prev_position
596 prev_position = position
620 character,
intent(inout) :: buffer(:)
622 integer,
intent(inout) :: position
627 integer :: prev_position
629 prev_position = position
658 integer,
intent(in) :: ncid
660 integer,
intent(out) :: dimid_aero_species
662 integer :: status, i_spec
663 integer :: varid_aero_species
664 integer :: aero_species_centers(aero_data_n_spec(aero_data))
665 character(len=(AERO_NAME_LEN * aero_data_n_spec(aero_data))) :: &
669 status = nf90_inq_dimid(ncid,
"aero_species", dimid_aero_species)
670 if (status == nf90_noerr)
return
677 aero_data_n_spec(aero_data), dimid_aero_species))
678 aero_species_names =
""
679 do i_spec = 1,aero_data_n_spec(aero_data)
680 aero_species_names((len_trim(aero_species_names) + 1):) &
681 = trim(aero_data%name(i_spec))
682 if (i_spec < aero_data_n_spec(aero_data))
then
683 aero_species_names((len_trim(aero_species_names) + 1):) =
","
686 call pmc_nc_check(nf90_def_var(ncid,
"aero_species", nf90_int, &
687 dimid_aero_species, varid_aero_species))
688 call pmc_nc_check(nf90_put_att(ncid, varid_aero_species,
"names", &
690 call pmc_nc_check(nf90_put_att(ncid, varid_aero_species,
"description", &
691 "dummy dimension variable (no useful value) - read species names " &
692 //
"as comma-separated values from the 'names' attribute"))
696 do i_spec = 1,aero_data_n_spec(aero_data)
697 aero_species_centers(i_spec) = i_spec
699 call pmc_nc_check(nf90_put_var(ncid, varid_aero_species, &
700 aero_species_centers))
715 integer,
intent(in) :: ncid
717 integer,
intent(out) :: dimid_aero_source
719 integer :: status, i_source
720 integer :: varid_aero_source
721 integer :: aero_source_centers(aero_data_n_source(aero_data))
722 character(len=(AERO_SOURCE_NAME_LEN * aero_data_n_source(aero_data))) &
726 status = nf90_inq_dimid(ncid,
"aero_source", dimid_aero_source)
727 if (status == nf90_noerr)
return
734 aero_data_n_source(aero_data), dimid_aero_source))
735 aero_source_names =
""
736 do i_source = 1,aero_data_n_source(aero_data)
737 aero_source_names((len_trim(aero_source_names) + 1):) &
738 = trim(aero_data%source_name(i_source))
739 if (i_source < aero_data_n_source(aero_data))
then
740 aero_source_names((len_trim(aero_source_names) + 1):) =
","
743 call pmc_nc_check(nf90_def_var(ncid,
"aero_source", nf90_int, &
744 dimid_aero_source, varid_aero_source))
745 call pmc_nc_check(nf90_put_att(ncid, varid_aero_source,
"names", &
747 call pmc_nc_check(nf90_put_att(ncid, varid_aero_source,
"description", &
748 "dummy dimension variable (no useful value) - read source names " &
749 //
"as comma-separated values from the 'names' attribute"))
753 do i_source = 1,aero_data_n_source(aero_data)
754 aero_source_centers(i_source) = i_source
756 call pmc_nc_check(nf90_put_var(ncid, varid_aero_source, &
757 aero_source_centers))
767 dimid_aero_weight_class)
772 integer,
intent(in) :: ncid
774 integer,
intent(out) :: dimid_aero_weight_class
776 integer :: status, i_source, i_class
777 integer :: varid_aero_weight_class
778 integer :: aero_weight_class_centers(aero_data_n_weight_class(aero_data))
779 character(len=(AERO_SOURCE_NAME_LEN &
* aero_data_n_weight_class(aero_data))) :: aero_weight_class_names
782 status = nf90_inq_dimid(ncid,
"aero_weight_class", dimid_aero_weight_class)
783 if (status == nf90_noerr)
return
789 call pmc_nc_check(nf90_def_dim(ncid,
"aero_weight_class", &
790 aero_data_n_weight_class(aero_data), dimid_aero_weight_class))
791 aero_weight_class_names =
""
792 do i_class = 1,aero_data_n_weight_class(aero_data)
793 aero_weight_class_names((len_trim(aero_weight_class_names) + 1):) &
794 = trim(aero_data%weight_class_name(i_class))
795 if (i_class < aero_data_n_weight_class(aero_data))
then
796 aero_weight_class_names((len_trim(aero_weight_class_names) + 1):) &
800 call pmc_nc_check(nf90_def_var(ncid,
"aero_data_weight_class", nf90_int, &
801 dimid_aero_weight_class, varid_aero_weight_class))
802 call pmc_nc_check(nf90_put_att(ncid, varid_aero_weight_class,
"names", &
803 aero_weight_class_names))
804 call pmc_nc_check(nf90_put_att(ncid, varid_aero_weight_class, &
806 "dummy dimension variable (no useful value) - read source names " &
807 //
"as comma-separated values from the 'names' attribute"))
811 do i_class = 1,aero_data_n_weight_class(aero_data)
812 aero_weight_class_centers(i_class) = i_class
814 call pmc_nc_check(nf90_put_var(ncid, varid_aero_weight_class, &
815 aero_weight_class_centers))
825 dimid_optical_wavelengths)
830 integer,
intent(in) :: ncid
832 integer,
intent(out) :: dimid_optical_wavelengths
837 status = nf90_inq_dimid(ncid,
"optical_wavelengths", &
838 dimid_optical_wavelengths)
839 if (status == nf90_noerr)
return
845 call pmc_nc_check(nf90_def_dim(ncid,
"optical_wavelengths", &
855 subroutine aero_data_output_netcdf(aero_data, ncid)
860 integer,
intent(in) :: ncid
862 integer :: dimid_aero_species, dimid_aero_source, &
863 dimid_aero_weight_classes, dimid_optical_wavelengths
898 dimid_aero_weight_classes)
900 dimid_optical_wavelengths)
903 "aero_mosaic_index", (/ dimid_aero_species /), &
904 long_name=
"MOSAIC indices of aerosol species")
906 "aero_optical_wavelengths", (/ dimid_optical_wavelengths /), &
907 long_name=
"wavelength of optical calculations", unit=
"m")
909 "aero_density", (/ dimid_aero_species /), unit=
"kg/m^3", &
910 long_name=
"densities of aerosol species")
912 "aero_num_ions", (/ dimid_aero_species /), &
913 long_name=
"number of ions after dissociation of aerosol species")
915 "aero_molec_weight", (/ dimid_aero_species /), unit=
"kg/mol", &
916 long_name=
"molecular weights of aerosol species")
918 "aero_kappa", (/ dimid_aero_species /), unit=
"1", &
919 long_name=
"hygroscopicity parameters (kappas) of aerosol species")
921 "aero_i_water", long_name=
"Index of aerosol water or " &
922 //
"0 if water does not exist.")
923 call fractal_output_netcdf(aero_data%fractal, ncid)
925 end subroutine aero_data_output_netcdf
935 integer,
intent(in) :: ncid
937 integer,
parameter :: MAX_SPECIES = 1000
938 integer,
parameter :: MAX_SOURCES = 1000
940 character(len=1000) :: name
941 integer :: dimid_aero_species, n_spec, varid_aero_species, i_spec, i
942 integer :: dimid_aero_source, n_source, varid_aero_source, i_source
943 integer :: dimid_aero_weight_class, n_class, i_class, &
944 varid_aero_weight_class
945 character(len=((AERO_NAME_LEN + 2) * MAX_SPECIES)) :: aero_species_names
946 character(len=:),
allocatable :: aero_source_names
947 character(len=:),
allocatable :: aero_weight_class_names
949 call pmc_nc_check(nf90_inq_dimid(ncid,
"aero_species", &
952 dimid_aero_species, name, n_spec))
953 call assert(141013948, n_spec < max_species)
958 dimid_aero_source, name, n_source))
959 call assert(739238793, n_source < max_sources)
964 "aero_optical_wavelengths")
970 call pmc_nc_check(nf90_inq_varid(ncid,
"aero_species", &
972 call pmc_nc_check(nf90_get_att(ncid, varid_aero_species,
"names", &
978 do while ((aero_species_names(i:i) /=
" ") &
979 .and. (aero_species_names(i:i) /=
","))
982 call assert(852937292, i > 1)
983 aero_data%name(i_spec) = aero_species_names(1:(i-1))
984 aero_species_names = aero_species_names((i+1):)
986 call assert(729138192, aero_species_names ==
"")
990 allocate(
character(len=((AERO_SOURCE_NAME_LEN + 2) * MAX_SOURCES)) &
991 :: aero_source_names)
992 call pmc_nc_check(nf90_get_att(ncid, varid_aero_source,
"names", &
998 do while ((aero_source_names(i:i) /=
" ") &
999 .and. (aero_source_names(i:i) /=
","))
1002 call assert(840982478, i > 1)
1003 aero_data%source_name(i_source) = aero_source_names(1:(i-1))
1004 aero_source_names = aero_source_names((i+1):)
1006 call assert(377166446, aero_source_names ==
"")
1008 call pmc_nc_check(nf90_inq_varid(ncid,
"aero_data_weight_class", &
1009 varid_aero_weight_class))
1010 allocate(
character(len=((AERO_SOURCE_NAME_LEN + 2) * MAX_SOURCES)) &
1011 :: aero_weight_class_names)
1012 call pmc_nc_check(nf90_get_att(ncid, varid_aero_weight_class,
"names", &
1013 aero_weight_class_names))
1015 call pmc_nc_check(nf90_inq_dimid(ncid,
"aero_weight_class", &
1016 dimid_aero_weight_class))
1018 dimid_aero_weight_class, name, n_class))
1020 do i_class = 1,n_class
1022 do while ((aero_weight_class_names(i:i) /=
" ") &
1023 .and. (aero_weight_class_names(i:i) /=
","))
1026 call assert(840982472, i > 1)
1027 aero_data%weight_class_name(i_class) = aero_weight_class_names(1:(i-1))
1028 aero_weight_class_names = aero_weight_class_names((i+1):)
1030 call assert(377166448, aero_weight_class_names ==
"")
1043 subroutine aero_data_initialize(aero_data, camp_core)
1048 type(camp_core_t),
intent(in) :: camp_core
1050 character(len=:),
allocatable :: rep_name, prop_name, str_val
1051 type(string_t),
allocatable :: spec_names(:), tmp_spec_names(:)
1052 integer :: num_spec, i_spec, spec_type
1053 type(chem_spec_data_t),
pointer :: chem_spec_data
1054 type(property_t),
pointer :: property_set
1056 rep_name =
"PartMC single particle"
1057 if (.not. camp_core%get_aero_rep(rep_name, aero_data%aero_rep_ptr))
then
1058 call die_msg(418509983,
"Missing 'PartMC single particle' aerosol "// &
1062 call assert_msg(935419266, camp_core%get_chem_spec_data(chem_spec_data), &
1063 "No chemical species data in camp_core.")
1066 spec_names = aero_data%aero_rep_ptr%unique_names()
1067 allocate(tmp_spec_names(
size(spec_names)))
1069 do i_spec = 1,
size(spec_names)
1070 call assert(496388827, chem_spec_data%get_type( &
1071 aero_data%aero_rep_ptr%spec_name(spec_names(i_spec)%string), &
1073 if (spec_type /= chem_spec_variable .and. &
1074 spec_type /= chem_spec_constant .and. &
1075 spec_type /= chem_spec_pssa) cycle
1076 if (spec_names(i_spec)%string(1:3) /=
"P1.")
exit
1077 num_spec = num_spec + 1
1078 tmp_spec_names(num_spec)%string = spec_names(i_spec)%string(4:)
1080 deallocate(spec_names)
1081 allocate(spec_names(num_spec))
1082 spec_names(:) = tmp_spec_names(1:num_spec)
1083 deallocate(tmp_spec_names)
1085 allocate(aero_data%name(num_spec))
1086 allocate(aero_data%mosaic_index(num_spec))
1087 allocate(aero_data%wavelengths(
n_swbands))
1088 allocate(aero_data%density(num_spec))
1089 allocate(aero_data%num_ions(num_spec))
1090 allocate(aero_data%molec_weight(num_spec))
1091 allocate(aero_data%kappa(num_spec))
1092 allocate(aero_data%camp_particle_spec_id(num_spec))
1095 aero_data%i_water = 0
1097 do i_spec = 1, num_spec
1098 aero_data%name(i_spec) = spec_names(i_spec)%string
1099 if (.not. chem_spec_data%get_property_set( &
1100 aero_data%aero_rep_ptr%spec_name(
"P1." &
1101 // spec_names(i_spec)%string), property_set))
then
1102 call die_msg(934844845,
"Missing property set for aerosol species " &
1103 // spec_names(i_spec)%string)
1105 prop_name =
"density [kg m-3]"
1106 if (.not. property_set%get_real(prop_name, &
1107 aero_data%density(i_spec)))
then
1108 call die_msg(547508215,
"Missing density for aerosol species " &
1109 // spec_names(i_spec)%string)
1111 prop_name =
"num_ions"
1112 if (.not. property_set%get_int(prop_name, &
1113 aero_data%num_ions(i_spec)))
then
1114 call die_msg(324777059,
"Missing num_ions for aerosol species " &
1115 // spec_names(i_spec)%string)
1117 prop_name =
"molecular weight [kg mol-1]"
1118 if (.not. property_set%get_real(prop_name, &
1119 aero_data%molec_weight(i_spec)))
then
1120 call die_msg(549413749,
"Missing molec_weight for aerosol species " &
1121 // spec_names(i_spec)%string)
1124 if (.not. property_set%get_real(prop_name, &
1125 aero_data%kappa(i_spec)))
then
1126 call die_msg(944207343,
"Missing kappa for aerosol species " &
1127 // spec_names(i_spec)%string)
1129 prop_name =
"PartMC name"
1130 if (property_set%get_string(prop_name, str_val))
then
1131 if (str_val ==
"H2O")
then
1132 call assert_msg(227489086, aero_data%i_water == 0, &
1133 "Multiple aerosol water species")
1134 aero_data%i_water = i_spec
1137 aero_data%camp_particle_spec_id(i_spec) = &
1138 aero_data%aero_rep_ptr%spec_state_id(
"P1." &
1139 // spec_names(i_spec)%string)
1142 select type( aero_rep => aero_data%aero_rep_ptr)
1143 type is(aero_rep_single_particle_t)
1146 aero_data%camp_particle_state_size = aero_rep%per_particle_size()
1149 call die_msg(281737350,
"Wrong aerosol representation type")
1152 end subroutine aero_data_initialize