diff --git a/bld/build-namelist b/bld/build-namelist index 99e0b2d390..6806da76c0 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -696,24 +696,16 @@ if (defined $nl->get_value('sim_year')) { $sim_year =~ s/['"]//g; #"' } -# sim_year_start -# If sim_year is input as a range of years, then select the first year -# to use with some datasets -my $sim_year_start = $sim_year; -if ($sim_year =~ /(\d+)-(\d+)/) { - $sim_year_start = $1; -} - # Setup default ndep streams only if not simple_phys or aqua_mode and # the chemistry cannot produce the nitrogen depostion fluxes if (!($simple_phys or $aqua_mode)) { my $chem_nitrodep = chem_has_species($cfg, 'NO') and chem_has_species($cfg, 'NH3'); if ((!$chem_nitrodep) or ($chem =~ /geoschem/)) { - add_default($nl, 'stream_ndep_mesh_filename'); - add_default($nl, 'stream_ndep_data_filename', 'sim_year'=>$sim_year); - add_default($nl, 'stream_ndep_year_first', 'sim_year'=>$sim_year); - add_default($nl, 'stream_ndep_year_last', 'sim_year'=>$sim_year); - add_default($nl, 'stream_ndep_year_align', 'sim_year'=>$sim_year); + add_default($nl, 'stream_ndep_mesh_filename' , 'sim_year'=>$sim_year); + add_default($nl, 'stream_ndep_data_filename' , 'sim_year'=>$sim_year); + add_default($nl, 'stream_ndep_year_first' , 'sim_year'=>$sim_year); + add_default($nl, 'stream_ndep_year_last' , 'sim_year'=>$sim_year); + add_default($nl, 'stream_ndep_year_align' , 'sim_year'=>$sim_year); } } @@ -809,22 +801,32 @@ if ($cfg->get('cosp')) { # Carbon cycle constituents my $co2_cycle = $cfg->get('co2_cycle'); -if ($co2_cycle) { - - # co2_flag turns on the co2_cycle code in CAM - add_default($nl, 'co2_flag', 'val'=>'.true.'); +# co2_flag turns on the co2_cycle code in CAM +# defaults set in namelist_defaults_cam.xml +add_default($nl, 'co2_flag'); +if ($co2_cycle) { - # Supply a fossil fuel dataset if the co2_cycle is active and it's a - # transient run ... - if ($sim_year =~ /(\d+)-(\d+)/) { + # Supply a fossil fuel and aircraft emission datasets if the + # co2_cycle is active and it is a transient run + if ($sim_year =~ /(\d+)-(\d+)/ || $sim_year =~ /(\d+)/) { add_default($nl, 'co2_readflux_fuel', 'val'=>'.true.'); # Check whether user has explicitly turned off reading the fossil fuel dataset. # (user specification has higher precedence than the true value set above) if ($nl->get_value('co2_readflux_fuel') =~ /$TRUE/io) { - add_default($nl, 'co2flux_fuel_file', 'sim_year'=>$sim_year); + if ($sim_year =~ /(\d+)-(\d+)/) { + add_default($nl, 'co2flux_fuel_taxmode', 'val'=>'extend' ); + } else { + add_default($nl, 'co2flux_fuel_taxmode', 'val'=>'cycle' ); + } + add_default($nl, 'co2flux_fuel_tintalgo'); + add_default($nl, 'co2flux_fuel_meshfile'); + add_default($nl, 'co2flux_fuel_datafile'); + add_default($nl, 'co2flux_fuel_year_first', 'sim_year'=>$sim_year); + add_default($nl, 'co2flux_fuel_year_last' , 'sim_year'=>$sim_year); + add_default($nl, 'co2flux_fuel_year_align', 'sim_year'=>$sim_year); } add_default($nl, 'co2_readflux_aircraft', 'val'=>'.true.'); @@ -832,19 +834,26 @@ if ($co2_cycle) { # Check whether user has explicitly turned off reading the aircraft CO2 dataset. # (user specification has higher precedence than the true value set above) if ($nl->get_value('co2_readflux_aircraft') =~ /$TRUE/io) { - - my $rel_filepath = get_default_value('ac_CO2_emis'); - my $emisval = quote_string('ac_CO2 -> ' . $rel_filepath); - add_default($nl, 'aircraft_specifier', 'val'=>$emisval); - - add_default($nl, 'aircraft_datapath'); - add_default($nl, 'aircraft_type'); - # This should be the same file as the one in the aircraft_specifier file. - # This is a workaround to get this filepath into the cam.input_data_list file - # to allow the CESM scripts to obtain all required data for a run. - add_default($nl, 'aircraft_co2_file'); + if ($sim_year =~ /(\d+)-(\d+)/) { + add_default($nl, 'aircraft_co2_taxmode', 'val'=>'extend' ); + } else { + add_default($nl, 'aircraft_co2_taxmode', 'val'=>'cycle' ); + } + add_default($nl, 'aircraft_co2_tintalgo'); + add_default($nl, 'aircraft_co2_meshfile'); + add_default($nl, 'aircraft_co2_datafile'); + add_default($nl, 'aircraft_co2_year_first', 'sim_year'=>$sim_year); + add_default($nl, 'aircraft_co2_year_last' , 'sim_year'=>$sim_year); + add_default($nl, 'aircraft_co2_year_align', 'sim_year'=>$sim_year); } - } + } else { + # This case probably should not happen (co2_cycle on with no sim year) + add_default($nl, 'co2_readflux_fuel', 'val'=>'.false.'); + add_default($nl, 'co2_readflux_aircraft', 'val'=>'.false.'); + } +} else { + add_default($nl, 'co2_readflux_fuel', 'val'=>'.false.'); + add_default($nl, 'co2_readflux_aircraft', 'val'=>'.false.'); } # By default the prognostic co2_cycle CO2 will be radiative active, unless the diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 75983f00d8..9b8683127f 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -262,12 +262,6 @@ atm/cam/topo/USGS-gtopo30_4x5_remap_c050520.nc atm/cam/topo/fv_10x15_nc0540_Nsw042_Nrs008_Co060_Fi001_20171220.nc -atm/cam/topo/fv3_C24_nc3000_Co180_Fi001_MulG_PF_nullRR_Nsw127_c200625.nc -atm/cam/topo/fv3_C48_nc3000_Co120_Fi001_MulG_PF_nullRR_Nsw085_c200625.nc -atm/cam/topo/fv3_C96_nc3000_Co060_Fi001_MulG_PF_nullRR_Nsw042_c200625.nc -atm/cam/topo/fv3_C192_nc3000_Co030_Fi001_MulG_PF_Nsw021_c200625.nc -atm/cam/topo/fv3_C384_nc3000_Co015_Fi001_MulG_PF_nullRR_Nsw011_c200625.nc - atm/cam/topo/se/ne5np4_nc3000_Co360_Fi001_MulG_PF_nullRR_Nsw064_20170515.nc atm/cam/topo/se/ne16np4_nc3000_Co120_Fi001_PF_nullRR_Nsw084_20171012.nc atm/cam/topo/se/ne30np4_nc3000_Co060_Fi001_PF_nullRR_Nsw042_20171020.nc @@ -658,32 +652,34 @@ atm/waccm/lb/LBC_1765-2100_1.9x2.5_CCMI_RCP60_za_RNOCStrend_c141002.nc atm/waccm/lb/LBC_17500116-20150116_CMIP6_0p5degLat_c180905.nc +.false. +.true. + -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_1.9x2.5_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_1.9x2.5_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_1.9x2.5_c20181011.nc - -ac_CO2_filelist_175001-201512_fv_0.9x1.25_c20181011.txt -ac_CO2_filelist_175001-201512_fv_1.9x2.5_c20181011.txt -ac_CO2_filelist_175001-201512_fv_0.9x1.25_c20181011.txt -ac_CO2_filelist_175001-201512_fv_0.9x1.25_c20181011.txt -ac_CO2_filelist_175001-201512_fv_0.9x1.25_c20181011.txt -ac_CO2_filelist_175001-201512_fv_1.9x2.5_c20181011.txt -ac_CO2_filelist_175001-201512_fv_1.9x2.5_c20181011.txt -atm/cam/ggas -SERIAL - -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_1.9x2.5_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_0.9x1.25_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_1.9x2.5_c20181011.nc -atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_1.9x2.5_c20181011.nc +atm/cam/ggas/emissions-cmip6_CO2_anthro_surface_175001-201512_fv_0.9x1.25_c20181011.nc +share/meshes/fv0.9x1.25_141008_polemod_ESMFmesh.nc +2000 +2000 +1 +1850 +1850 +1 +1850 +2015 +1850 + + +atm/cam/ggas/emissions-cmip6_CO2_anthro_ac_175001-201512_fv_0.9x1.25_c20181011.nc +share/meshes/fv0.9x1.25_141008_polemod_ESMFmesh.nc +2000 +2000 +1 +1850 +1850 +1 +1850 +2015 +1850 atm/cam/ggas/noaamisc.r8.nc @@ -1981,12 +1977,6 @@ atm/cam/chem/trop_mam/atmsrf_ne30x4_ARCTIC_191011.nc atm/cam/chem/trop_mam/atmsrf_ne30x8_ARCTICGRIS_191212.nc -atm/cam/chem/trop_mam/atmsrf_C24_c200625.nc -atm/cam/chem/trop_mam/atmsrf_C48_c200625.nc -atm/cam/chem/trop_mam/atmsrf_C96_c200625.nc -atm/cam/chem/trop_mam/atmsrf_C192_c200625.nc -atm/cam/chem/trop_mam/atmsrf_C384_c200625.nc - atm/cam/chem/trop_mam/atmsrf_mpasa120_c090720.nc atm/cam/chem/trop_mam/atmsrf_mpasa480_c090720.nc @@ -2026,8 +2016,9 @@ share/meshes/fv0.9x1.25_141008_polemod_ESMFmesh.nc +lnd/clm2/ndepdata/fndep_clm_hist_b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensmean_1849-2015_monthly_0.9x1.25_c180926.nc -lnd/clm2/ndepdata/fndep_clm_hist_b.e21.BWHIST.f09_g17.CMIP6-historical-WACCM.ensmean_1849-2015_monthly_0.9x1.25_c180926.nc +share/meshes/fv0.9x1.25_141008_polemod_ESMFmesh.nc lnd/clm2/ndepdata/fndep_clm_WACCM6_CMIP6piControl001_y21-50avg_1850monthly_0.95x1.25_c180802.nc 2000 diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 3eaa0106cc..a39843f2d4 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -1778,34 +1778,133 @@ If TRUE turn on CO2 code. Default: set by build-namelist - -If TRUE read co2 fuel flux from file. -Default: set by build-namelist +If TRUE read co2 aircraft flux from file and use all the settings for +aircraft_co2_xxx namelist variables. +Default: TRUE for transient model runs, FALSE otherwise. - -If TRUE read co2 ocn flux from file. -Default: FALSE +If TRUE read co2 fuel flux from file and use all the settings for +co2flux_fuel_xxx namelist variables. +Default: TRUE for transient model runs, FALSE otherwise. - -If TRUE read co2 aircraft flux from file. + + + +Full filepath for dataset containing CO2 flux from fossil fuel. +Default: set by build-namelist, see namelist_defaults_cam.xml + + + +Full filepath for ESMF Mesh corresponding to co2flux_fuel_datafile File. Default: set by build-namelist - -Filepath for dataset containing CO2 flux from ocn. -Default: none + +First year to use in co2flux_fuel_datafile. +Default: set by build-namelist based on sim_year. + + +Last year to use in co2flux_fuel_datafile. +Default: set by build-namelist based on sim_year. - -Filepath for dataset containing CO2 flux from fossil fuel. -Default: none + +The simulation year corresponding to co2flux_fuel_year_first. +A common usage is to set this to the first year of the model run +(corresponding to the xml variable RUN_STARTDATE). With this setting, +the forcing in the first year of the run will be the forcing of year +yearFirst. +Another usage is to align the calendar of transient forcing with the +model calendar. For example, setting yearAlign = yearFirst will lead +to the forcing calendar being the same as the model calendar. The +forcing for a given model year would be the forcing of the same +year. This would be appropriate in transient runs where the model +calendar is setup to span the same year range as the forcing data. +Default set by build-namelist based on sim_year. + + + +Time interpolation algorithm to use for co2flux_fuel_datafile. +If not set by user, then set to 'nearest' if the variable 'time_bnds' is in +the forcing file otherwise, set to 'linear' + + + +Time extrapolation algorithm to use for co2flux_fuel_datafile. +- cycle: Simply cycle through the available model data (i.e., start over when the end is reached) +- extend: Use the final value when end of data is reached. +- limit: Halt the model run if the model time exceeds the available data. +Default: extend for transient, otherwise cycle + + + + + +Full pathname of the data file containing aircraft co2 emissions. +Default: set by build-namelist based on sim_year. + + + +Full pathname of the ESMF mesh file corresponding to aircraft_co2_datafile. +Default: set by build-namelist. + + + +First year of the aircraft_co2_datafile to use. +Default: set by build-namelist based on sim_year. + + + +Last year of the aircraft_co2_datafile to use. +Default: set by build-namelist based on sim_year. + + + +The simulation year corresponding to aircraft_co2_year_first. +A common usage is to set this to the first year of the model run +(corresponding to the xml variable RUN_STARTDATE). With this setting, +the forcing in the first year of the run will be the forcing of year +yearFirst. +Another usage is to align the calendar of transient forcing with the +model calendar. For example, setting yearAlign = yearFirst will lead +to the forcing calendar being the same as the model calendar. The +forcing for a given model year would be the forcing of the same +year. This would be appropriate in transient runs where the model +calendar is setup to span the same year range as the forcing data. +Default: set by build-namelist based on sim_year + + + +Time interpolation algorithm to use for aircraft_co2_datafile. +If not set by user, then set to 'nearest' if the variable 'time_bnds' is in +the forcing file otherwise, set to 'linear' + + + +Time extrapolation algorithm to use for aircraft_co2_datafile. +- cycle: Simply cycle through the available model data (i.e., start over when the end is reached) +- extend: Use the final value when end of data is reached. +- limit: Halt the model run if the model time exceeds the available data. +Default: extend for transient, otherwise cycle @@ -6936,35 +7035,127 @@ Full pathname of dataset of O2 and 03 column densities above the model for look- Default: set by build-namelist. - + + -Filename of file that contains aircraft input file lists. The filenames in the files are relative -to the directory specified by {{ hilight }}aircraft_datapath{{ closehilight }}. -Default: set by build-namelist. +Full pathname of the data file containing aircraft h2o emissions. - -Full pathname of the directory that contains the files specified in -{{ hilight }}aircraft_specifier{{ closehilight }}. -Default: set by build-namelist. +Full pathname of the ESMF mesh file corresponding to aircraft_h2o_datafile. - -Type of time interpolation for data in aircraft aerosol files. -Default: 'CYCLICAL_LIST' + + First year of the aircraft_h2o_datafile to use. + Default: NONE - -Full pathname of the ac_CO2 file specified in the filelist in -{{ hilight }}aircraft_specifier{{ closehilight }}. This is only to -get this name into the cam.input_data_list for the CESM scripts. -Default: set by build-namelist. + + Last year of the aircraft_h2o_datafile to use. + Default: NONE + +The simulation year corresponding to aircraft_h2o_year_first. +A common usage is to set this to the first year of the model run +(corresponding to the xml variable RUN_STARTDATE). With this setting, +the forcing in the first year of the run will be the forcing of year +yearFirst. +Another usage is to align the calendar of transient forcing with the +model calendar. For example, setting yearAlign = yearFirst will lead +to the forcing calendar being the same as the model calendar. The +forcing for a given model year would be the forcing of the same +year. This would be appropriate in transient runs where the model +calendar is setup to span the same year range as the forcing data. +Default: NONE + + +Time interpolation algorithm to use for aircraft_h2o_datafile. +If not set by user, then set to 'nearest' if the variable 'time_bnds' is in +the forcing file otherwise, set to 'linear' + + + +Time extrapolation algorithm to use for aircraft_h2o_datafile. +- cycle: Simply cycle through the available model data (i.e., start over when the end is reached) +- extend: Use the final value when end of data is reached. +- limit: Halt the model run if the model time exceeds the available data. +Default: extend for transient, otherwise cycle + + + + + + Full pathname of the data file containing ac_SLANT_DIST (by default). + Default: NONE + + + + Full pathname of the ESMF mesh file corresponding to aircraft_slant_dist. + Default: NONE + + + + Full pathname of the ESMF mesh file corresponding to aircraft_slant_dist. + Default: NONE + + + + First year of the aircraft_slant_dist to use. + Default: NONE + + + + Last year of the aircraft_slant_dist to use. + Default: NONE + + + +The simulation year corresponding to aircraft_slant_dist_year_first. +A common usage is to set this to the first year of the model run +(corresponding to the xml variable RUN_STARTDATE). With this setting, +the forcing in the first year of the run will be the forcing of year +yearFirst. +Another usage is to align the calendar of transient forcing with the +model calendar. For example, setting yearAlign = yearFirst will lead +to the forcing calendar being the same as the model calendar. The +forcing for a given model year would be the forcing of the same +year. This would be appropriate in transient runs where the model +calendar is setup to span the same year range as the forcing data. +Default: NONE + + + +Time interpolation algorithm to use for aircraft_slant_dist_datafile. +If not set by user, then set to 'nearest' if the variable 'time_bnds' is in +the forcing file otherwise, set to 'linear' + + + +Time extrapolation algorithm to use for aircraft_slant_dist. +- cycle: Simply cycle through the available model data (i.e., start over when the end is reached) +- extend: Use the final value when end of data is reached. +- limit: Halt the model run if the model time exceeds the available data. +Default: extend for transient, otherwise cycle + + + @@ -7024,7 +7215,7 @@ if {{ hilight }}gcr_ionization_type{{ closehilight }} is 'FIXED'. Default: 0 seconds - + @@ -8005,17 +8196,31 @@ Default: FALSE -Year first to use in nitrogen deposition stream data. + First year to use in nitrogen deposition stream data. + Default: 2000 -Year last to use in nitrogen deposition stream data. + Last year last to use in nitrogen deposition stream data. + Default: 2000 -Model year to align with stream_ndep_year_first. +The simulation year corresponding to stream_ndep_year_first. +A common usage is to set this to the first year of the model run +(corresponding to the xml variable RUN_STARTDATE). With this setting, +the forcing in the first year of the run will be the forcing of year +yearFirst. +Another usage is to align the calendar of transient forcing with the +model calendar. For example, setting yearAlign = yearFirst will lead +to the forcing calendar being the same as the model calendar. The +forcing for a given model year would be the forcing of the same +year. This would be appropriate in transient runs where the model +calendar is setup to span the same year range as the forcing data. +Set by build-namelist. +Default: 1 This varible is only used internally by build-namelist to determine appropriate defaults for climatological or transient forcing datasets. -Default: set by build-namelist. +Default: 2000 @@ -9900,38 +10105,6 @@ If TRUE perform temp tendency scaling before send to fv3 dynamics Default: FALSE - - - - - - - - -Stream filename(s) for Nitrogen Deposition data - - - -Stream meshfile for Nitrogen Deposition data - - - -First year to loop over for Nitrogen Deposition data - - - -Last year to loop over for Nitrogen Deposition data - - - -Simulation year that aligns with stream_year_first_ndep value - - diff --git a/bld/namelist_files/use_cases/1850_camnor_lt_ghgosloaero.xml b/bld/namelist_files/use_cases/1850_camnor_lt_ghgosloaero.xml index d8898106dc..5c352f3c35 100644 --- a/bld/namelist_files/use_cases/1850_camnor_lt_ghgosloaero.xml +++ b/bld/namelist_files/use_cases/1850_camnor_lt_ghgosloaero.xml @@ -98,4 +98,7 @@ 4 .true. + +1850 + diff --git a/bld/namelist_files/use_cases/1850_camnor_lt_osloaero.xml b/bld/namelist_files/use_cases/1850_camnor_lt_osloaero.xml index 0167860ad5..c414bb45d2 100644 --- a/bld/namelist_files/use_cases/1850_camnor_lt_osloaero.xml +++ b/bld/namelist_files/use_cases/1850_camnor_lt_osloaero.xml @@ -99,4 +99,7 @@ 4 .true. + +1850 + diff --git a/bld/namelist_files/use_cases/1850_camnor_lt_tropmam4.xml b/bld/namelist_files/use_cases/1850_camnor_lt_tropmam4.xml index 0dd5289c14..abea1bc3c0 100644 --- a/bld/namelist_files/use_cases/1850_camnor_lt_tropmam4.xml +++ b/bld/namelist_files/use_cases/1850_camnor_lt_tropmam4.xml @@ -78,4 +78,7 @@ 4 .true. + +1850 + diff --git a/bld/namelist_files/use_cases/2000_camnor_osloaero.xml b/bld/namelist_files/use_cases/2000_camnor_osloaero.xml index 9ba14f40b7..27e6ca6304 100644 --- a/bld/namelist_files/use_cases/2000_camnor_osloaero.xml +++ b/bld/namelist_files/use_cases/2000_camnor_osloaero.xml @@ -71,7 +71,7 @@ 4.60D0 -Leung_2023 +Leung_2023 @@ -91,4 +91,7 @@ 4 .true. + +2000 + diff --git a/bld/namelist_files/use_cases/2000_camnor_tropmam4.xml b/bld/namelist_files/use_cases/2000_camnor_tropmam4.xml index 08a4a04f6b..310ddb7673 100644 --- a/bld/namelist_files/use_cases/2000_camnor_tropmam4.xml +++ b/bld/namelist_files/use_cases/2000_camnor_tropmam4.xml @@ -101,4 +101,7 @@ 4 .true. + +2000 + diff --git a/bld/namelist_files/use_cases/hist_camnor_lt_osloaero.xml b/bld/namelist_files/use_cases/hist_camnor_lt_osloaero.xml index ef7f962161..5694cd424a 100644 --- a/bld/namelist_files/use_cases/hist_camnor_lt_osloaero.xml +++ b/bld/namelist_files/use_cases/hist_camnor_lt_osloaero.xml @@ -72,7 +72,7 @@ 4.60D0 -Leung_2023 +Leung_2023 @@ -92,4 +92,7 @@ 4 .true. + +1850-2015 + diff --git a/bld/namelist_files/use_cases/hist_camnor_lt_tropmam4.xml b/bld/namelist_files/use_cases/hist_camnor_lt_tropmam4.xml index 25e9e17d13..094971eeb3 100644 --- a/bld/namelist_files/use_cases/hist_camnor_lt_tropmam4.xml +++ b/bld/namelist_files/use_cases/hist_camnor_lt_tropmam4.xml @@ -60,4 +60,7 @@ 4 .true. + +1850-2015 + diff --git a/src/chemistry/utils/aircraft_emit.F90 b/src/chemistry/utils/aircraft_emit.F90 index f2a242aa63..2788fd0490 100644 --- a/src/chemistry/utils/aircraft_emit.F90 +++ b/src/chemistry/utils/aircraft_emit.F90 @@ -1,458 +1,693 @@ module aircraft_emit -!----------------------------------------------------------------------- -! -! Purpose: -! Manages reading and interpolation of aircraft aerosols -! -! Authors: Chih-Chieh (Jack) Chen and Cheryl Craig -- February 2010 -! -!----------------------------------------------------------------------- - use perf_mod, only : t_startf, t_stopf - - use shr_kind_mod, only : r8 => shr_kind_r8 - use cam_abortutils, only : endrun - use spmd_utils, only : masterproc - use tracer_data, only : trfld, trfile - use cam_logfile, only : iulog - - implicit none - private - save - - public :: aircraft_emit_init - public :: aircraft_emit_adv - public :: aircraft_emit_register - public :: aircraft_emit_readnl - public :: get_aircraft - - type :: forcing_air - real(r8) :: mw - character(len=256) :: filelist - character(len=256) :: filename - real(r8), pointer :: times(:) - real(r8), pointer :: levi(:) - character(len=11) :: species - character(len=8) :: units - integer :: nsectors - character(len=32),pointer :: sectors(:) - type(trfld),pointer :: fields(:) - type(trfile) :: file - end type forcing_air - - type(forcing_air),allocatable :: forcings_air(:) - - integer, parameter :: N_AERO = 13 - character(len=13) :: aero_names(N_AERO) = (/'ac_SLANT_DIST','ac_TRACK_DIST','ac_HC ','ac_NOX ','ac_PMNV ',& - 'ac_PMSO ','ac_PMFO ','ac_FUELBURN ','ac_CO2 ','ac_H2O ',& - 'ac_SOX ','ac_CO ','ac_BC '/) - - real(r8), parameter :: molmass(N_AERO) = 1._r8 - - logical :: advective_tracer(N_AERO) = .false. - character(len=3) :: mixtype(N_AERO) = 'wet' - - real(r8) :: cptmp = 666.0_r8 - real(r8) :: qmin = 0.0_r8 - logical :: cam_outfld = .false. - - integer :: index_map(N_AERO) - character(len=256) :: air_specifier(N_AERO)='' - character(len=256) :: air_datapath='' - character(len=24) :: air_type = 'CYCLICAL_LIST' ! 'CYCLICAL_LIST' - - logical :: rmv_file = .false. - - integer :: number_flds - - integer :: aircraft_cnt = 0 - character(len=16) :: spc_name_list(N_AERO) - character(len=256) :: spc_flist(N_AERO),spc_fname(N_AERO) - logical :: dist(N_AERO) - + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Manages reading and interpolation of aircraft aerosols + ! + ! Authors: + ! Chih-Chieh (Jack) Chen and Cheryl Craig -- February 2010 + ! Mariana Vertenstein (Refactored using CDEPS in-line functionality) -- November 2025 + ! + !----------------------------------------------------------------------- + + use perf_mod, only : t_startf, t_stopf + use shr_kind_mod, only : r8 => shr_kind_r8, cl=>shr_kind_cl, cs=>shr_kind_cs + use cam_abortutils, only : endrun + use cam_logfile, only : iulog + use spmd_utils, only : masterproc, iam + use dshr_strdata_mod, only : shr_strdata_type + + implicit none + private + + public :: aircraft_emit_init + public :: aircraft_emit_adv + public :: aircraft_emit_register + public :: aircraft_emit_readnl + public :: get_aircraft + + private :: get_vertical_dimension + private :: interpz_conserve + private :: chkrc + + type :: forcing_type + type(shr_strdata_type) :: sdat + character(len=cs) :: fldname = 'unset ' + character(len=cs) :: fldunits = 'unset' + character(len=cl) :: datafile = 'unset' + character(len=cl) :: meshfile = 'unset' + character(len=cs) :: mapalgo = 'consf' + character(len=cs) :: tintalgo = 'lower' + character(len=cs) :: taxmode = 'unset' + integer :: year_first = -999 + integer :: year_last = -999 + integer :: year_align = -999 + integer :: nilev = -1 + integer :: nlev = -1 + integer :: pbuf_index = -1 + real(r8), pointer :: altitude_int(:) + real(r8), pointer :: altitude_lev(:) + end type forcing_type + + integer, parameter :: N_AERO = 3 + type(forcing_type) :: forcing(N_AERO) + real(r8), parameter :: molmass(N_AERO) = 1._r8 + + logical :: first_update = .true. + + character(len=*),parameter :: u_FILE_u = __FILE__ + +!============================================================================ contains - - subroutine get_aircraft(cnt, spc_name_list_out) - integer, intent(out) :: cnt - character(len=16), optional, intent(out) :: spc_name_list_out(N_AERO) - integer :: i - - spc_name_list_out = '' - - cnt = aircraft_cnt - if( cnt>0 ) then - do i=1,cnt - spc_name_list_out(i) = spc_name_list(i) - end do - end if - - end subroutine get_aircraft - - subroutine aircraft_emit_register() - -!------------------------------------------------------------------ -! **** Add the aircraft aerosol data to the physics buffer **** -!------------------------------------------------------------------ - use ppgrid, only: pver,pcols - use physics_buffer, only : pbuf_add_field, dtype_r8 - use tracer_data, only: incr_filename - use constituents, only: cnst_add - - integer :: i,idx, mm, ind, n - character(len=16) :: spc_name - character(len=256) :: filelist, curr_filename - character(len=128) :: long_name - logical :: has_fixed_ubc=.false. - logical :: read_iv=.false. - - !------------------------------------------------------------------ - ! Return if air_specifier is blank (no aircraft data to process) - !------------------------------------------------------------------ - dist(:) = .false. - aircraft_cnt = 0 - if (air_specifier(1) == "") return - -! count aircraft emission species used in the simulation - count_emis: do n=1,N_AERO - - if( len_trim(air_specifier(n) ) == 0 ) then - exit count_emis - endif - - i = scan(air_specifier(n),'->') - spc_name = trim(adjustl(air_specifier(n)(:i-1))) - filelist = trim(adjustl(air_specifier(n)(i+2:))) - - mm = get_aircraft_ndx(spc_name) - if( mm < 1 ) then - call endrun('aircraft_emit_register: '//trim(spc_name)//' is not in the aircraft emission dataset') - endif - - if (trim(spc_name) == 'ac_SLANT_DIST'.or. trim(spc_name) == 'ac_TRACK_DIST') dist(n) = .true. - - aircraft_cnt = aircraft_cnt + 1 - call pbuf_add_field(aero_names(mm),'physpkg',dtype_r8,(/pcols,pver/),idx) - - spc_flist(aircraft_cnt) = filelist - spc_name_list(aircraft_cnt) = spc_name - index_map(aircraft_cnt) = mm - - curr_filename='' - spc_fname(aircraft_cnt) = incr_filename( curr_filename, filenames_list=spc_flist(aircraft_cnt), & - datapath=air_datapath) - - if( advective_tracer(mm) ) then - long_name = 'aircraft_'//trim(spc_name) - call cnst_add(aero_names(mm),molmass(mm),cptmp,qmin,ind,longname=long_name,readiv=read_iv, & - mixtype=mixtype(mm),cam_outfld=cam_outfld,fixed_ubc=has_fixed_ubc) - endif - - enddo count_emis -! count aircraft emission species used in the simulation - - endsubroutine aircraft_emit_register - - subroutine aircraft_emit_init() -!------------------------------------------------------------------- -! **** Initialize the aircraft aerosol data handling **** -!------------------------------------------------------------------- - use cam_history, only: addfld, add_default - use tracer_data, only: trcdata_init - use phys_control, only: phys_getopts - - implicit none - - character(len=16) :: spc_name - - integer :: astat, m - - logical :: history_chemistry - - call phys_getopts(history_chemistry_out=history_chemistry) - - !------------------------------------------------------------------ - ! Return if aircraft_cnt is zero (no aircraft data to process) - !------------------------------------------------------------------ - if (aircraft_cnt == 0 ) return - - if (masterproc) write(iulog,*) ' ' - - !----------------------------------------------------------------------- - ! allocate forcings type array - !----------------------------------------------------------------------- - allocate( forcings_air(aircraft_cnt), stat=astat ) - if( astat/= 0 ) then - write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air array; error = ',astat - call endrun - end if - - !----------------------------------------------------------------------- - ! setup the forcings_air type array - !----------------------------------------------------------------------- - species_loop : do m = 1,aircraft_cnt - - allocate( forcings_air(m)%sectors(1), stat=astat ) - if( astat/= 0 ) then - write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%sectors array; error = ',astat - call endrun - end if - - allocate( forcings_air(m)%fields(1), stat=astat ) - if( astat/= 0 ) then - write(iulog,*) 'aircraft_emit_init: failed to allocate forcings_air%fields array; error = ',astat - call endrun - end if - - spc_name = spc_name_list(m) - !----------------------------------------------------------------------- - ! default settings - !----------------------------------------------------------------------- - forcings_air(m)%file%stepTime = .true. ! Aircraft data is not to be interpolated in time - forcings_air(m)%file%cyclical_list = .true. ! Aircraft data cycles over the filename list - forcings_air(m)%file%weight_by_lat = .true. ! Aircraft data - interpolated with latitude weighting - forcings_air(m)%file%conserve_column = .true. ! Aircraft data - vertically interpolated to conserve the total column - forcings_air(m)%file%dist = dist(m) - forcings_air(m)%species = spc_name - forcings_air(m)%sectors = spc_name ! Only one species per file for aircraft data - forcings_air(m)%nsectors = 1 - forcings_air(m)%filelist = spc_flist(m) -! forcings_air(m)%file%curr_filename = spc_fname(m) - forcings_air(m)%filename = spc_fname(m) - end do species_loop - - if (masterproc) then - !----------------------------------------------------------------------- - ! diagnostics - !----------------------------------------------------------------------- - write(iulog,*) ' ' - write(iulog,*) 'aircraft_emit_init: diagnostics' - write(iulog,*) ' ' - write(iulog,*) 'aircraft_emit timing specs' - write(iulog,*) 'type = ',air_type - write(iulog,*) ' ' - write(iulog,*) 'there are ',aircraft_cnt,' species of aircraft emission' - do m = 1,aircraft_cnt - write(iulog,*) ' ' - write(iulog,*) 'forcing type ',m - write(iulog,*) 'species = ',trim(forcings_air(m)%species) - write(iulog,*) 'filelist= ',trim(forcings_air(m)%filelist) - end do - write(iulog,*) ' ' - endif - - !------------------------------------------------------------------ - ! Initialize the aircraft file processing - !------------------------------------------------------------------ - do m=1,aircraft_cnt - - allocate (forcings_air(m)%file%in_pbuf(size(forcings_air(m)%sectors))) - forcings_air(m)%file%in_pbuf(:) = .true. - call trcdata_init( forcings_air(m)%sectors, forcings_air(m)%filename, forcings_air(m)%filelist, air_datapath, & - forcings_air(m)%fields, forcings_air(m)%file, rmv_file, 0, 0, 0, air_type) - - - number_flds = 0 - if (associated(forcings_air(m)%fields)) number_flds = size( forcings_air(m)%fields ) - - if( number_flds < 1 ) then - if ( masterproc ) then - write(iulog,*) 'There are no aircraft aerosols' - write(iulog,*) ' ' - call endrun - endif - end if - - spc_name = spc_name_list(m) - call addfld( trim(spc_name), (/ 'lev' /), 'A', forcings_air(m)%fields(1)%units, & - 'aircraft emission '//trim(spc_name) ) - if (history_chemistry) then - call add_default( trim(spc_name), 1, ' ' ) - end if - end do - - - end subroutine aircraft_emit_init - - - - subroutine aircraft_emit_adv( state, pbuf2d) -!------------------------------------------------------------------- -! **** Advance to the next aircraft data **** -!------------------------------------------------------------------- - - use tracer_data, only : advance_trcdata - use physics_types,only : physics_state - use ppgrid, only : begchunk, endchunk - use ppgrid, only : pcols, pver - use string_utils, only : to_lower, GLC - use cam_history, only : outfld - use physconst, only : mwdry ! molecular weight dry air ~ kg/kmole - use physconst, only : boltz ! J/K/molecule -! C.-C. Chen -! use physconst, only : gravit, rearth - use phys_grid, only : get_wght_all_p - - use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk - - implicit none - - type(physics_state), intent(in) :: state(begchunk:endchunk) - - type(physics_buffer_desc), pointer :: pbuf2d(:,:) - - type(physics_buffer_desc), pointer :: pbuf_chnk(:) - integer :: ind,c,ncol,i,caseid,m,n - real(r8) :: to_mmr(pcols,pver) - real(r8),pointer :: tmpptr(:,:) - -! C.-C. Chen - real(r8) :: wght(pcols) - - !------------------------------------------------------------------ - ! Return if aircraft_cnt is zero (no aircraft data to process) - !------------------------------------------------------------------ - if (aircraft_cnt == 0 ) return - call t_startf('All_aircraft_emit_adv') - - !------------------------------------------------------------------- - ! For each field, read more data if needed and interpolate it to the current model time - !------------------------------------------------------------------- - do m = 1, aircraft_cnt - call advance_trcdata( forcings_air(m)%fields, forcings_air(m)%file, state, pbuf2d) - - !------------------------------------------------------------------- - ! set the tracer fields with the correct units - !------------------------------------------------------------------- - do i = 1,number_flds - -! C.-C. Chen, adding case 4 for kg/sec - select case ( to_lower(trim(forcings_air(m)%fields(i)%units(:GLC(forcings_air(m)%fields(i)%units)))) ) - case ("molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3") - caseid = 1 - case ('kg/kg','mmr') - caseid = 2 - case ('mol/mol','mole/mole','vmr','fraction') - caseid = 3 - case ('kg/kg/sec') - caseid = 4 - case ('kg m-2 s-1') - caseid = 5 - case ('m/sec' ) - caseid = 6 - case default - print*, 'aircraft_emit_adv: units = ',trim(forcings_air(m)%fields(i)%units) ,' are not recognized' - call endrun('aircraft_emit_adv: units are not recognized') - end select - - ind = index_map(i) - -!$OMP PARALLEL DO PRIVATE (C, NCOL, TO_MMR, tmpptr, pbuf_chnk) - do c = begchunk,endchunk - ncol = state(c)%ncol - -! C.-C. Chen, turning emission data to mixing ratio - call get_wght_all_p(c,ncol,wght(:ncol)) - - if (caseid == 1) then - to_mmr(:ncol,:) = (molmass(ind)*1.e6_r8*boltz*state(c)%t(:ncol,:))/(mwdry*state(c)%pmiddry(:ncol,:)) - elseif (caseid == 2) then - to_mmr(:ncol,:) = 1._r8 - elseif (caseid == 4) then -! do n=1,ncol -! to_mmr(n,:) = 1.0_r8/(rearth*rearth*wght(n)*state(c)%pdel(n,:)/gravit) -! end do - to_mmr(:ncol,:) = 1.0_r8 - elseif (caseid == 5) then - to_mmr(:ncol,:) = 1.0_r8 - elseif (caseid == 6) then - to_mmr(:ncol,:) = 1.0_r8 - else - to_mmr(:ncol,:) = molmass(ind)/mwdry - endif - pbuf_chnk => pbuf_get_chunk(pbuf2d, c) - call pbuf_get_field(pbuf_chnk, forcings_air(m)%fields(i)%pbuf_ndx, tmpptr ) - - tmpptr(:ncol,:) = tmpptr(:ncol,:)*to_mmr(:ncol,:) - - call outfld( forcings_air(m)%fields(i)%fldnam, & - tmpptr(:ncol,:), ncol, state(c)%lchnk ) - - enddo - enddo - enddo - - call t_stopf('All_aircraft_emit_adv') - end subroutine aircraft_emit_adv - - subroutine aircraft_emit_readnl(nlfile) -!------------------------------------------------------------------- -! **** Read in the aircraft_emit namelist ***** -!------------------------------------------------------------------- - use namelist_utils, only: find_group_name - use units, only: getunit, freeunit - use mpishorthand - - character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - - ! Local variables - integer :: unitn, ierr - character(len=*), parameter :: subname = 'aircraft_emit_readnl' - - character(len=256) :: aircraft_specifier(N_AERO) - character(len=256) :: aircraft_datapath - character(len=24) :: aircraft_type - - namelist /aircraft_emit_nl/ aircraft_specifier, aircraft_datapath, aircraft_type - !----------------------------------------------------------------------------- - - ! Initialize namelist variables from local module variables. - aircraft_specifier= air_specifier - aircraft_datapath = air_datapath - aircraft_type = air_type - - ! Read namelist - if (masterproc) then - unitn = getunit() - open( unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'aircraft_emit_nl', status=ierr) - if (ierr == 0) then - read(unitn, aircraft_emit_nl, iostat=ierr) - if (ierr /= 0) then - call endrun(subname // ':: ERROR reading namelist') +!============================================================================ + + subroutine aircraft_emit_readnl(nlfile) + + !------------------------------------------------------------------- + ! **** Read in the aircraft_emit namelist ***** + !------------------------------------------------------------------- + + use namelist_utils, only: find_group_name + use spmd_utils, only: mpicom, masterprocid + use spmd_utils, only: mpi_integer, mpi_logical, mpi_character + use co2_cycle, only: co2_readflux_aircraft + use cam_pio_utils, only: cam_pio_openfile + use string_utils, only: int2str + use pio, only: PIO_BCAST_ERROR, PIO_NOERR, PIO_NOWRITE + use pio, only: file_desc_t, pio_seterrorhandling, pio_inq_varid + use pio, only: pio_closefile + + ! Arguments + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! Local variables + integer :: nf, ni + integer :: unitn, ierr + type(file_desc_t) :: fileid + integer :: err_handling + integer :: varid + logical :: use_time_bnds + + character(len=cs) :: aircraft_co2_fldname = 'ac_CO2' + character(len=cl) :: aircraft_co2_datafile = 'unset' + character(len=cl) :: aircraft_co2_meshfile = 'unset' + character(len=cs) :: aircraft_co2_taxmode = 'unset' + character(len=cs) :: aircraft_co2_tintalgo = 'unset' + integer :: aircraft_co2_year_first = -999 + integer :: aircraft_co2_year_last = -999 + integer :: aircraft_co2_year_align = -999 + + character(len=cs) :: aircraft_h2o_fldname = 'ac_H2O' + character(len=cl) :: aircraft_h2o_datafile = 'unset' + character(len=cl) :: aircraft_h2o_meshfile = 'unset' + character(len=cs) :: aircraft_h2o_taxmode = 'unset' + character(len=cs) :: aircraft_h2o_tintalgo = 'unset' + integer :: aircraft_h2o_year_first = -999 + integer :: aircraft_h2o_year_last = -999 + integer :: aircraft_h2o_year_align = -999 + + character(len=cs) :: aircraft_slant_dist_fldname = 'ac_SLANT_DIST' + character(len=cl) :: aircraft_slant_dist_datafile = 'unset' + character(len=cl) :: aircraft_slant_dist_meshfile = 'unset' + character(len=cs) :: aircraft_slant_dist_tintalgo = 'unset' + character(len=cs) :: aircraft_slant_dist_taxmode = 'unset' + integer :: aircraft_slant_dist_year_first= -999 + integer :: aircraft_slant_dist_year_last = -999 + integer :: aircraft_slant_dist_year_align= -999 + + character(len=*), parameter :: subname = 'aircraft_emit_readnl' + + namelist /aircraft_emit_nl/ & + aircraft_co2_datafile, aircraft_co2_meshfile, & + aircraft_co2_year_first, aircraft_co2_year_last, aircraft_co2_year_align, & + aircraft_co2_taxmode, aircraft_co2_tintalgo, & + aircraft_h2o_datafile, aircraft_h2o_meshfile, & + aircraft_h2o_year_first, aircraft_h2o_year_last, aircraft_h2o_year_align, & + aircraft_h2o_taxmode, aircraft_h2o_tintalgo, & + aircraft_slant_dist_datafile, aircraft_slant_dist_meshfile, & + aircraft_slant_dist_year_first, aircraft_slant_dist_year_last, aircraft_slant_dist_year_align, & + aircraft_slant_dist_taxmode, aircraft_slant_dist_tintalgo + !----------------------------------------------------------------------------- + + ! Read namelist + if (masterproc) then + + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'aircraft_emit_nl', status=ierr) + if (ierr == 0) then + read(unitn, aircraft_emit_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading namelist') + end if + end if + close(unitn) + + ! Note - the following call assumes that co2_readflux_aircraft is + ! set in co2_cycle_readnl and this is called before this routine in + ! runtime_opts.F90. If co2_readflux_aircraft is .false. then, the + ! forcing(nf)%datafile = 'unset' and this logic will be triggered + ! in the other routines in this module + if (co2_readflux_aircraft) then + if (trim(aircraft_co2_datafile) /= 'unset') then + nf = 1 + forcing(nf)%fldname = aircraft_co2_fldname + forcing(nf)%datafile = aircraft_co2_datafile + forcing(nf)%meshfile = aircraft_co2_meshfile + forcing(nf)%year_first = aircraft_co2_year_first + forcing(nf)%year_last = aircraft_co2_year_last + forcing(nf)%year_align = aircraft_co2_year_align + forcing(nf)%taxmode = aircraft_co2_taxmode + forcing(nf)%tintalgo = aircraft_co2_tintalgo + end if + end if + if (trim(aircraft_h2o_datafile) /= 'unset') then + nf = 2 + forcing(nf)%datafile = aircraft_h2o_datafile + forcing(nf)%fldname = aircraft_h2o_fldname + forcing(nf)%meshfile = aircraft_h2o_meshfile + forcing(nf)%year_first = aircraft_h2o_year_first + forcing(nf)%year_last = aircraft_h2o_year_last + forcing(nf)%year_align = aircraft_h2o_year_align + forcing(nf)%taxmode = aircraft_h2o_taxmode + forcing(nf)%tintalgo = aircraft_h2o_tintalgo + end if + if (trim(aircraft_slant_dist_datafile) /= 'unset') then + nf = 3 + forcing(nf)%datafile = aircraft_slant_dist_datafile + forcing(nf)%fldname = aircraft_slant_dist_fldname + forcing(nf)%meshfile = aircraft_slant_dist_meshfile + forcing(nf)%year_first = aircraft_slant_dist_year_first + forcing(nf)%year_last = aircraft_slant_dist_year_last + forcing(nf)%year_align = aircraft_slant_dist_year_align + forcing(nf)%taxmode = aircraft_slant_dist_taxmode + forcing(nf)%tintalgo = aircraft_slant_dist_tintalgo end if + end if - close(unitn) - call freeunit(unitn) - end if -#ifdef SPMD - ! Broadcast namelist variables - call mpibcast(aircraft_specifier,len(aircraft_specifier(1))*N_AERO, mpichar, 0, mpicom) - call mpibcast(aircraft_datapath, len(aircraft_datapath), mpichar, 0, mpicom) - call mpibcast(aircraft_type, len(aircraft_type), mpichar, 0, mpicom) -#endif + n_aero_loop: do nf = 1,N_AERO + + ! Broadcast namelist variables + call mpi_bcast(forcing(nf)%datafile, len(forcing(nf)%datafile), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%datapath") + call mpi_bcast(forcing(nf)%fldname,len(forcing(nf)%fldname), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%fldname") + call mpi_bcast(forcing(nf)%meshfile, len(forcing(nf)%meshfile), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%meshfile") + call mpi_bcast(forcing(nf)%year_first, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%year_first") + call mpi_bcast(forcing(nf)%year_last, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%year_last") + call mpi_bcast(forcing(nf)%year_align, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%year_align") + call mpi_bcast(forcing(nf)%tintalgo, len(forcing(nf)%tintalgo), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%year_tintalgo") + call mpi_bcast(forcing(nf)%taxmode, len(forcing(nf)%taxmode), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: forcing("//int2str(nf)//"%year_taxmode") + + datafile_isnot_unset: if (trim(forcing(nf)%datafile) /= 'unset') then + ! overwrite mapalgo for ac_SLANT_DIST + if ( trim(forcing(nf)%fldname) == 'ac_SLANT_DIST') then + forcing(nf)%mapalgo = 'nn' + end if + + ! Overwrite forcing(nf)%tintalgo if it is set to 'unset' + ! Check if the data file has a time_bnds variable and if so set the time interpolation + ! type to 'nearest' otherwise set it to 'linear' + + if (trim(forcing(nf)%tintalgo) == 'unset') then + call cam_pio_openfile( fileid, forcing(nf)%datafile, PIO_NOWRITE ) + call pio_seterrorhandling( fileid, PIO_BCAST_ERROR, oldmethod=err_handling ) + ierr = pio_inq_varid( fileid, 'time_bnds', varid ) + call pio_seterrorhandling( fileid, err_handling) + use_time_bnds = (ierr == PIO_NOERR) + if (use_time_bnds) then + forcing(nf)%tintalgo = 'nearest' + else + forcing(nf)%tintalgo = 'linear' + end if + call pio_closefile( fileid ) + end if + + ! diagnostics + if (masterproc) then + write(iulog,*) ' ' + write(iulog,'(2a)' ) ' aircraft init settings for: ',trim(forcing(nf)%fldname) + write(iulog,'(2a)' ) ' aircraft datafile = ',trim(forcing(nf)%datafile) + write(iulog,'(2a)' ) ' aircraft meshfile = ',trim(forcing(nf)%meshfile) + write(iulog,'(2a)' ) ' aircraft mapalgo = ',trim(forcing(nf)%mapalgo) + write(iulog,'(2a)' ) ' aircraft tintalgo = ',trim(forcing(nf)%tintalgo) + write(iulog,'(2a)' ) ' aircraft taxmode = ',trim(forcing(nf)%taxmode) + write(iulog,'(a,i0)')' aircraft year_first = ',forcing(nf)%year_first + write(iulog,'(a,i0)')' aircraft year_last = ',forcing(nf)%year_last + write(iulog,'(a,i0)')' aircraft year_align = ',forcing(nf)%year_align + write(iulog,*) ' ' + end if + end if datafile_isnot_unset + + end do n_aero_loop + + end subroutine aircraft_emit_readnl + + !========================================================================= + subroutine aircraft_emit_register() + + !------------------------------------------------------------------ + ! **** Add the aircraft aerosol data to the physics buffer **** + !------------------------------------------------------------------ + use ppgrid, only: pver, pcols + use physics_buffer, only: pbuf_add_field, dtype_r8 + + ! Local variables + integer :: nf + !-------------------------------------------- + + do nf = 1,N_AERO + if (trim(forcing(nf)%datafile) /= 'unset') then + ! Add fldname to pbuf and obtain pbuf_index + call pbuf_add_field(forcing(nf)%fldname, 'physpkg', dtype_r8, (/pcols,pver/), & + forcing(nf)%pbuf_index) + end if + end do + + end subroutine aircraft_emit_register + + !========================================================================= + subroutine aircraft_emit_init() + + !------------------------------------------------------------------- + ! **** Initialize the aircraft aerosol data handling **** + !------------------------------------------------------------------- + use cam_history, only: addfld, add_default + use phys_control, only: phys_getopts + use cam_pio_utils, only: cam_pio_openfile, cam_pio_closefile + use pio, only: file_desc_t, var_desc_t + use pio, only: pio_inq_varid, pio_get_att + use pio, only: PIO_NOERR, PIO_NOWRITE + + ! Local variables + type(file_desc_t) :: pioid + type(var_desc_t) :: varid + integer :: ierr + integer :: klev + integer :: nf + logical :: history_chemistry + character(len=*), parameter :: subname = 'aircraft_emit_init' + !----------------------------------------------- + + call phys_getopts(history_chemistry_out=history_chemistry) + + loop_n_aero: do nf = 1,N_AERO + if (trim(forcing(nf)%datafile) /= 'unset') then + + ! Open file + call cam_pio_openfile( pioid, forcing(nf)%datafile, PIO_NOWRITE) + + ! Determine units + ierr = pio_inq_varid( pioid, forcing(nf)%fldname, varid ) + if (ierr/=pio_noerr) then + call endrun(trim(subname)//' Cannot find variable '//trim(forcing(nf)%fldname)// & + ' in file '//trim(forcing(nf)%datafile)) + endif + ierr = pio_get_att( pioid, varid, 'units', forcing(nf)%fldunits) + if (ierr/=pio_noerr) then + call endrun(trim(subname)//' Cannot get attribute units from '//trim(forcing(nf)%fldname)// & + ' in file '//trim(forcing(nf)%datafile)) + endif + + ! Determine vertical levels - altitude_int and altitude_lev + call get_vertical_dimension(fid=pioid, dname='altitude_int', dsize=forcing(nf)%nilev, & + data=forcing(nf)%altitude_int) + call get_vertical_dimension(fid=pioid, dname='altitude' , dsize=forcing(nf)%nlev , & + data=forcing(nf)%altitude_lev) + + ! Write out info to log file + if (masterproc) then + write(iulog,'(a)') trim(subname)// ' file: ',trim(forcing(nf)%datafile) + write(iulog,'(a)')' variable '//trim(forcing(nf)%fldname) + write(iulog,'(a)')' units '//trim(forcing(nf)%fldunits) + write(iulog,'(a)')' altitude levels ' + do klev=1,size(forcing(nf)%altitude_lev) + write(iulog,*)' ',klev,forcing(nf)%altitude_lev(klev) + end do + write(iulog,'(a)')' altitude interfaces ' + do klev=1,size(forcing(nf)%altitude_int) + write(iulog,*)' ',klev,forcing(nf)%altitude_int(klev) + end do + end if + + ! Close file + call cam_pio_closefile(pioid) + + ! Add field to cam history output + call addfld(trim(forcing(nf)%fldname), (/ 'lev' /), 'A', trim(forcing(nf)%fldunits), & + 'aircraft emission '//trim(forcing(nf)%fldname)) + if (history_chemistry) then + call add_default( trim(forcing(nf)%fldname), 1, ' ' ) + end if + + end if + end do loop_n_aero + + end subroutine aircraft_emit_init + + !========================================================================= + subroutine aircraft_emit_adv( state, pbuf2d ) + + !------------------------------------------------------------------- + ! **** Advance to the aircraft data **** + !------------------------------------------------------------------- + + use dshr_methods_mod , only: dshr_fldbun_getfldptr + use dshr_strdata_mod , only: shr_strdata_init_from_inline, shr_strdata_advance + use cam_esmf_mod, only: model_mesh, model_clock + use physics_types, only: physics_state + use ppgrid, only: begchunk, endchunk, pcols, pver, pverp + use string_utils, only: to_lower, GLC + use cam_history, only: outfld + use physconst, only: mwdry ! molecular weight dry air ~ kg/kmole + use physconst, only: boltz ! J/K/molecule + use phys_grid, only: get_wght_all_p, get_ncols_p + use physics_buffer, only: physics_buffer_desc, pbuf_get_field + use physics_buffer, only: pbuf_get_chunk + use string_utils, only: int2str + use time_manager, only: get_curr_date + + ! Arguments + type(physics_state), intent(in) :: state(begchunk:endchunk) + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + ! Local variables + integer :: gcell, nf + integer :: lchnk, icol, klev, ncol + integer :: caseid + integer :: year, mon, day, sec + integer :: mcdate + real(r8) :: to_mmr(pcols,pver) + real(r8) :: wght(pcols) + real(r8), pointer :: tmpptr(:,:) + real(r8), pointer :: data_out(:,:) + real(r8), pointer :: dataptr2d(:,:) + real(r8), allocatable :: datain3d(:,:,:) + real(r8) :: data_col(pver) + real(r8) :: model_z(pverp) + character(len=cs) :: units + integer :: rc + type(physics_buffer_desc), pointer :: pbuf_chnk(:) + real(r8), parameter :: m2km = 1.e-3_r8 + character(len=*), parameter :: subname = 'aircraft_emit_adv' + !------------------------------------------------------------------ + + call t_startf('All_aircraft_emit_adv') + + !------------------------------------------------------------------ + ! The stream initialization must be called after the cam initailization + !------------------------------------------------------------------ + n_aero_loop: do nf = 1,N_AERO + unset_file: if (trim(forcing(nf)%datafile) /= 'unset') then + + first_call: if (first_update) then + ! Initialize forcing%sdat + call shr_strdata_init_from_inline(forcing(nf)%sdat, & + my_task = iam, & + logunit = iulog, & + compname = 'ATM', & + model_clock = model_clock, & + model_mesh = model_mesh, & + stream_meshfile = trim(forcing(nf)%meshfile), & + stream_filenames = (/forcing(nf)%datafile/), & + stream_yearFirst = forcing(nf)%year_first, & + stream_yearLast = forcing(nf)%year_last, & + stream_yearAlign = forcing(nf)%year_align, & + stream_fldlistFile = (/forcing(nf)%fldname/), & + stream_fldListModel = (/forcing(nf)%fldname/), & + stream_lev_dimname = 'altitude', & + stream_mapalgo = trim(forcing(nf)%mapalgo), & + stream_offset = 0, & + stream_taxmode = trim(forcing(nf)%taxmode), & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = trim(forcing(nf)%tintalgo), & + stream_name = 'Aircraft forcing data ', & + rc = rc) + call chkrc(rc,__LINE__,u_FILE_u) + + first_update = .false. + end if first_call + + !------------------------------------------------------------------- + ! For each field, interpolate data in time and to model horizontal grid + !------------------------------------------------------------------- + + ! Extract YMD from model_update_next_time + call get_curr_date(year, mon, day, sec) + mcdate = year*10000 + mon*100 + day + + ! Advance sdat streams + call shr_strdata_advance(forcing(nf)%sdat, ymd=mcdate, tod=sec, logunit=iulog, & + istr='aircraft_stream', rc=rc) + call chkrc(rc,__LINE__,u_FILE_u) + + ! Get pointer to horizontally interpolated data + call dshr_fldbun_getFldPtr(forcing(nf)%sdat%pstrm(1)%fldbun_model, trim(forcing(nf)%fldname), & + fldptr2=dataptr2d, rc=rc) + call chkrc(rc,__LINE__,u_FILE_u) + + ! Obtain datain on model horizontal grid but the same vertical levels as the forcing dataset + allocate(datain3d(pcols,pver,begchunk:endchunk), stat=rc) + if ( rc /= 0 ) then + call endrun(trim(subname)//': failed to allocate datain3d, error = '//int2str(rc)) + end if + do klev = 1, forcing(nf)%nlev !nlev is the number of levels in the forcing data + gcell = 1 + do lchnk = begchunk,endchunk + ncol = get_ncols_p(lchnk) + do icol = 1,ncol + datain3d(icol,klev,lchnk) = dataptr2d(klev,gcell) + gcell = gcell + 1 + end do + end do + end do + + ! Do vertical interpolation - aircraft data is vertically + ! interpolated to conserve the total column + do lchnk = begchunk,endchunk + call pbuf_get_field(pbuf2d, lchnk, forcing(nf)%pbuf_index, data_out) + ncol = get_ncols_p(lchnk) + do icol = 1,ncol + model_z(1:pverp) = m2km * state(lchnk)%zi(icol,pverp:1:-1) + call interpz_conserve( forcing(nf)%nlev, pver, forcing(nf)%altitude_int, model_z, & + datain3d(icol,:,lchnk), data_col(:) ) + data_out(icol,:) = data_col(pver:1:-1) + end do + end do + deallocate(datain3d) + + !------------------------------------------------------------------- + ! set the tracer fields with the correct units + !------------------------------------------------------------------- + + ! GLC is position of last significant character in string. + units = to_lower(trim(forcing(nf)%fldunits(:GLC(forcing(nf)%fldunits)))) + select case (trim(units)) + case ("molec/cm3","/cm3","molecules/cm3","cm^-3","cm**-3") + caseid = 1 + case ('kg/kg','mmr') + caseid = 2 + case ('mol/mol','mole/mole','vmr','fraction') + caseid = 3 + case ('kg/kg/sec') + caseid = 4 + case ('kg m-2 s-1') + caseid = 5 + case ('m/sec' ) + caseid = 6 + case default + if (masterproc) then + write(iulog,*)'aircraft_emit_adv: units = '//trim(units)//' are not recognized' + end if + call endrun(trim(subname)//' aircraft_emit_adv: units are not recognized') + end select + + !$OMP PARALLEL DO PRIVATE (lchnk, ncol, to_mmr, tmpptr, pbuf_chnk, wght) + do lchnk = begchunk,endchunk + ncol = state(lchnk)%ncol + + ! Turn emission data to mixing ratio + call get_wght_all_p(lchnk, ncol, wght(:ncol)) + + if (caseid == 1) then + to_mmr(:ncol,:) = (molmass(nf)*1.e6_r8*boltz*state(lchnk)%t(:ncol,:)) & + /(mwdry*state(lchnk)%pmiddry(:ncol,:)) + elseif (caseid == 2) then + to_mmr(:ncol,:) = 1._r8 + elseif (caseid == 4) then + to_mmr(:ncol,:) = 1.0_r8 + elseif (caseid == 5) then + to_mmr(:ncol,:) = 1.0_r8 + elseif (caseid == 6) then + to_mmr(:ncol,:) = 1.0_r8 + else + to_mmr(:ncol,:) = molmass(nf)/mwdry + endif + + pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk) + call pbuf_get_field(pbuf_chnk, forcing(nf)%pbuf_index, tmpptr) + tmpptr(:ncol,:) = tmpptr(:ncol,:)*to_mmr(:ncol,:) + call outfld( forcing(nf)%fldname, tmpptr(:ncol,:), ncol, state(lchnk)%lchnk ) + enddo + + end if unset_file + end do n_aero_loop + + call t_stopf('All_aircraft_emit_adv') + + end subroutine aircraft_emit_adv + + !========================================================================= + subroutine interpz_conserve( nsrc, ndst, src_x, dst_x, src, dst) + ! Note, this routine is an edited version of the tracer_data version + + ! Arguments + integer, intent(in) :: nsrc ! dimension source array + integer, intent(in) :: ndst ! dimension target array + real(r8), intent(in) :: src_x(nsrc+1) ! source coordinates + real(r8), intent(in) :: dst_x(ndst+1) ! target coordinates + real(r8), intent(in) :: src(nsrc) ! source array + real(r8), intent(out) :: dst(ndst) ! target array + + ! local variables + integer :: i, j + integer :: isrc + real(r8) :: y + real(r8) :: bot, top + !--------------------------------------------------------------- + + do i = 1, ndst + if ( (dst_x(i)src_x(1)) ) then + do isrc = 1,nsrc + if ( (dst_x(i)-src_x(isrc))*(dst_x(i)-src_x(isrc+1))<=0.0_r8 ) then + exit + end if + end do + + if ( dst_x(i)src_x(j+1) ) then + y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j)) + bot = src_x(j+1) + else + y = y+(top-bot)*src(j)/(src_x(j+1)-src_x(j)) + exit + endif + enddo + dst(i) = y + else + dst(i) = 0.0_r8 + end if + end do + + if ( dst_x(1)>src_x(1) ) then + top = dst_x(1) + bot = src_x(1) + y = 0.0_r8 + do j = 1, nsrc + if ( top>src_x(j+1) ) then + y = y+(src_x(j+1)-bot)*src(j)/(src_x(j+1)-src_x(j)) + bot = src_x(j+1) + else + y = y+(top-bot)*src(j)/(src_x(j+1)-src_x(j)) + exit + endif + end do + dst(1) = dst(1)+y + end if - ! Update module variables with user settings. - air_specifier = aircraft_specifier - air_datapath = aircraft_datapath - air_type = aircraft_type + end subroutine interpz_conserve - end subroutine aircraft_emit_readnl + !========================================================================= + subroutine get_aircraft(cnt, spc_name_list_out) - integer function get_aircraft_ndx( name ) + ! Arguments + integer, intent(out) :: cnt + character(len=*), intent(out) :: spc_name_list_out(:) - implicit none - character(len=*), intent(in) :: name + ! Local variables + integer :: nf + !------------------------------------------------------------------ - integer :: i + cnt = 0 + spc_name_list_out(:) = '' - get_aircraft_ndx = 0 - do i = 1,N_AERO - if ( trim(name) == trim(aero_names(i)) ) then - get_aircraft_ndx = i - return + do nf = 1,N_AERO + if (trim(forcing(nf)%datafile) /= 'unset') then + cnt = cnt + 1 + spc_name_list_out(nf) = trim(forcing(nf)%fldname) + end if + end do + + end subroutine get_aircraft + + !========================================================================= + subroutine get_vertical_dimension( fid, dname, dsize, data ) + + use pio, only: file_desc_t, pio_seterrorhandling + use pio, only: pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, pio_get_var + use pio, only: PIO_BCAST_ERROR, PIO_NOERR + use string_utils, only: int2str + + ! Arguments + type(file_desc_t), intent(inout) :: fid + character(*), intent(in) :: dname + integer, intent(out) :: dsize + real(r8), pointer :: data(:) + + ! Local variables + integer :: vid, ierr, id + integer :: err_handling + character(len=*), parameter :: subname = 'get_vertical_dimension' + !------------------------------------------------------------------ + + call pio_seterrorhandling(fid, PIO_BCAST_ERROR, oldmethod=err_handling) + ierr = pio_inq_dimid( fid, dname, id ) + if ( ierr == PIO_NOERR ) then + ierr = pio_inq_dimlen( fid, id, dsize ) + if (ierr /= PIO_NOERR) then + call endrun(trim(subname)//': failed on pio_inq_dimid') + end if + allocate( data(dsize), stat=ierr ) + if ( ierr /= 0 ) then + call endrun(trim(subname)//': failed to allocate data array, error = '//int2str(ierr)) + end if + ierr = pio_inq_varid( fid, dname, vid ) + if (ierr /= PIO_NOERR) then + call endrun(trim(subname)//': failed on pio_inq_varid') + end if + ierr = pio_get_var( fid, vid, data ) + if (ierr /= PIO_NOERR) then + call endrun(trim(subname)//': failed on pio_get_var') + end if endif - enddo + call pio_seterrorhandling(fid, err_handling) + + end subroutine get_vertical_dimension + + !================================================================ + subroutine chkrc(rc, line, file) + use ESMF, only: ESMF_LOGMSG_ERROR, ESMF_SUCCESS, ESMF_LogWrite - end function get_aircraft_ndx + ! Arguments + integer , intent(in) :: rc + integer , intent(in) :: line + character(len=*) , intent(in) :: file + + if ( rc /= ESMF_SUCCESS ) then + call ESMF_LogWrite('ERROR:', ESMF_LOGMSG_ERROR, line=line, file=file) + call endrun('chkrc') + end if + end subroutine chkrc end module aircraft_emit diff --git a/src/chemistry/utils/tracer_data.F90 b/src/chemistry/utils/tracer_data.F90 index ce6843eafa..578e3b0951 100644 --- a/src/chemistry/utils/tracer_data.F90 +++ b/src/chemistry/utils/tracer_data.F90 @@ -26,11 +26,11 @@ module tracer_data pio_get_var, pio_get_att, pio_nowrite, pio_inq_dimlen, & pio_inq_vardimid, pio_inq_dimlen, pio_closefile, & pio_inquire_variable + use string_utils, only: int2str implicit none private ! all unless made public - save public :: trfld, input3d, input2d, trfile public :: trcdata_init @@ -170,8 +170,6 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & use horizontal_interpolate, only : xy_interp_init use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8, mpi_integer - implicit none - character(len=*), intent(in) :: specifier(:) character(len=*), intent(in) :: filename character(len=*), intent(in) :: filelist @@ -237,8 +235,8 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & file%cyc_yr = data_cycle_yr case( 'SERIAL' ) case default - write(iulog,*) 'trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename) - write(iulog,*) 'trcdata_init: valid data types: SERIAL | CYCLICAL | CYCLICAL_LIST | FIXED | INTERP_MISSING_MONTHS ' + write(iulog,'(4a)') 'trcdata_init: invalid data type: ', trim(data_type), ' file: ', trim(filename) + write(iulog,'(a)') 'trcdata_init: valid data types: SERIAL | CYCLICAL | CYCLICAL_LIST | FIXED | INTERP_MISSING_MONTHS ' call endrun('trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename)) endselect @@ -254,7 +252,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & end if if (masterproc) then - write(iulog,*) 'trcdata_init: data type: '//trim(data_type)//' file: '//trim(filename) + write(iulog,'(4a)') 'trcdata_init: data type: ', trim(data_type), ' file: ', trim(filename) endif ! if there is no list of files (len_trim(file%filenames_list)<1) then @@ -376,7 +374,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & endif if (masterproc) then - write(iulog,*) 'trcdata_init: file%has_ps = ' , file%has_ps + write(iulog,'(a,l4)') 'trcdata_init: file%has_ps = ' , file%has_ps endif ! masterproc if (file%alt_data) then @@ -397,14 +395,12 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & allocate( file%hyam(file%nlev), file%hybm(file%nlev), stat=astat ) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%hyam,file%hybm allocation error = ',astat - call endrun('trcdata_init: failed to allocate file%hyam and file%hybm arrays') + call endrun('trcdata_init: failed to allocate file%hyam and file%hybm arrays, error code = '//int2str(astat)) end if allocate( file%hyai(file%nlev+1), file%hybi(file%nlev+1), stat=astat ) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%hyai,file%hybi allocation error = ',astat - call endrun('trcdata_init: failed to allocate file%hyai and file%hybi arrays') + call endrun('trcdata_init: failed to allocate file%hyai and file%hybi arrays, error code = '//int2str(astat)) end if call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR, oldmethod=err_handling) @@ -429,24 +425,20 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & allocate( file%ps_in(1)%data(pcols,begchunk:endchunk), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(1)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate file%ps_in(1)%data array; error = '//int2str(astat)) end if allocate( file%ps_in(2)%data(pcols,begchunk:endchunk), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(2)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate file%ps_in(2)%data array; error = '//int2str(astat)) end if if( file%fill_in_months ) then allocate( file%ps_in(3)%data(pcols,begchunk:endchunk), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(3)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate file%ps_in(3)%data array; error = '//int2str(astat)) end if allocate( file%ps_in(4)%data(pcols,begchunk:endchunk), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(4)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate file%ps_in(4)%data array; error = '//int2str(astat)) end if end if endif @@ -459,8 +451,7 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & ! get netcdf variable id for the field ierr = pio_inq_varid( file%curr_fileid, flds(f)%srcnam, flds(f)%var_id ) if (ierr/=pio_noerr) then - call endrun('trcdata_init: Cannot find var "'//trim(flds(f)%srcnam)// & - '" in file "'//trim(file%curr_filename)//'"') + call endrun('trcdata_init: Cannot find var "'//trim(flds(f)%srcnam)//'" in file "'//trim(file%curr_filename)//'"') endif ! determine if the field has a vertical dimension @@ -482,21 +473,19 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & allocate( flds(f)%data(pcols,pver,begchunk:endchunk), stat=astat ) endif if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate flds(f)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate flds(f)%data array; error = '//int2str(astat)) end if else flds(f)%pbuf_ndx = pbuf_get_index(flds(f)%fldnam,errcode) endif if (flds(f)%srf_fld) then - allocate( flds(f)%input(1)%data(pcols,1,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(1)%data(pcols,1,begchunk:endchunk), stat=astat) else - allocate( flds(f)%input(1)%data(pcols,file%nlev,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(1)%data(pcols,file%nlev,begchunk:endchunk), stat=astat) endif if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(1)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate flds(f)%input(1)%data array; error = '//int2str(astat)) end if if (flds(f)%srf_fld) then allocate( flds(f)%input(2)%data(pcols,1,begchunk:endchunk), stat=astat ) @@ -504,28 +493,25 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & allocate( flds(f)%input(2)%data(pcols,file%nlev,begchunk:endchunk), stat=astat ) endif if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(2)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate flds(f)%input(2)%data array; error = '//int2str(astat)) end if if( file%fill_in_months ) then if (flds(f)%srf_fld) then - allocate( flds(f)%input(3)%data(pcols,1,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(3)%data(pcols,1,begchunk:endchunk), stat=astat) else - allocate( flds(f)%input(3)%data(pcols,file%nlev,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(3)%data(pcols,file%nlev,begchunk:endchunk), stat=astat) endif if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(3)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate flds(f)%input(3)%data array; error = '//int2str(astat)) end if if (flds(f)%srf_fld) then - allocate( flds(f)%input(4)%data(pcols,1,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(4)%data(pcols,1,begchunk:endchunk), stat=astat) else - allocate( flds(f)%input(4)%data(pcols,file%nlev,begchunk:endchunk), stat=astat ) + allocate( flds(f)%input(4)%data(pcols,file%nlev,begchunk:endchunk), stat=astat) endif if( astat/= 0 ) then - write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(4)%data array; error = ',astat - call endrun + call endrun('trcdata_init: failed to allocate flds(f)%input(4)%data array; error = '//int2str(astat)) end if endif @@ -625,36 +611,30 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & deallocate(phi,lam) -! weight_x & weight_y are weighting function for x & y interpolation + ! weight_x & weight_y are weighting function for x & y interpolation allocate(file%weight_x(plon,file%nlon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%weight_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate weight_x array') + call endrun('trcdata_init: file%weight_x allocation error = '//int2str(astat)) end if allocate(file%weight_y(plat,file%nlat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%weight_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate weight_y array') + call endrun('trcdata_init: file%weight_y allocation error = '//int2str(astat)) end if allocate(file%count_x(plon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%count_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate count_x array') + call endrun('trcdata_init: file%count_x allocation error = '//int2str(astat)) end if allocate(file%count_y(plat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%count_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate count_y array') + call endrun('trcdata_init: file%count_y allocation error = '//int2str(astat)) end if allocate(file%index_x(plon,file%nlon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%index_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate index_x array') + call endrun('trcdata_init: file%index_x allocation error = '//int2str(astat)) end if allocate(file%index_y(plat,file%nlat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%index_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate index_y array') + call endrun('trcdata_init: file%index_y allocation error = '//int2str(astat)) end if file%weight_x(:,:) = 0.0_r8 file%weight_y(:,:) = 0.0_r8 @@ -666,33 +646,27 @@ subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & if( file%dist ) then allocate(file%weight0_x(plon,file%nlon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%weight0_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate weight0_x array') + call endrun('trcdata_init: file%weight0_x allocation error = '//int2str(astat)) end if allocate(file%weight0_y(plat,file%nlat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%weight0_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate weight0_y array') + call endrun('trcdata_init: file%weight0_y allocation error = '//int2str(astat)) end if allocate(file%count0_x(plon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%count0_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate count0_x array') + call endrun('trcdata_init: file%count0_x allocation error = '//int2str(astat)) end if allocate(file%count0_y(plat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%count0_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate count0_y array') + call endrun('trcdata_init: file%count0_y allocation error = '//int2str(astat)) end if allocate(file%index0_x(plon,file%nlon), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%index0_x allocation error = ',astat - call endrun('trcdata_init: failed to allocate index0_x array') + call endrun('trcdata_init: file%index0_x allocation error = '//int2str(astat)) end if allocate(file%index0_y(plat,file%nlat), stat=astat) if( astat /= 0 ) then - write(iulog,*) 'trcdata_init: file%index0_y allocation error = ',astat - call endrun('trcdata_init: failed to allocate index0_y array') + call endrun('trcdata_init: file%index0_y allocation error = '//int2str(astat)) end if file%weight0_x(:,:) = 0.0_r8 file%weight0_y(:,:) = 0.0_r8 @@ -794,8 +768,6 @@ end subroutine trcdata_init subroutine advance_trcdata( flds, file, state, pbuf2d ) use physics_types,only : physics_state - implicit none - type(trfile), intent(inout) :: file type(trfld), intent(inout) :: flds(:) type(physics_state), intent(in) :: state(begchunk:endchunk) @@ -824,7 +796,7 @@ subroutine advance_trcdata( flds, file, state, pbuf2d ) call t_startf('read_next_trcdata') call read_next_trcdata( flds, file ) call t_stopf('read_next_trcdata') - if(masterproc) write(iulog,*) 'READ_NEXT_TRCDATA ', flds%fldnam + if(masterproc) write(iulog,'(2a)') 'READ_NEXT_TRCDATA ',flds%fldnam end if endif @@ -845,9 +817,6 @@ end subroutine advance_trcdata !------------------------------------------------------------------- subroutine get_fld_data( flds, field_name, data, ncol, lchnk, pbuf ) - - implicit none - type(trfld), intent(inout) :: flds(:) character(len=*), intent(in) :: field_name real(r8), intent(out) :: data(:,:) @@ -879,8 +848,6 @@ end subroutine get_fld_data !------------------------------------------------------------------- subroutine get_fld_ndx( flds, field_name, idx ) - implicit none - type(trfld), intent(in) :: flds(:) character(len=*), intent(in) :: field_name integer, intent(out) :: idx @@ -901,7 +868,7 @@ end subroutine get_fld_ndx !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ subroutine get_model_time(file) - implicit none + type(trfile), intent(inout) :: file integer yr, mon, day, ncsec ! components of a date @@ -918,8 +885,6 @@ end subroutine get_model_time !------------------------------------------------------------------------------ subroutine check_files( file, fids, itms, times_found) - implicit none - type(trfile), intent(inout) :: file type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs integer, optional, intent(out) :: itms(2) @@ -997,9 +962,6 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ !----------------------------------------------------------------------- use string_utils, only : incstr - use shr_file_mod, only : shr_file_getunit, shr_file_freeunit - - implicit none character(len=*), intent(in) :: filename ! present dynamical dataset filename character(len=*), optional, intent(in) :: filenames_list @@ -1034,15 +996,13 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ !----------------------------------------------------------------------- pos = len_trim( filename ) fn_new = filename(:pos) - if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(fn_new) + if ( masterproc ) write(iulog,'(2a)') 'incr_flnm: old filename = ', trim(fn_new) if( fn_new(pos-2:) == '.nc' ) then pos = pos - 3 end if istat = incstr( fn_new(:pos), 1 ) if( istat /= 0 ) then - write(iulog,*) 'incr_flnm: incstr returned ', istat - write(iulog,*) ' while trying to decrement ',trim( fn_new ) - call endrun + call endrun('incr_flnm: incstr returned '//int2str(istat)//' while trying to decrement '//trim(fn_new)) end if else @@ -1050,16 +1010,17 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ !------------------------------------------------------------------- ! ... open filenames_list !------------------------------------------------------------------- - if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(filename) - if ( masterproc ) write(iulog,*) 'incr_flnm: open filenames_list : ',trim(filenames_list) - unitnumber = shr_file_getUnit() + if ( masterproc ) then + write(iulog,'(2a)') 'incr_flnm: old filename = ', trim(filename) + write(iulog,'(2a)') 'incr_flnm: open filenames_list : ', trim(filenames_list) + end if if ( present(datapath) ) then filepath = trim(datapath) //'/'// trim(filenames_list) else filepath = trim(filenames_list) endif - open( unit=unitnumber, file=filepath, iostat=ios, status="OLD") + open( newunit=unitnumber, file=filepath, iostat=ios, status="OLD") if (ios /= 0) then call endrun('not able to open file: '//trim(filepath)) endif @@ -1075,8 +1036,8 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ fn_new = 'NOT_FOUND' incr_filename = trim(fn_new) return - endif - endif + end if + end if !------------------------------------------------------------------- ! If current filename is '', then initialize with the first filename read in @@ -1096,9 +1057,9 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ fn_new = 'NOT_FOUND' incr_filename = trim(fn_new) return - endif - endif - enddo + end if + end if + end do !------------------------------------------------------------------- ! Read next filename @@ -1142,14 +1103,13 @@ function incr_filename( filename, filenames_list, datapath, cyclical_list, list_ fn_new = trim(line) close(unit=unitnumber) - call shr_file_freeUnit(unitnumber) - endif + end if !--------------------------------------------------------------------------------- ! return the current filename !--------------------------------------------------------------------------------- incr_filename = trim(fn_new) - if ( masterproc ) write(iulog,*) 'incr_flnm: new filename = ',trim(incr_filename) + if ( masterproc ) write(iulog,'(2a)') 'incr_flnm: new filename = ', trim(incr_filename) end function incr_filename @@ -1159,8 +1119,6 @@ subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found use intp_util, only: findplb - implicit none - type(trfile), intent(in) :: file real(r8), intent(out) :: datatimem, datatimep @@ -1186,8 +1144,7 @@ subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found allocate( all_data_times( all_tsize ), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'find_times: failed to allocate all_data_times array; error = ',astat - call endrun + call endrun('find_times: failed to allocate all_data_times array; error = '//int2str(astat)) end if all_data_times(:curr_tsize) = file%curr_data_times(:) @@ -1238,7 +1195,7 @@ subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found if ( .not. times_found ) then if (masterproc) then write(iulog,*)'FIND_TIMES: Failed to find dates bracketing desired time =', time - write(iulog,*) 'filename = '//trim(file%curr_filename) + write(iulog,'(2a)') 'filename = ', trim(file%curr_filename) write(iulog,*)' datatimem = ',file%datatimem write(iulog,*)' datatimep = ',file%datatimep endif @@ -1247,8 +1204,7 @@ subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found deallocate( all_data_times, stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'find_times: failed to deallocate all_data_times array; error = ',astat - call endrun + call endrun('find_times: failed to deallocate all_data_times array; error = '//int2str(astat)) end if if ( .not. file%cyclical ) then @@ -1273,7 +1229,6 @@ end subroutine find_times !------------------------------------------------------------------------ !------------------------------------------------------------------------ subroutine read_next_trcdata( flds, file ) - implicit none type (trfile), intent(inout) :: file type (trfld),intent(inout) :: flds(:) @@ -1456,7 +1411,6 @@ subroutine read_2d_trc( fid, vid, loc_arr, strt, cnt, file, order ) use polar_avg, only: polar_average use horizontal_interpolate, only : xy_interp - implicit none type(file_desc_t), intent(in) :: fid type(var_desc_t), intent(in) :: vid integer, intent(in) :: strt(:), cnt(:), order(2) @@ -1476,19 +1430,16 @@ subroutine read_2d_trc( fid, vid, loc_arr, strt, cnt, file, order ) nullify(wrk2d_in) allocate( wrk2d(cnt(1),cnt(2)), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr - call endrun + call endrun('read_2d_trc: wrk2d allocation error = '//int2str(ierr)) end if if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlon .or. cnt(2)/=file%nlat) then allocate( wrk2d_in(file%nlon, file%nlat), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr - call endrun + call endrun('read_2d_trc: wrk2d_in allocation error = '//int2str(ierr)) end if end if - ierr = pio_get_var( fid, vid, strt, cnt, wrk2d ) if(associated(wrk2d_in)) then wrk2d_in = reshape( wrk2d(:,:),(/file%nlon,file%nlat/), order=order ) @@ -1574,7 +1525,6 @@ subroutine read_za_trc( fid, vid, loc_arr, strt, cnt, file, order ) use ppgrid, only : pcols, begchunk, endchunk use phys_grid, only : get_ncols_p, get_rlat_all_p - implicit none type(file_desc_t), intent(in) :: fid type(var_desc_t), intent(in) :: vid integer, intent(in) :: strt(:), cnt(:) @@ -1591,19 +1541,16 @@ subroutine read_za_trc( fid, vid, loc_arr, strt, cnt, file, order ) nullify(wrk2d_in) allocate( wrk2d(cnt(1),cnt(2)), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr - call endrun + call endrun('read_2d_trc: wrk2d allocation error = '//int2str(ierr)) end if if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlat .or. cnt(2)/=file%nlev) then allocate( wrk2d_in(file%nlat, file%nlev), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr - call endrun + call endrun('read_2d_trc: wrk2d_in allocation error = '//int2str(ierr)) end if end if - ierr = pio_get_var( fid, vid, strt, cnt, wrk2d ) if(associated(wrk2d_in)) then wrk2d_in = reshape( wrk2d(:,:),(/file%nlat,file%nlev/), order=order ) @@ -1707,8 +1654,6 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) use polar_avg, only : polar_average use horizontal_interpolate, only : xy_interp - implicit none - type(file_desc_t), intent(in) :: fid type(var_desc_t), intent(in) :: vid integer, intent(in) :: strt(:), cnt(:), order(3) @@ -1731,8 +1676,7 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) nullify(wrk3d_in) allocate(wrk3d(cnt(1),cnt(2),cnt(3)), stat=ierr) if( ierr /= 0 ) then - write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr - call endrun + call endrun('read_3d_trc: wrk3d allocation error = '//int2str(ierr)) end if ierr = pio_get_var( fid, vid, strt, cnt, wrk3d ) @@ -1741,8 +1685,7 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) cnt(1)/=file%nlon.or.cnt(2)/=file%nlat.or.cnt(3)/=file%nlev) then allocate(wrk3d_in(file%nlon,file%nlat,file%nlev),stat=ierr) if( ierr /= 0 ) then - write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr - call endrun + call endrun('read_3d_trc: wrk3d allocation error = '//int2str(ierr)) end if wrk3d_in = reshape( wrk3d(:,:,:),(/file%nlon,file%nlat,file%nlev/), order=order ) deallocate(wrk3d) @@ -1802,8 +1745,7 @@ subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) deallocate( wrk3d_in, stat=astat ) end if if( astat/= 0 ) then - write(iulog,*) 'read_3d_trc: failed to deallocate wrk3d array; error = ',astat - call endrun + call endrun('read_3d_trc: failed to deallocate wrk3d array; error = '//int2str(astat)) endif if(dycore_is('LR')) call polar_average(file%nlev, loc_arr) end subroutine read_3d_trc @@ -1815,8 +1757,6 @@ subroutine interpolate_trcdata( state, flds, file, pbuf2d ) use physics_types,only : physics_state use physconst, only : cday, rga - implicit none - type(physics_state), intent(in) :: state(begchunk:endchunk) type (trfld), intent(inout) :: flds(:) type (trfile), intent(inout) :: file @@ -1997,7 +1937,7 @@ end subroutine interpolate_trcdata !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine get_dimension( fid, dname, dsize, dimid, data ) - implicit none + type(file_desc_t), intent(inout) :: fid character(*), intent(in) :: dname integer, intent(out) :: dsize @@ -2024,14 +1964,12 @@ subroutine get_dimension( fid, dname, dsize, dimid, data ) if ( associated(data) ) then deallocate(data, stat=ierr) if( ierr /= 0 ) then - write(iulog,*) 'get_dimension: data deallocation error = ',ierr - call endrun('get_dimension: failed to deallocate data array') + call endrun('get_dimension: failed to deallocate data array, error = '//int2str(ierr)) end if endif allocate( data(dsize), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'get_dimension: data allocation error = ',ierr - call endrun('get_dimension: failed to allocate data array') + call endrun('get_dimension: failed to allocate data array, error = '//int2str(ierr)) end if ierr = pio_inq_varid( fid, dname, vid ) @@ -2050,8 +1988,6 @@ end subroutine get_dimension !----------------------------------------------------------------------- subroutine set_cycle_indices( fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr ) - implicit none - type(file_desc_t), intent(inout) :: fileid integer, intent(out) :: cyc_ndx_beg integer, intent(out) :: cyc_ndx_end @@ -2065,8 +2001,7 @@ subroutine set_cycle_indices( fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr ) allocate( dates(timesize), stat=astat ) if( astat/= 0 ) then - write(*,*) 'set_cycle_indices: failed to allocate dates array; error = ',astat - call endrun + call endrun('set_cycle_indices: failed to allocate dates array; error = '//int2str(astat)) end if ierr = pio_inq_varid( fileid, 'date', dateid ) @@ -2083,12 +2018,10 @@ subroutine set_cycle_indices( fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr ) enddo deallocate( dates, stat=astat ) if( astat/= 0 ) then - write(*,*) 'set_cycle_indices: failed to deallocate dates array; error = ',astat - call endrun + call endrun('set_cycle_indices: failed to deallocate dates array; error = '//int2str(astat)) end if if (cyc_ndx_beg < 0) then - write(*,*) 'set_cycle_indices: cycle year not found : ' , cyc_yr - call endrun('set_cycle_indices: cycle year not found') + call endrun('set_cycle_indices: cycle year not found : '//int2str(cyc_yr)) endif end subroutine set_cycle_indices @@ -2100,8 +2033,6 @@ subroutine open_trc_datafile( fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_ use ioFileMod, only: getfil use cam_pio_utils, only: cam_pio_openfile - implicit none - character(*), intent(in) :: fname character(*), intent(in) :: path type(file_desc_t), intent(inout) :: piofile @@ -2129,32 +2060,28 @@ subroutine open_trc_datafile( fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_ ! call getfil( filepath, filen, 0 ) call cam_pio_openfile( piofile, filen, PIO_NOWRITE) - if(masterproc) write(iulog,*)'open_trc_datafile: ',trim(filen) + if(masterproc) write(iulog,'(2a)')'open_trc_datafile: ', trim(filen) call get_dimension(piofile, 'time', timesize) if ( associated(times) ) then deallocate(times, stat=ierr) if( ierr /= 0 ) then - write(iulog,*) 'open_trc_datafile: data deallocation error = ',ierr - call endrun('open_trc_datafile: failed to deallocate data array') + call endrun('open_trc_datafile: failed to deallocate data array, error = '//int2str(ierr)) end if endif allocate( times(timesize), stat=ierr ) if( ierr /= 0 ) then - write(iulog,*) 'open_trc_datafile: data allocation error = ',ierr - call endrun('open_trc_datafile: failed to allocate data array') + call endrun('open_trc_datafile: failed to allocate data array, error = '//int2str(ierr)) end if allocate( dates(timesize), stat=astat ) if( astat/= 0 ) then - if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate dates array; error = ',astat - call endrun + call endrun('open_trc_datafile: failed to allocate dates array; error = '//int2str(astat)) end if allocate( datesecs(timesize), stat=astat ) if( astat/= 0 ) then - if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate datesec array; error = ',astat - call endrun + call endrun('open_trc_datafile: failed to allocate datesec array; error = '//int2str(astat)) end if ierr = pio_inq_varid( piofile, 'date', dateid ) @@ -2191,18 +2118,15 @@ subroutine open_trc_datafile( fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_ deallocate( dates, stat=astat ) if( astat/= 0 ) then - if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate dates array; error = ',astat - call endrun + call endrun('open_trc_datafile: failed to deallocate dates array; error = '//int2str(astat)) end if deallocate( datesecs, stat=astat ) if( astat/= 0 ) then - if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate datesec array; error = ',astat - call endrun + call endrun('open_trc_datafile: failed to deallocate datesec array; error = '//int2str(astat)) end if if ( present(cyc_yr) .and. present(cyc_ndx_beg) ) then if (cyc_ndx_beg < 0) then - write(iulog,*) 'open_trc_datafile: cycle year not found : ' , cyc_yr call endrun('open_trc_datafile: cycle year not found '//trim(filepath)) endif endif @@ -2213,8 +2137,6 @@ end subroutine open_trc_datafile !-------------------------------------------------------------------------- subroutine specify_fields( specifier, fields ) - implicit none - character(len=*), intent(in) :: specifier(:) type(trfld), pointer, dimension(:) :: fields @@ -2228,8 +2150,7 @@ subroutine specify_fields( specifier, fields ) allocate(fld_name(nflds), src_name(nflds), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'specify_fields: failed to allocate fld_name, src_name arrays; error = ',astat - call endrun + call endrun('specify_fields: failed to allocate fld_name, src_name arrays; error = '//int2str(astat)) end if fld_cnt = 0 @@ -2266,8 +2187,7 @@ subroutine specify_fields( specifier, fields ) !----------------------------------------------------------------------- allocate( fields(fld_cnt), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'specify_fields: failed to allocate fields array; error = ',astat - call endrun + call endrun('specify_fields: failed to allocate fields array; error = '//int2str(astat)) end if do i = 1,fld_cnt @@ -2283,7 +2203,6 @@ end subroutine specify_fields subroutine init_trc_restart( whence, piofile, tr_file ) - implicit none character(len=*), intent(in) :: whence type(file_desc_t), intent(inout) :: piofile type(trfile), intent(inout) :: tr_file @@ -2327,8 +2246,6 @@ end subroutine init_trc_restart !------------------------------------------------------------------------- subroutine write_trc_restart( piofile, tr_file ) - implicit none - type(file_desc_t), intent(inout) :: piofile type(trfile), intent(inout) :: tr_file @@ -2350,8 +2267,6 @@ end subroutine write_trc_restart !------------------------------------------------------------------------- subroutine read_trc_restart( whence, piofile, tr_file ) - implicit none - character(len=*), intent(in) :: whence type(file_desc_t), intent(inout) :: piofile type(trfile), intent(inout) :: tr_file @@ -2388,8 +2303,6 @@ end subroutine read_trc_restart !------------------------------------------------------------------------------ subroutine interpz_conserve( nsrc, ntrg, src_x, trg_x, src, trg) - implicit none - integer, intent(in) :: nsrc ! dimension source array integer, intent(in) :: ntrg ! dimension target array real(r8), intent(in) :: src_x(nsrc+1) ! source coordinates @@ -2457,8 +2370,6 @@ end subroutine interpz_conserve !------------------------------------------------------------------------------ subroutine vert_interp_mixrat( ncol, nsrc, ntrg, trg_x, src, trg, p0, ps, hyai, hybi, use_flight_distance) - implicit none - integer, intent(in) :: ncol integer, intent(in) :: nsrc ! dimension source array integer, intent(in) :: ntrg ! dimension target array @@ -2597,7 +2508,6 @@ subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout ) ! ! Interpolate data from current time-interpolated values to model levels !-------------------------------------------------------------------------- - implicit none ! Arguments ! integer, intent(in) :: ncol ! number of atmospheric columns @@ -2674,7 +2584,6 @@ subroutine vert_interp_ub( ncol, nlevs, plevs, datain, dataout ) ! Interpolate data from current time-interpolated values to top interface pressure ! -- from mo_tgcm_ubc.F90 !-------------------------------------------------------------------------- - implicit none ! Arguments ! integer, intent(in) :: ncol @@ -2783,8 +2692,6 @@ subroutine advance_file(file) use shr_sys_mod, only: shr_sys_system use ioFileMod, only: getfil - implicit none - type(trfile), intent(inout) :: file !----------------------------------------------------------------------- @@ -2804,10 +2711,10 @@ subroutine advance_file(file) !----------------------------------------------------------------------- if( file%remove_trc_file ) then call getfil( file%curr_filename, loc_fname, 0 ) - write(iulog,*) 'advance_file: removing file = ',trim(loc_fname) + write(iulog,'(2a)') 'advance_file: removing file = ',trim(loc_fname) ctmp = 'rm -f ' // trim(loc_fname) - write(iulog,*) 'advance_file: fsystem issuing command - ' - write(iulog,*) trim(ctmp) + write(iulog,'(a)') 'advance_file: fsystem issuing command - ' + write(iulog,'(a)') trim(ctmp) call shr_sys_system( ctmp, istat ) end if @@ -2822,13 +2729,11 @@ subroutine advance_file(file) !----------------------------------------------------------------------- deallocate( file%curr_data_times, stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'advance_file: failed to deallocate file%curr_data_times array; error = ',astat - call endrun + call endrun('advance_file: failed to deallocate file%curr_data_times array; error = '//int2str(astat)) end if allocate( file%curr_data_times( size( file%next_data_times ) ), stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'advance_file: failed to allocate file%curr_data_times array; error = ',astat - call endrun + call endrun('advance_file: failed to allocate file%curr_data_times array; error = '//int2str(astat)) end if file%curr_data_times(:) = file%next_data_times(:) @@ -2839,8 +2744,7 @@ subroutine advance_file(file) deallocate( file%next_data_times, stat=astat ) if( astat/= 0 ) then - write(iulog,*) 'advance_file: failed to deallocate file%next_data_times array; error = ',astat - call endrun + call endrun('advance_file: failed to deallocate file%next_data_times array; error = '//int2str(astat)) end if nullify( file%next_data_times ) diff --git a/src/control/cam_esmf_mod.F90 b/src/control/cam_esmf_mod.F90 new file mode 100644 index 0000000000..f1312f2894 --- /dev/null +++ b/src/control/cam_esmf_mod.F90 @@ -0,0 +1,165 @@ +module cam_esmf_mod + + use shr_kind_mod , only : r8=>shr_kind_r8 + use ESMF , only : ESMF_Mesh, ESMF_Clock + use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_VMGetCurrent + use ESMF , only : ESMF_SUCCESS, ESMF_REDUCE_SUM + use nuopc_shr_methods , only : chkerr + use error_messages , only : alloc_err + use cam_logfile , only : iulog + + implicit none + private + + public :: cam_esmf_set_clock + public :: cam_esmf_set_mesh + public :: cam_esmf_set_areas + public :: cam_esmf_global_sum + + type(ESMF_Mesh) , public, protected :: model_mesh ! model mesh + type(ESMF_Clock), public, protected :: model_clock ! model clock + + real(r8), allocatable, public, protected :: model_areas(:) + real(r8), allocatable, public, protected :: mesh_areas(:) + + logical :: check_global_areas = .false. + + character(len=*), parameter :: u_FILE_u = __FILE__ + +!===================================================================== +contains +!===================================================================== + + subroutine cam_esmf_set_clock(clock_in, rc) + use ESMF, only: ESMF_Clock, ESMF_ClockCreate, ESMF_ClockIsCreated + use cam_abortutils, only: endrun + + ! Arguments + type(ESMF_Clock), intent(in) :: clock_in + integer , intent(out) :: rc + !--------------------------------------- + + rc = ESMF_SUCCESS + + if (ESMF_ClockIsCreated(model_clock, rc=rc)) then + call endrun('cam_esmf_set_clock: model_clock already set') + else + model_clock = ESMF_ClockCreate(clock_in, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + model_clock = clock_in + end if + + end subroutine cam_esmf_set_clock + + !===================================================================== + subroutine cam_esmf_set_mesh(mesh_in) + type(ESMF_Mesh) , intent(in) :: mesh_in + + model_mesh = mesh_in + + end subroutine cam_esmf_set_mesh + + !===================================================================== + subroutine cam_esmf_set_areas(model_areas_in, mesh_areas_in, rc) + use spmd_utils, only: masterproc + + ! Arguments + real(r8), intent(in) :: model_areas_in(:) + real(r8), intent(in) :: mesh_areas_in(:) + integer , intent(out) :: rc + + ! Local variables + type(ESMF_VM) :: vm + integer :: locsize + integer :: ng + integer :: istat + real(r8) :: global_mesh_area(1) + real(r8) :: global_model_area(1) + real(r8) :: local_mesh_area(1) + real(r8) :: local_model_area(1) + !--------------------------------------- + + rc = ESMF_SUCCESS + + locsize = size(model_areas_in) + + allocate(model_areas(locsize), stat=istat) + call alloc_err(istat,'cam_esmf_set_areas','model_areas',locsize) + model_areas(:) = model_areas_in(:) + + allocate(mesh_areas(locsize), stat=istat) + call alloc_err(istat,'cam_esmf_set_areas','mesh_areas',locsize) + mesh_areas(:) = mesh_areas_in(:) + + ! Compare global sum of model_areas and mesh_areas + local_model_area(1) = 0._r8 + local_mesh_area(1) = 0._r8 + do ng = 1,locsize + local_model_area(1) = local_model_area(1) + model_areas(ng) + local_mesh_area(1) = local_mesh_area(1) + mesh_areas(ng) + end do + + if (check_global_areas) then + call ESMF_VMGetCurrent(vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_VMAllreduce(vm, senddata=local_model_area, recvdata=global_model_area, & + count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call ESMF_VMAllreduce(vm, senddata=local_mesh_area, recvdata=global_mesh_area, & + count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + if (masterproc) then + write(iulog,'(a,d13.5)') ' global mesh area = ',global_mesh_area(1) + write(iulog,'(a,d13.5)') ' global model area = ',global_model_area(1) + end if + end if + + end subroutine cam_esmf_set_areas + + !===================================================================== + subroutine cam_esmf_global_sum(fldname, flddata, global_sum_model, global_sum_mesh, rc) + + ! Arguments + character(len=*), intent(in) :: fldname + real(r8), intent(in) :: flddata(:) + real(r8), intent(out) :: global_sum_model + real(r8), intent(out) :: global_sum_mesh + integer , intent(out) :: rc + + ! local variables + type(ESMF_VM) :: vm + integer :: ng + real(r8) :: local_sum_model(1) + real(r8) :: local_sum_mesh(1) + real(r8) :: global_model(1) + real(r8) :: global_mesh(1) + !--------------------------------------- + + rc = ESMF_SUCCESS + + local_sum_model(1) = 0._r8 + local_sum_mesh(1) = 0._r8 + do ng=1,size(flddata) + local_sum_model(1) = local_sum_model(1) + flddata(ng) * model_areas(ng) + local_sum_mesh(1) = local_sum_mesh(1) + flddata(ng) * mesh_areas(ng) + end do + + call ESMF_VMGetCurrent(vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_sum_model, recvdata=global_model, & + count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=local_sum_mesh, recvdata=global_mesh, & + count=1, reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + global_sum_model = global_model(1) + global_sum_mesh = global_mesh(1) + + end subroutine cam_esmf_global_sum + +end module cam_esmf_mod diff --git a/src/cpl/nuopc/atm_comp_nuopc.F90 b/src/cpl/nuopc/atm_comp_nuopc.F90 index 10150e4dc2..b75e187e45 100644 --- a/src/cpl/nuopc/atm_comp_nuopc.F90 +++ b/src/cpl/nuopc/atm_comp_nuopc.F90 @@ -69,6 +69,7 @@ module atm_comp_nuopc use pio , only : pio_noerr, pio_bcast_error, pio_internal_error, pio_seterrorhandling use pio , only : pio_def_var, pio_get_var, pio_put_var, PIO_INT use ioFileMod + use cam_esmf_mod , only : cam_esmf_set_clock, cam_esmf_set_mesh !$use omp_lib , only : omp_set_num_threads implicit none @@ -129,9 +130,6 @@ module atm_comp_nuopc real(R8) , parameter :: grid_tol = 1.e-2_r8 ! tolerance for calculated lat/lon vs read in - type(ESMF_Mesh) :: model_mesh ! model_mesh - type(ESMF_Clock) :: model_clock ! model_clock - !=============================================================================== contains !=============================================================================== @@ -335,6 +333,8 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) integer, intent(out) :: rc ! local variables + type(ESMF_Mesh) :: model_mesh + type(ESMF_Clock) :: model_clock type(ESMF_VM) :: vm type(ESMF_Time) :: currTime ! Current time type(ESMF_Time) :: startTime ! Start time @@ -623,6 +623,10 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call shr_sys_abort( subname//'ERROR:: bad calendar for ESMF' ) end if + ! Create model_clock as a variable in cam_esmf_mod.F90 - needed for generating streams + call cam_esmf_set_clock(clock, rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Initialize module orbital values and update orbital call cam_orbital_init(gcomp, iulog, masterproc, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -761,11 +765,11 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) call realize_fields(gcomp, model_mesh, flds_scalar_name, flds_scalar_num, single_column, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Create model_clock as a module variable - needed for generating streams - model_clock = clock + ! Set module variables in src/control/cam_esmf_mod.F90 (must be done before call to export_fields) + call cam_esmf_set_mesh(model_mesh) ! Create cam export array and set the state scalars - call export_fields( gcomp, model_mesh, model_clock, cam_out, rc=rc ) + call export_fields( gcomp, cam_out, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return call get_horiz_grid_dim_d(hdim1_d, hdim2_d) @@ -920,7 +924,7 @@ subroutine DataInitialize(gcomp, rc) call import_fields( gcomp, cam_in, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return call cam_run1 ( cam_in, cam_out ) - call export_fields( gcomp, model_mesh, model_clock, cam_out, rc=rc ) + call export_fields( gcomp, cam_out, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return else call cam_read_srfrest( gcomp, clock, rc=rc ) @@ -928,7 +932,7 @@ subroutine DataInitialize(gcomp, rc) call import_fields( gcomp, cam_in, restart_init=.true., rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return call cam_run1 ( cam_in, cam_out ) - call export_fields( gcomp, model_mesh, model_clock, cam_out, rc=rc ) + call export_fields( gcomp, cam_out, rc=rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if @@ -1175,7 +1179,7 @@ subroutine ModelAdvance(gcomp, rc) if (mediator_present) then ! Set export fields call t_startf ('CAM_export') - call export_fields( gcomp, model_mesh, model_clock, cam_out, rc ) + call export_fields( gcomp, cam_out, rc ) if (ChkErr(rc,__LINE__,u_FILE_u)) return call t_stopf ('CAM_export') diff --git a/src/cpl/nuopc/atm_import_export.F90 b/src/cpl/nuopc/atm_import_export.F90 index 907c4d2f1e..59ab67d4b0 100644 --- a/src/cpl/nuopc/atm_import_export.F90 +++ b/src/cpl/nuopc/atm_import_export.F90 @@ -13,7 +13,7 @@ module atm_import_export use shr_mpi_mod , only : shr_mpi_min, shr_mpi_max use nuopc_shr_methods , only : chkerr use cam_logfile , only : iulog - use cam_history , only: outfld + use cam_history , only : outfld use spmd_utils , only : masterproc, mpicom use srf_field_check , only : set_active_Sl_ram1 use srf_field_check , only : set_active_Sl_fv @@ -27,6 +27,7 @@ module atm_import_export use atm_stream_ndep , only : ndep_stream_active use chemistry , only : chem_has_ndep_flx use cam_control_mod , only : aqua_planet, simple_phys + use cam_esmf_mod , only : cam_esmf_set_areas implicit none private ! except @@ -497,6 +498,9 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, singl end do deallocate(area) + call cam_esmf_set_areas(model_areas, mesh_areas, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + ! Determine flux correction factors (module variables) do n = 1,numOwnedElements mod2med_areacor(n) = model_areas(n) / mesh_areas(n) @@ -540,9 +544,9 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) use phys_grid , only : get_ncols_p use ppgrid , only : begchunk, endchunk use shr_const_mod , only : shr_const_stebol - use co2_cycle , only : c_i, co2_readFlux_ocn, co2_readFlux_fuel - use co2_cycle , only : co2_transport, co2_time_interp_ocn, co2_time_interp_fuel - use co2_cycle , only : data_flux_ocn, data_flux_fuel + use co2_cycle , only : c_i, co2_readFlux_fuel + use co2_cycle , only : co2_transport + use co2_data_flux , only : data_flux_fuel, co2_data_flux_advance use physconst , only : mwco2 use time_manager , only : is_first_step, get_nstep @@ -882,11 +886,6 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) g = g + 1 end do end do - else - ! Consistency check - if (co2_readFlux_ocn) then - call shr_sys_abort(subname // ':: co2_readFlux_ocn and x2a_Faoo_fco2_ocn cannot both be active') - end if end if call state_getfldptr(importState, 'Faoo_fdms_ocn', fldptr=fldptr1d, exists=exists, rc=rc) @@ -950,16 +949,13 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) if (co2_transport() .and. overwrite_flds) then ! Interpolate in time for flux data read in - if (co2_readFlux_ocn) then - call co2_time_interp_ocn - end if if (co2_readFlux_fuel) then - call co2_time_interp_fuel + call co2_data_flux_advance() end if - ! from ocn : data read in or from coupler or zero ! from fuel: data read in or zero - ! from lnd : through coupler or zero + ! from ocn : from mediator or zero + ! from lnd : from mediator or zero ! all co2 fluxes in unit kgCO2/m2/s do c=begchunk,endchunk @@ -968,10 +964,6 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) ! co2 flux from ocn if (exists_fco2_ocn) then cam_in(c)%cflx(i,c_i(1)) = cam_in(c)%fco2_ocn(i) - else if (co2_readFlux_ocn) then - ! convert from molesCO2/m2/s to kgCO2/m2/s - cam_in(c)%cflx(i,c_i(1)) = & - -data_flux_ocn%co2flx(i,c)*(1._r8- cam_in(c)%landfrac(i))*mwco2*1.0e-3_r8 else cam_in(c)%cflx(i,c_i(1)) = 0._r8 end if @@ -1012,7 +1004,7 @@ end subroutine import_fields !=============================================================================== - subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) + subroutine export_fields(gcomp, cam_out, rc) ! ----------------------------------------------------- ! Set field pointers in export set @@ -1031,15 +1023,12 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) ! input/output variables type(ESMF_GridComp) :: gcomp - type(ESMF_Mesh) , intent(in) :: model_mesh - type(ESMF_Clock), intent(in) :: model_clock type(cam_out_t) , intent(inout) :: cam_out(begchunk:endchunk) integer , intent(out) :: rc ! local variables type(ESMF_State) :: exportState type(ESMF_State) :: importState - type(ESMF_Clock) :: clock integer :: i,m,c,n,g ! indices integer :: nstep logical :: exists @@ -1256,7 +1245,7 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) ! The ndep_stream_nl namelist group is read in stream_ndep_init. This sets whether ! or not the stream will be used. if (.not. stream_ndep_is_initialized) then - call stream_ndep_init(model_mesh, model_clock, rc) + call stream_ndep_init(rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return stream_ndep_is_initialized = .true. end if diff --git a/src/cpl/nuopc/atm_stream_ndep.F90 b/src/cpl/nuopc/atm_stream_ndep.F90 index f54509b269..b9ec64a5e9 100644 --- a/src/cpl/nuopc/atm_stream_ndep.F90 +++ b/src/cpl/nuopc/atm_stream_ndep.F90 @@ -17,6 +17,7 @@ module atm_stream_ndep use spmd_utils , only : mpi_character, mpi_integer use cam_logfile , only : iulog use cam_abortutils , only : endrun + use cam_esmf_mod , only : model_clock, model_mesh implicit none private @@ -34,7 +35,6 @@ module atm_stream_ndep type(shr_strdata_type) :: sdat_ndep ! input data stream logical, public :: stream_ndep_is_initialized = .false. character(len=CS) :: stream_varlist_ndep(2) - type(ESMF_Clock) :: model_clock character(len=*), parameter :: sourcefile = __FILE__ @@ -132,23 +132,22 @@ subroutine stream_ndep_readnl(nlfile) end subroutine stream_ndep_readnl - subroutine stream_ndep_init(model_mesh, model_clock, rc) + !================================================================ + subroutine stream_ndep_init(rc) use dshr_strdata_mod, only: shr_strdata_init_from_inline ! input/output variables - type(ESMF_CLock), intent(in) :: model_clock - type(ESMF_Mesh) , intent(in) :: model_mesh - integer , intent(out) :: rc + integer, intent(out) :: rc ! local variables character(*), parameter :: subName = "('stream_ndep_init')" + !----------------------------------------------------------------------- rc = ESMF_SUCCESS if (.not.ndep_stream_active) then return end if - ! - ! Initialize data stream information. + ! Read in units call stream_ndep_check_units(stream_ndep_data_filename) @@ -250,7 +249,6 @@ subroutine stream_ndep_interp(cam_out, rc) ! NDEP read from forcing is expected to be in units of gN/m2/sec - but the mediator ! expects units of kgN/m2/sec real(r8), parameter :: scale_ndep = .001_r8 - !----------------------------------------------------------------------- ! Advance sdat stream diff --git a/src/physics/cam/co2_cycle.F90 b/src/physics/cam/co2_cycle.F90 index d3b4ee5e7e..cffbc5618d 100644 --- a/src/physics/cam/co2_cycle.F90 +++ b/src/physics/cam/co2_cycle.F90 @@ -1,22 +1,20 @@ module co2_cycle -!------------------------------------------------------------------------------- -! -! Purpose: -! Provides distributions of CO2_LND, CO2_OCN, CO2_FF, CO2 -! Surface flux from CO2_LND and CO2_OCN can be provided by the flux coupler. -! Surface flux from CO2_FFF and CO2_OCN can be read from a file. -! -! Author: Jeff Lee, Keith Lindsay -! -!------------------------------------------------------------------------------- - - use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl - use co2_data_flux, only: co2_data_flux_type - use srf_field_check, only: active_Faoo_fco2_ocn + !------------------------------------------------------------------------------- + ! + ! Purpose: + ! Provides distributions of CO2_LND, CO2_OCN, CO2_FF, CO2 + ! Surface flux from CO2_LND and CO2_OCN provided by the mediator. + ! Surface flux from CO2_FFF can be read from a file. + ! + ! Author: Jeff Lee, Keith Lindsay + ! Mariana Vertenstein, Refactored for NUOPC Stream functionality + ! + !------------------------------------------------------------------------------- - implicit none + use shr_kind_mod, only: r8=>shr_kind_r8 + implicit none private ! Public interfaces @@ -26,402 +24,312 @@ module co2_cycle public co2_implements_cnst ! returns true if consituent is implemented by this package public co2_init_cnst ! initialize mixing ratios if not read from initial file public co2_init ! initialize (history) variables - public co2_time_interp_ocn ! time interpolate co2 flux - public co2_time_interp_fuel ! time interpolate co2 flux public co2_cycle_set_ptend ! set tendency from aircraft emissions - ! Public data - public data_flux_ocn ! data read in for co2 flux from ocn - public data_flux_fuel ! data read in for co2 flux from fuel - - type(co2_data_flux_type) :: data_flux_ocn - type(co2_data_flux_type) :: data_flux_fuel - - public c_i ! global index for new constituents - public co2_readFlux_ocn ! read ocn co2 flux from data file - public co2_readFlux_fuel ! read fuel co2 flux from data file - ! Namelist variables - logical :: co2_flag = .false. ! true => turn on co2 code, namelist variable - logical :: co2_readFlux_ocn = .false. ! true => read ocn co2 flux from date file, namelist variable - logical :: co2_readFlux_fuel = .false. ! true => read fuel co2 flux from date file, namelist variable - logical :: co2_readFlux_aircraft = .false. ! true => read aircraft co2 flux from date file, namelist variable - character(len=cl) :: co2flux_ocn_file = 'unset' ! co2 flux from ocn - character(len=cl) :: co2flux_fuel_file = 'unset' ! co2 flux from fossil fuel + logical :: co2_flag = .false. ! true => turn on co2 code, namelist variable + logical, public, protected :: co2_readFlux_fuel = .false. ! true => read fuel co2 flux from date file, namelist variable + logical, public, protected :: co2_readFlux_aircraft = .false. ! true => read aircraft co2 flux from date file, namelist variable !------------------------------------------------------------------------------- ! new constituents !------------------------------------------------------------------------------- - integer, parameter :: ncnst=4 ! number of constituents implemented + integer, parameter :: ncnst=4 ! number of constituents implemented + integer, public, protected :: c_i(ncnst) ! global index for new constituents character(len=7), dimension(ncnst), parameter :: & ! constituent names - c_names = (/'CO2_OCN', 'CO2_FFF', 'CO2_LND', 'CO2 '/) - - integer :: co2_ocn_glo_ind ! global index of 'CO2_OCN' - integer :: co2_fff_glo_ind ! global index of 'CO2_FFF' - integer :: co2_lnd_glo_ind ! global index of 'CO2_LND' - integer :: co2_glo_ind ! global index of 'CO2' + c_names = (/'CO2_OCN', 'CO2_FFF', 'CO2_LND', 'CO2 '/) - integer, dimension(ncnst) :: c_i ! global index + integer :: co2_fff_glo_ind = -1 ! global index of 'CO2_FFF' + integer :: co2_glo_ind = -1 ! global index of 'CO2' + integer :: idx_ac_CO2 = -1 ! pbuf index of aircraft CO2 field !=============================================================================== contains !=============================================================================== -subroutine co2_cycle_readnl(nlfile) - -!------------------------------------------------------------------------------- -! Purpose: Read co2_cycle_nl namelist group. -!------------------------------------------------------------------------------- - - use namelist_utils, only: find_group_name - use units, only: getunit, freeunit - use spmd_utils, only: masterproc - use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_logical, mpi_character - use cam_logfile, only: iulog - use cam_abortutils, only: endrun - - ! Arguments - character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - - ! Local variables - integer :: unitn, ierr - character(len=256) :: msg - character(len=*), parameter :: subname = 'co2_cycle_readnl' - - namelist /co2_cycle_nl/ co2_flag, co2_readFlux_ocn, co2_readFlux_fuel, co2_readFlux_aircraft, & - co2flux_ocn_file, co2flux_fuel_file - !---------------------------------------------------------------------------- - - if (masterproc) then - unitn = getunit() - open( unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'co2_cycle_nl', status=ierr) - if (ierr == 0) then - read(unitn, co2_cycle_nl, iostat=ierr) - if (ierr /= 0) then - call endrun(subname // ':: ERROR reading namelist') + subroutine co2_cycle_readnl(nlfile) + + !-------------------------------------------- + ! Purpose: Read co2_cycle_nl namelist group. + !-------------------------------------------- + + use namelist_utils, only: find_group_name + use spmd_utils, only: masterproc, mpicom, masterprocid + use spmd_utils, only: mpi_logical, mpi_character, mpi_integer + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use co2_data_flux, only: co2_data_flux_readnl + + ! Arguments + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! Local variables + integer :: unitn, ierr + character(len=256) :: msg + character(len=*), parameter :: subname = 'co2_cycle_readnl' + + namelist /co2_cycle_nl/ & + co2_flag, & + co2_readFlux_aircraft, & ! if true, read aircraft data + co2_readFlux_fuel ! if true, read fuel data + !---------------------------------------------------------------------------- + + if (masterproc) then + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'co2_cycle_nl', status=ierr) + if (ierr == 0) then + read(unitn, co2_cycle_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading co2_cycle_nl namelist') + end if end if + close(unitn) end if - close(unitn) - call freeunit(unitn) - end if - - ! Broadcast namelist variables - call mpi_bcast(co2_flag, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_flag") - call mpi_bcast(co2_readFlux_ocn, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_readFlux_ocn") - call mpi_bcast(co2_readFlux_fuel, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_readFlux_fuel") - call mpi_bcast(co2_readFlux_aircraft, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_readFlux_aircraft") - call mpi_bcast(co2flux_ocn_file, len(co2flux_ocn_file), mpi_character, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_ocn_file") - call mpi_bcast(co2flux_fuel_file, len(co2flux_fuel_file), mpi_character, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_file") - - ! Consistency check - if (co2_readFlux_ocn .and. active_Faoo_fco2_ocn) then - msg = subname//': ERROR: reading ocn flux dataset is enabled, but coupler is setting'& - //' the ocn co2 flux. Cannot do both.' - write(iulog,*) trim(msg) - call endrun(trim(msg)) - end if - -end subroutine co2_cycle_readnl -!=============================================================================== + call mpi_bcast(co2_flag, 1, mpi_logical, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_flag") + call mpi_bcast(co2_readFlux_aircraft, 1, mpi_logical, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_readFlux_aircraft") + call mpi_bcast(co2_readFlux_fuel, 1, mpi_logical, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2_readFlux_fuel") + + if (masterproc) then + write(iulog, '(a, l4)') "co2_flag = ", co2_flag + write(iulog, '(a, l4)')"co2_readFlux_aircraft = ", co2_readFlux_aircraft + write(iulog, '(a, l4)')"co2_readFlux_fuel = ", co2_readFlux_fuel + end if -subroutine co2_register + if (co2_readFlux_fuel) then + call co2_data_flux_readnl(nlfile) + end if -!------------------------------------------------------------------------------- -! Purpose: register advected constituents -!------------------------------------------------------------------------------- + end subroutine co2_cycle_readnl - use physconst, only: mwco2, cpair - use constituents, only: cnst_add + !=============================================================================== - ! Local variables - real(r8), dimension(ncnst) :: & - c_mw, &! molecular weights - c_cp, &! heat capacities - c_qmin ! minimum mmr + subroutine co2_register - integer :: i + !------------------------------------------------------------------------------- + ! Purpose: register advected constituents + !------------------------------------------------------------------------------- - !---------------------------------------------------------------------------- + use physconst, only: mwco2, cpair + use constituents, only: cnst_add - if (.not. co2_flag) return + ! Local variables + real(r8), dimension(ncnst) :: & + c_mw, &! molecular weights + c_cp, &! heat capacities + c_qmin ! minimum mmr - c_mw = (/ mwco2, mwco2, mwco2, mwco2 /) - c_cp = (/ cpair, cpair, cpair, cpair /) - c_qmin = (/ 1.e-20_r8, 1.e-20_r8, 1.e-20_r8, 1.e-20_r8 /) + integer :: icnst + !---------------------------------------------------------------------------- - ! register CO2 constiuents as dry tracers, set indices + if (.not. co2_flag) return - do i = 1, ncnst - call cnst_add(c_names(i), c_mw(i), c_cp(i), c_qmin(i), c_i(i), longname=c_names(i), mixtype='dry') + c_mw = (/ mwco2, mwco2, mwco2, mwco2 /) + c_cp = (/ cpair, cpair, cpair, cpair /) + c_qmin = (/ 1.e-20_r8, 1.e-20_r8, 1.e-20_r8, 1.e-20_r8 /) - select case (trim(c_names(i))) - case ('CO2_OCN') - co2_ocn_glo_ind = c_i(i) - case ('CO2_FFF') - co2_fff_glo_ind = c_i(i) - case ('CO2_LND') - co2_lnd_glo_ind = c_i(i) - case ('CO2') - co2_glo_ind = c_i(i) - end select - end do + ! register CO2 constiuents as dry tracers, set indices -end subroutine co2_register + do icnst = 1, ncnst + call cnst_add(c_names(icnst), c_mw(icnst), c_cp(icnst), c_qmin(icnst), c_i(icnst), & + longname=c_names(icnst), mixtype='dry') -!=============================================================================== + select case (trim(c_names(icnst))) + case ('CO2_FFF') + co2_fff_glo_ind = c_i(icnst) + case ('CO2') + co2_glo_ind = c_i(icnst) + end select + end do -function co2_transport() + end subroutine co2_register -!------------------------------------------------------------------------------- -! Purpose: return true if this package is active -!------------------------------------------------------------------------------- + !=============================================================================== - ! Return value - logical :: co2_transport + logical function co2_transport() - !---------------------------------------------------------------------------- + !------------------------------------------------------------------------------- + ! Purpose: return true if this package is active + !------------------------------------------------------------------------------- - co2_transport = co2_flag + !---------------------------------------------------------------------------- -end function co2_transport + co2_transport = co2_flag -!=============================================================================== + end function co2_transport -function co2_implements_cnst(name) + !=============================================================================== -!------------------------------------------------------------------------------- -! Purpose: return true if specified constituent is implemented by this package -!------------------------------------------------------------------------------- + logical function co2_implements_cnst(name) - ! Return value - logical :: co2_implements_cnst + !------------------------------------------------------------------------------- + ! Purpose: return true if specified constituent is implemented by this package + !------------------------------------------------------------------------------- - ! Arguments - character(len=*), intent(in) :: name ! constituent name + ! Arguments + character(len=*), intent(in) :: name ! constituent name - ! Local variables - integer :: m + ! Local variables + integer :: mind + !---------------------------------------------------------------------------- - !---------------------------------------------------------------------------- + co2_implements_cnst = .false. - co2_implements_cnst = .false. + if (.not. co2_flag) return - if (.not. co2_flag) return + do mind = 1, ncnst + if (name == c_names(mind)) then + co2_implements_cnst = .true. + return + end if + end do - do m = 1, ncnst - if (name == c_names(m)) then - co2_implements_cnst = .true. + end function co2_implements_cnst + + !=============================================================================== + + subroutine co2_init_cnst(name, latvals, lonvals, mask, q) + + !------------------------------------------------------------------------------- + ! Purpose: + ! Set initial values of CO2_OCN, CO2_FFF, CO2_LND, CO2 + ! Need to be called from process_inidat in inidat.F90 + ! (or, initialize co2 in co2_timestep_init) + !------------------------------------------------------------------------------- + + use chem_surfvals, only: chem_surfvals_get + + ! Arguments + character(len=*), intent(in) :: name ! constituent name + real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol) + real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol) + logical, intent(in) :: mask(:) ! Only initialize where .true. + real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev) + + ! Local variables + integer :: kindx + !---------------------------------------------------------------------------- + + if (.not. co2_flag) return + + do kindx = 1, size(q, 2) + select case (name) + case ('CO2_OCN') + where(mask) + q(:, kindx) = chem_surfvals_get('CO2MMR') + end where + case ('CO2_FFF') + where(mask) + q(:, kindx) = chem_surfvals_get('CO2MMR') + end where + case ('CO2_LND') + where(mask) + q(:, kindx) = chem_surfvals_get('CO2MMR') + end where + case ('CO2') + where(mask) + q(:, kindx) = chem_surfvals_get('CO2MMR') + end where + end select + end do + + end subroutine co2_init_cnst + + !=============================================================================== + + subroutine co2_init + + !------------------------------------------------------------------------------- + ! Purpose: initialize co2, + ! declare history variables, + ! read co2 flux form fuel, as data_flux_fuel + !------------------------------------------------------------------------------- + + use cam_history, only: addfld, add_default, horiz_only + use constituents, only: cnst_name, cnst_longname, sflxnam + use physics_buffer, only: pbuf_get_index + + ! Local variables + integer :: m, mm + !---------------------------------------------------------------------------- + + if (.not. co2_flag) return + + ! Add constituents and fluxes to history file + do m = 1, ncnst + mm = c_i(m) + + call addfld(trim(cnst_name(mm))//'_BOT', horiz_only, 'A', 'kg/kg', trim(cnst_longname(mm))//', Bottom Layer') + call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm)) + call addfld(sflxnam(mm), horiz_only, 'A', 'kg/m2/s', trim(cnst_name(mm))//' surface flux') + + call add_default(cnst_name(mm), 1, ' ') + call add_default(sflxnam(mm), 1, ' ') + + ! The addfld call for the 'TM*' fields are made by default in the + ! constituent_burden module. + call add_default('TM'//trim(cnst_name(mm)), 1, ' ') + end do + + ! Find and store the aircraft CO2 index + idx_ac_CO2 = pbuf_get_index('ac_CO2') + + end subroutine co2_init + + !=============================================================================== + subroutine co2_cycle_set_ptend(state, pbuf, ptend) + + !------------------------------------------------------------------------------- + ! Purpose: + ! Set ptend, using aircraft CO2 emissions in ac_CO2 from pbuf + !------------------------------------------------------------------------------- + + use physics_types, only: physics_state, physics_ptend, physics_ptend_init + use physics_buffer, only: physics_buffer_desc, pbuf_get_field + use constituents, only: pcnst + use ppgrid, only: pver + use physconst, only: gravit + + ! Arguments + type(physics_state), intent(in) :: state + type(physics_buffer_desc), pointer :: pbuf(:) + type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies + + ! Local variables + logical :: lq(pcnst) + integer :: ifld, ncol, k + real(r8), pointer :: ac_CO2(:,:) + !---------------------------------------------------------------------------- + + if (.not. co2_flag .or. .not. co2_readFlux_aircraft) then + call physics_ptend_init(ptend, state%psetcols, 'none') return end if - end do - -end function co2_implements_cnst - -!=============================================================================== - -subroutine co2_init_cnst(name, latvals, lonvals, mask, q) - -!------------------------------------------------------------------------------- -! Purpose: -! Set initial values of CO2_OCN, CO2_FFF, CO2_LND, CO2 -! Need to be called from process_inidat in inidat.F90 -! (or, initialize co2 in co2_timestep_init) -!------------------------------------------------------------------------------- - - use chem_surfvals, only: chem_surfvals_get - - ! Arguments - character(len=*), intent(in) :: name ! constituent name - real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol) - real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol) - logical, intent(in) :: mask(:) ! Only initialize where .true. - real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev) - - ! Local variables - integer :: k - - !---------------------------------------------------------------------------- - - if (.not. co2_flag) return - - do k = 1, size(q, 2) - select case (name) - case ('CO2_OCN') - where(mask) - q(:, k) = chem_surfvals_get('CO2MMR') - end where - case ('CO2_FFF') - where(mask) - q(:, k) = chem_surfvals_get('CO2MMR') - end where - case ('CO2_LND') - where(mask) - q(:, k) = chem_surfvals_get('CO2MMR') - end where - case ('CO2') - where(mask) - q(:, k) = chem_surfvals_get('CO2MMR') - end where - end select - end do - -end subroutine co2_init_cnst - -!=============================================================================== - -subroutine co2_init - -!------------------------------------------------------------------------------- -! Purpose: initialize co2, -! declare history variables, -! read co2 flux form ocn, as data_flux_ocn -! read co2 flux form fule, as data_flux_fuel -!------------------------------------------------------------------------------- - - use cam_history, only: addfld, add_default, horiz_only - use co2_data_flux, only: co2_data_flux_init - use constituents, only: cnst_name, cnst_longname, sflxnam - - ! Local variables - integer :: m, mm - !---------------------------------------------------------------------------- + ! aircraft fluxes are added to 'CO2_FFF' and 'CO2' tendencies + lq(:) = .false. + lq(co2_fff_glo_ind) = .true. + lq(co2_glo_ind) = .true. - if (.not. co2_flag) return + call physics_ptend_init(ptend, state%psetcols, 'co2_cycle_ac', lq=lq) - ! Add constituents and fluxes to history file - do m = 1, ncnst - mm = c_i(m) - - call addfld(trim(cnst_name(mm))//'_BOT', horiz_only, 'A', 'kg/kg', trim(cnst_longname(mm))//', Bottom Layer') - call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm)) - call addfld(sflxnam(mm), horiz_only, 'A', 'kg/m2/s', trim(cnst_name(mm))//' surface flux') - - call add_default(cnst_name(mm), 1, ' ') - call add_default(sflxnam(mm), 1, ' ') - - ! The addfld call for the 'TM*' fields are made by default in the - ! constituent_burden module. - call add_default('TM'//trim(cnst_name(mm)), 1, ' ') - end do - - ! Read flux data - if (co2_readFlux_ocn) then - call co2_data_flux_init ( co2flux_ocn_file, 'CO2_flux', data_flux_ocn ) - end if - - if (co2_readFlux_fuel) then - call co2_data_flux_init ( co2flux_fuel_file, 'CO2_flux', data_flux_fuel ) - end if - -end subroutine co2_init - -!=============================================================================== + if (idx_ac_CO2 > 0) then + call pbuf_get_field(pbuf, idx_ac_CO2, ac_CO2) -subroutine co2_time_interp_ocn - -!------------------------------------------------------------------------------- -! Purpose: Time interpolate co2 flux to current time. -! Read in new monthly data if necessary -!------------------------------------------------------------------------------- - - use time_manager, only: is_first_step - use co2_data_flux, only: co2_data_flux_advance - - !---------------------------------------------------------------------------- - - if (.not. co2_flag) return - - if (co2_readFlux_ocn) then - call co2_data_flux_advance ( data_flux_ocn ) - endif - -end subroutine co2_time_interp_ocn - -!=============================================================================== - -subroutine co2_time_interp_fuel - -!------------------------------------------------------------------------------- -! Purpose: Time interpolate co2 flux to current time. -! Read in new monthly data if necessary -!------------------------------------------------------------------------------- - - use time_manager, only: is_first_step - use co2_data_flux, only: co2_data_flux_advance - - !---------------------------------------------------------------------------- - - if (.not. co2_flag) return - - if (co2_readFlux_fuel) then - call co2_data_flux_advance ( data_flux_fuel ) - endif - -end subroutine co2_time_interp_fuel - -!=============================================================================== - -subroutine co2_cycle_set_ptend(state, pbuf, ptend) - -!------------------------------------------------------------------------------- -! Purpose: -! Set ptend, using aircraft CO2 emissions in ac_CO2 from pbuf -!------------------------------------------------------------------------------- - - use physics_types, only: physics_state, physics_ptend, physics_ptend_init - use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field - use constituents, only: pcnst - use ppgrid, only: pver - use physconst, only: gravit - - ! Arguments - type(physics_state), intent(in) :: state - type(physics_buffer_desc), pointer :: pbuf(:) - type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies - - ! Local variables - logical :: lq(pcnst) - integer :: ifld, ncol, k - real(r8), pointer :: ac_CO2(:,:) - - !---------------------------------------------------------------------------- - - if (.not. co2_flag .or. .not. co2_readFlux_aircraft) then - call physics_ptend_init(ptend, state%psetcols, 'none') - return - end if - - ! aircraft fluxes are added to 'CO2_FFF' and 'CO2' tendencies - lq(:) = .false. - lq(co2_fff_glo_ind) = .true. - lq(co2_glo_ind) = .true. - - call physics_ptend_init(ptend, state%psetcols, 'co2_cycle_ac', lq=lq) - - ifld = pbuf_get_index('ac_CO2') - call pbuf_get_field(pbuf, ifld, ac_CO2) - - ! [ac_CO2] = 'kg m-2 s-1' - ! [ptend%q] = 'kg kg-1 s-1' - ncol = state%ncol - do k = 1, pver - ptend%q(:ncol,k,co2_fff_glo_ind) = gravit * state%rpdeldry(:ncol,k) * ac_CO2(:ncol,k) - ptend%q(:ncol,k,co2_glo_ind) = gravit * state%rpdeldry(:ncol,k) * ac_CO2(:ncol,k) - end do - -end subroutine co2_cycle_set_ptend + ! [ac_CO2] = 'kg m-2 s-1' + ! [ptend%q] = 'kg kg-1 s-1' + ncol = state%ncol + do k = 1, pver + ptend%q(:ncol,k,co2_fff_glo_ind) = gravit * state%rpdeldry(:ncol,k) * ac_CO2(:ncol,k) + ptend%q(:ncol,k,co2_glo_ind) = gravit * state%rpdeldry(:ncol,k) * ac_CO2(:ncol,k) + end do + end if -!=============================================================================== + end subroutine co2_cycle_set_ptend end module co2_cycle diff --git a/src/physics/cam/co2_data_flux.F90 b/src/physics/cam/co2_data_flux.F90 index 5e12525766..6ffcde55c6 100644 --- a/src/physics/cam/co2_data_flux.F90 +++ b/src/physics/cam/co2_data_flux.F90 @@ -1,164 +1,276 @@ module co2_data_flux -!------------------------------------------------------------------------------- -! utilities for reading and interpolating co2 surface fluxes -!------------------------------------------------------------------------------- - - use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl - use input_data_utils, only: time_coordinate + !------------------------------------------------------------------------------- + ! read and interpolate co2 fossil fuel surface flux + !------------------------------------------------------------------------------- + + use shr_kind_mod, only: r8=>shr_kind_r8, cl=>shr_kind_cl, cs=> shr_kind_cs + use ESMF, only: ESMF_Mesh, ESMF_Finalize, ESMF_LogFoundError + use ESMF, only: ESMF_END_ABORT, ESMF_LOGERR_PASSTHRU, ESMF_END_ABORT + use cam_logfile, only: iulog use cam_abortutils, only: endrun + use spmd_utils, only: iam, masterproc + use dshr_strdata_mod, only: shr_strdata_type implicit none private - save ! Public interfaces - public co2_data_flux_type - - public co2_data_flux_init - public co2_data_flux_advance + public :: co2_data_flux_readnl + public :: co2_data_flux_advance -!------------------------------------------------------------------------------- + ! Private interfaces + private :: co2_data_flux_init type :: co2_data_flux_type - character(len=cl) :: filename - character(len=cl) :: varname - logical :: initialized - type(time_coordinate) :: time_coord - real(r8), pointer :: co2bdy(:,:,:) ! bracketing data (pcols,begchunk:endchunk,2) - real(r8), pointer :: co2flx(:,:) ! Interpolated output (pcols,begchunk:endchunk) + character(len=cs) :: varname = "CO2_flux" + real(r8), pointer :: co2flx(:,:) ! Interpolated output (pcols,begchunk:endchunk) end type co2_data_flux_type + type(co2_data_flux_type), public :: data_flux_fuel - ! dimension names for physics grid (physgrid) - logical :: dimnames_set = .false. - character(len=8) :: dim1name, dim2name - -!=============================================================================== -contains -!=============================================================================== - -subroutine co2_data_flux_init (input_file, varname, xin) - -!------------------------------------------------------------------------------- -! Initialize co2_data_flux_type instance -! including initial read of input and interpolation to the current timestep -!------------------------------------------------------------------------------- + type(shr_strdata_type) :: sdat_co2 - use ioFileMod, only: getfil - use ppgrid, only: begchunk, endchunk, pcols - use cam_grid_support, only: cam_grid_id, cam_grid_check - use cam_grid_support, only: cam_grid_get_dim_names + character(len=cl) :: co2flux_fuel_datafile = 'unset' ! co2 flux from fossil fuel + character(len=cl) :: co2flux_fuel_meshfile = 'unset' ! ESMF mesh corresponding to co2flux_fuel_datafile + integer :: co2flux_fuel_year_first = -999 ! first year in stream to use + integer :: co2flux_fuel_year_last = -999 ! last year in stream to use + integer :: co2flux_fuel_year_align = -999 ! align stream_year_first + character(len=cs) :: co2flux_fuel_tintalgo = 'unset' ! time interpolation [linear, lower, upper] + character(len=cs) :: co2flux_fuel_taxmode = 'unset' ! time extrapolation [cycle, extend or limit] - ! Arguments - character(len=*), intent(in) :: input_file - character(len=*), intent(in) :: varname - type(co2_data_flux_type), intent(inout) :: xin + logical :: debug = .false. - ! Local variables - integer :: grid_id - real(r8) :: dtime - character(len=*), parameter :: subname = 'co2_data_flux_init' - !---------------------------------------------------------------------------- + logical :: first_advance_call = .true. - if (.not. dimnames_set) then - grid_id = cam_grid_id('physgrid') - if (.not. cam_grid_check(grid_id)) then - call endrun(subname // ': ERROR: no "physgrid" grid') - endif - call cam_grid_get_dim_names(grid_id, dim1name, dim2name) - dimnames_set = .true. - end if - - call getfil(input_file, xin%filename) - xin%varname = varname - xin%initialized = .false. - - dtime = 1.0_r8 - 200.0_r8 / 86400.0_r8 - call xin%time_coord%initialize(input_file, delta_days=dtime) - - allocate( xin%co2bdy(pcols,begchunk:endchunk,2), & - xin%co2flx(pcols,begchunk:endchunk) ) - - call co2_data_flux_advance(xin) - - xin%initialized = .true. - -end subroutine co2_data_flux_init + character(len=*),parameter :: u_FILE_u = __FILE__ +!=============================================================================== +contains !=============================================================================== -subroutine co2_data_flux_advance (xin) - -!------------------------------------------------------------------------------- -! Advance the contents of a co2_data_flux_type instance -! including reading new data, if necessary -!------------------------------------------------------------------------------- - - use cam_pio_utils, only: cam_pio_openfile - use ncdio_atm, only: infld - use pio, only: file_desc_t, pio_nowrite, pio_closefile - use ppgrid, only: begchunk, endchunk, pcols - - ! Arguments - type(co2_data_flux_type), intent(inout) :: xin - - ! Local variables - character(len=*), parameter :: subname = 'co2_data_flux_advance' - logical :: read_data - integer :: indx2_pre_adv - type(file_desc_t) :: fh_co2_data_flux - logical :: found - - !---------------------------------------------------------------------------- - - read_data = xin%time_coord%read_more() .or. .not. xin%initialized - - indx2_pre_adv = xin%time_coord%indxs(2) - - call xin%time_coord%advance() - - if ( read_data ) then - - call cam_pio_openfile(fh_co2_data_flux, trim(xin%filename), PIO_NOWRITE) - - ! read time-level 1 - ! skip the read if the needed vals are present in time-level 2 - if (xin%initialized .and. xin%time_coord%indxs(1) == indx2_pre_adv) then - xin%co2bdy(:,:,1) = xin%co2bdy(:,:,2) - else - call infld(trim(xin%varname), fh_co2_data_flux, dim1name, dim2name, & - 1, pcols, begchunk, endchunk, xin%co2bdy(:,:,1), found, & - gridname='physgrid', timelevel=xin%time_coord%indxs(1)) - if (.not. found) then - call endrun(subname // ': ERROR: ' // trim(xin%varname) // ' not found') - endif - endif - - ! read time-level 2 - call infld(trim(xin%varname), fh_co2_data_flux, dim1name, dim2name, & - 1, pcols, begchunk, endchunk, xin%co2bdy(:,:,2), found, & - gridname='physgrid', timelevel=xin%time_coord%indxs(2)) - if (.not. found) then - call endrun(subname // ': ERROR: ' // trim(xin%varname) // ' not found') - endif - - call pio_closefile(fh_co2_data_flux) - endif - - ! interpolate between time-levels - ! If time:bounds is in the dataset, and the dataset calendar is compatible with CAM's, - ! then the time_coordinate class will produce time_coord%wghts(2) == 0.0, - ! generating fluxes that are piecewise constant in time. - - if (xin%time_coord%wghts(2) == 0.0_r8) then - xin%co2flx(:,:) = xin%co2bdy(:,:,1) - else - xin%co2flx(:,:) = xin%co2bdy(:,:,1) + & - xin%time_coord%wghts(2) * (xin%co2bdy(:,:,2) - xin%co2bdy(:,:,1)) - endif + subroutine co2_data_flux_readnl(nlfile) + + !-------------------------------------------- + ! Purpose: Read co2_ffuel_nl namelist group. + !-------------------------------------------- + use namelist_utils, only: find_group_name + use spmd_utils, only: masterproc, mpicom, masterprocid + use spmd_utils, only: mpi_logical, mpi_character, mpi_integer + use cam_logfile, only: iulog + use cam_abortutils, only: endrun + use cam_pio_utils, only: cam_pio_openfile + use string_utils, only: int2str + use pio, only: PIO_BCAST_ERROR, PIO_NOERR, PIO_NOWRITE + use pio, only: file_desc_t, pio_seterrorhandling, pio_inq_varid + use pio, only: pio_closefile + + ! Arguments + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! Local variables + integer :: unitn, ierr + type(file_desc_t) :: fileid + integer :: err_handling + integer :: varid + logical :: use_time_bnds + character(len=*), parameter :: subname = 'co2_cycle_readnl' + + namelist /co2_ffuel_nl/ & + co2flux_fuel_datafile, & ! input fuel dataset + co2flux_fuel_meshfile, & ! ESMF mesh file for input dataset + co2flux_fuel_year_first, & ! first year in stream to use + co2flux_fuel_year_last, & ! last year in stream to use + co2flux_fuel_year_align, & ! align stream_year_first + co2flux_fuel_tintalgo, & ! time interpolation [linear, lower, upper] + co2flux_fuel_taxmode ! time extraploation [cycle, extend or limit] + !-------------------------------------------- + + if (masterproc) then + open( newunit=unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'co2_ffuel_nl', status=ierr) + if (ierr == 0) then + read(unitn, co2_ffuel_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(subname // ':: ERROR reading co2_ffuel_nl namelist') + end if + end if + close(unitn) + end if + + call mpi_bcast(co2flux_fuel_datafile, len(co2flux_fuel_datafile), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_datafile "//trim(co2flux_fuel_datafile)) + call mpi_bcast(co2flux_fuel_meshfile, len(co2flux_fuel_meshfile), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_meshfile "//trim(co2flux_fuel_meshfile)) + call mpi_bcast(co2flux_fuel_year_first, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_year_first "//int2str(co2flux_fuel_year_first)) + call mpi_bcast(co2flux_fuel_year_last, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_year_last "//int2str(co2flux_fuel_year_last)) + call mpi_bcast(co2flux_fuel_year_align, 1, mpi_integer, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_year_align "//int2str(co2flux_fuel_year_align)) + call mpi_bcast(co2flux_fuel_tintalgo, len(co2flux_fuel_tintalgo), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_tintalgo "//trim(co2flux_fuel_tintalgo)) + call mpi_bcast(co2flux_fuel_taxmode, len(co2flux_fuel_taxmode), mpi_character, masterprocid, mpicom, ierr) + if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: co2flux_fuel_taxmode "//trim(co2flux_fuel_taxmode)) + + if (masterproc) then + write(iulog, '(2a)') "co2flux_fuel_datafile = ", trim(co2flux_fuel_datafile) + write(iulog, '(2a)') "co2flux_fuel_meshfile = ", trim(co2flux_fuel_meshfile) + write(iulog, '(a,i0)') "co2flux_fuel_year_first = ", co2flux_fuel_year_first + write(iulog, '(a,i0)') "co2flux_fuel_year_last = ", co2flux_fuel_year_last + write(iulog, '(a,i0)') "co2flux_fuel_year_align = ", co2flux_fuel_year_align + write(iulog, '(2a)') "co2flux_fuel_tintalgo = ", trim(co2flux_fuel_tintalgo) + write(iulog, '(2a)') "co2flux_fuel_taxmode = ", trim(co2flux_fuel_taxmode) + end if + + ! Overwrite co2flux_fuel_tintalgo if it is set to 'unset' + ! Check if the data file has a time_bnds variable and if so set the time interpolation + ! type to 'nearest' otherwise set it to 'linear' + + if (trim(co2flux_fuel_tintalgo) == 'unset') then + call cam_pio_openfile( fileid, co2flux_fuel_datafile, PIO_NOWRITE ) + call pio_seterrorhandling( fileid, PIO_BCAST_ERROR, oldmethod=err_handling ) + ierr = pio_inq_varid( fileid, 'time_bnds', varid ) + call pio_seterrorhandling( fileid, err_handling) + use_time_bnds = (ierr == PIO_NOERR) + if (use_time_bnds) then + co2flux_fuel_tintalgo = 'nearest' + else + co2flux_fuel_tintalgo = 'linear' + end if + call pio_closefile( fileid ) + end if + + end subroutine co2_data_flux_readnl + + !=============================================================================== + subroutine co2_data_flux_init () + + !-------------------------------------------- + ! Initialize co2_data_flux_type instance + ! including initial read of input and interpolation to the current timestep + !-------------------------------------------- + + use ppgrid, only: begchunk, endchunk, pcols + use cam_esmf_mod, only: model_mesh, model_clock + use error_messages, only: alloc_err + use dshr_strdata_mod, only: shr_strdata_init_from_inline + + ! Local variables + integer :: istat + integer :: rc + character(len=*), parameter :: subname = 'co2_data_flux_init' + !-------------------------------------------- + + ! Allocate data_flux_fuel%co2flx + allocate( data_flux_fuel%co2flx(pcols,begchunk:endchunk), stat=istat) + call alloc_err(istat, subname, 'data_flux_fuel%co2flx', pcols*(endchunk-begchunk+1)) + + ! Initialize sdat_co2 + call shr_strdata_init_from_inline(sdat_co2, & + my_task = iam, & + logunit = iulog, & + compname = 'ATM', & + model_clock = model_clock, & + model_mesh = model_mesh, & + stream_meshfile = trim(co2flux_fuel_meshfile), & + stream_filenames = (/co2flux_fuel_datafile/), & + stream_yearFirst = co2flux_fuel_year_first, & + stream_yearLast = co2flux_fuel_year_last, & + stream_yearAlign = co2flux_fuel_year_align, & + stream_fldlistFile = (/data_flux_fuel%varname/), & + stream_fldListModel = (/data_flux_fuel%varname/), & + stream_lev_dimname = 'null', & + stream_mapalgo = 'consf', & + stream_offset = 0, & + stream_taxmode = trim(co2flux_fuel_taxmode), & + stream_dtlimit = 1.0e30_r8, & + stream_tintalgo = trim(co2flux_fuel_tintalgo), & + stream_name = 'CO2 forcing data ', & + rc = rc) + call chkrc(rc,__LINE__,u_FILE_u) + + end subroutine co2_data_flux_init + + !=============================================================================== + subroutine co2_data_flux_advance() + + !------------------------------------------------------------------------------- + ! Advance the contents of a co2_data_flux_type sdat (map and interpolate in time) + !------------------------------------------------------------------------------- + + use dshr_methods_mod , only : dshr_fldbun_getfldptr + use dshr_strdata_mod , only : shr_strdata_advance + use ppgrid , only : begchunk, endchunk + use phys_grid , only : get_ncols_p + use time_manager , only : get_curr_date + use cam_esmf_mod , only : cam_esmf_global_sum + + ! Local variables + integer :: icol,lchnk,gindx + integer :: year ! year (0, ...) for nstep+1 + integer :: mon ! month (1, ..., 12) for nstep+1 + integer :: day ! day of month (1, ..., 31) for nstep+1 + integer :: sec ! seconds into current date for nstep+1 + integer :: mcdate ! Current model date (yyyymmdd) + integer :: rc + real(r8) :: global_sum_model, global_sum_mesh + real(r8), pointer :: dataptr1d(:) + !---------------------------------------------------------------------------- + + ! Initialize stream data type for fossil fuel read + if (first_advance_call) then + call co2_data_flux_init() + first_advance_call = .false. + end if + + ! Advance sdat stream + call get_curr_date(year, mon, day, sec) + mcdate = year*10000 + mon*100 + day + call shr_strdata_advance(sdat_co2, ymd=mcdate, tod=sec, logunit=iulog, istr='co2_advance', rc=rc) + call chkrc(rc,__LINE__,u_FILE_u) + + ! Get pointer for stream data that is time and spatially interpolated to model time and grid + call dshr_fldbun_getFldPtr(sdat_co2%pstrm(1)%fldbun_model, data_flux_fuel%varname, fldptr1=dataptr1d, rc=rc) + call chkrc(rc,__LINE__,u_FILE_u) + + if (debug) then + if (masterproc) then + write(iulog,*) + write(iulog,'(a)')'Calling cam_esmf_global_sum from co2_data_flux' + end if + call cam_esmf_global_sum(trim(data_flux_fuel%varname), dataptr1d, & + global_sum_model, global_sum_mesh, rc) + call chkrc(rc,__LINE__,u_FILE_u) + write(iulog,'(a)') 'Global sum for forcing field '//trim(data_flux_fuel%varname) + write(iulog,'(a,d20.10)') ' global sum with model areas = ',global_sum_model + write(iulog,'(a,d20.10)') ' global sum with mesh areas = ',global_sum_mesh + end if + + gindx = 1 + do lchnk = begchunk,endchunk + do icol = 1,get_ncols_p(lchnk) + data_flux_fuel%co2flx(icol,lchnk) = dataptr1d(gindx) + gindx = gindx + 1 + end do + end do + + end subroutine co2_data_flux_advance + + !================================================================ + subroutine chkrc(rc, line, file) + use ESMF, only: ESMF_LOGMSG_ERROR, ESMF_SUCCESS, ESMF_LogWrite + + ! Arguments + integer , intent(in) :: rc + integer , intent(in) :: line + character(len=*) , intent(in) :: file + + if ( rc /= ESMF_SUCCESS ) then + call ESMF_LogWrite('ERROR:', ESMF_LOGMSG_ERROR, line=line, file=file) + call endrun('chkrc: See ESMF log files (filenames begin with PET)') + end if + end subroutine chkrc -end subroutine co2_data_flux_advance -!=============================================================================== end module co2_data_flux diff --git a/src/utils/ioFileMod.F90 b/src/utils/ioFileMod.F90 index c013d8aa7f..01a4204b22 100644 --- a/src/utils/ioFileMod.F90 +++ b/src/utils/ioFileMod.F90 @@ -3,13 +3,12 @@ module ioFileMod ! ! Purpose: ! -! Input/Output file manipulations. Mind file on archival system, or local -! disk etc. +! Input/Output file manipulations. ! ! Author: Mariana Vertenstein ! !--------------------------------------------------------------------- - + use shr_kind_mod, only: r8 => shr_kind_r8 use cam_abortutils, only: endrun use spmd_utils, only: masterproc @@ -33,9 +32,9 @@ module ioFileMod !======================================================================= contains !======================================================================= - + subroutine getfil(fulpath, locfn, iflag, lexist) - + ! -------------------------------------------------------------------- ! Determine whether file is on local disk. ! . first check current working directory @@ -44,7 +43,7 @@ subroutine getfil(fulpath, locfn, iflag, lexist) ! to 1 overrides this behavior, and in that case the optional lexist ! arg is used to return status of whether the file was found or not. ! -------------------------------------------------------------------- - + ! ------------------------ arguments ----------------------------------- character(len=*), intent(in) :: fulpath ! full pathname on local disk character(len=*), intent(out) :: locfn ! local file name if found in working directory, @@ -54,7 +53,7 @@ subroutine getfil(fulpath, locfn, iflag, lexist) logical, optional, intent(out) :: lexist ! When iflag=1 then getfil will return whether the ! file is found or not. This flag is set .true. ! if the file is found, otherwise .false. - + ! ------------------------ local variables --------------------------- integer :: i ! loop index integer :: klen ! length of fulpath character string @@ -63,7 +62,7 @@ subroutine getfil(fulpath, locfn, iflag, lexist) logical :: lexist_in ! true if local file exists logical :: abort_on_failure ! -------------------------------------------------------------------- - + abort_on_failure = .true. if (present(iflag)) then if (iflag==1) abort_on_failure = .false. @@ -73,7 +72,7 @@ subroutine getfil(fulpath, locfn, iflag, lexist) ! first check if file is in current working directory. ! get local file name from full name: start at end. look for first "/" - + klen = len_trim(fulpath) i = index(fulpath, '/', back=.true.) @@ -81,7 +80,9 @@ subroutine getfil(fulpath, locfn, iflag, lexist) if (abort_on_failure) then call endrun('(GETFIL): local filename variable is too short for path length') else - if (masterproc) write(iulog,*) '(GETFIL): local filename variable is too short for path length',klen-i,maxlen + if (masterproc) then + write(iulog,'(a,i0,a,i0)') '(GETFIL): local filename variable is too short for path length: ',klen-i,' > ',maxlen + end if if (present(lexist)) lexist = .false. return end if @@ -91,23 +92,25 @@ subroutine getfil(fulpath, locfn, iflag, lexist) if (len_trim(locfn) == 0) then call endrun ('(GETFIL): local filename has zero length') else if (masterproc) then - write(iulog,*)'(GETFIL): attempting to find local file ', trim(locfn) + write(iulog,'(2a)')'(GETFIL): attempting to find local file ',trim(locfn) end if - + inquire(file=locfn, exist=lexist_in) if (present(lexist)) lexist = lexist_in if (lexist_in) then - if (masterproc) write(iulog,*) '(GETFIL): using ',trim(locfn), ' in current working directory' + if (masterproc) write(iulog,'(3a)') '(GETFIL): using ',trim(locfn),' in current working directory' return end if - + ! second check for full pathname on disk - + if (klen > maxlen) then if (abort_on_failure) then call endrun('(GETFIL): local filename variable is too short for path length') else - if (masterproc) write(iulog,*) '(GETFIL): local filename variable is too short for path length',klen,maxlen + if (masterproc) then + write(iulog,'(a,i0,a,i0)') '(GETFIL): local filename variable is too short for path length: ',klen,' > ',maxlen + end if if (present(lexist)) lexist = .false. return end if @@ -117,10 +120,10 @@ subroutine getfil(fulpath, locfn, iflag, lexist) inquire(file=locfn, exist=lexist_in) if (present(lexist)) lexist = lexist_in if (lexist_in) then - if (masterproc) write(iulog,*)'(GETFIL): using ',trim(fulpath) + if (masterproc) write(iulog,'(2a)')'(GETFIL): using ',trim(fulpath) return else - if (masterproc) write(iulog,*)'(GETFIL): all tries to get file have been unsuccessful: ',trim(fulpath) + if (masterproc) write(iulog,'(2a)')'(GETFIL): all tries to get file have been unsuccessful: ',trim(fulpath) if (abort_on_failure) then call endrun ('GETFIL: FAILED to get '//trim(fulpath)) else @@ -129,29 +132,30 @@ subroutine getfil(fulpath, locfn, iflag, lexist) endif end subroutine getfil - + !======================================================================= - - + + subroutine opnfil (locfn, iun, form, status) - + use string_utils, only: int2str + !----------------------------------------------------------------------- ! open file locfn in unformatted or formatted form on unit iun !----------------------------------------------------------------------- - + ! ------------------------ input variables --------------------------- character(len=*), intent(in):: locfn !file name integer, intent(in):: iun !fortran unit number character(len=1), intent(in):: form !file format: u = unformatted. f = formatted character(len=*), optional, intent(in):: status !file status ! -------------------------------------------------------------------- - + ! ------------------------ local variables --------------------------- integer ioe !error return from fortran open character(len=11) ft !format type: formatted. unformatted character(len=11) st !file status: old or unknown ! -------------------------------------------------------------------- - + if (len_trim(locfn) == 0) then call endrun ('(OPNFIL): local filename has zero length') endif @@ -167,16 +171,20 @@ subroutine opnfil (locfn, iun, form, status) end if open (unit=iun,file=locfn,status=st, form=ft,iostat=ioe) if (ioe /= 0) then - if(masterproc) write(iulog,*)'(OPNFIL): failed to open file ',trim(locfn), ' on unit ',iun,' ierr=',ioe - call endrun ('opnfil') + if(masterproc) then + write(iulog,'(3a,i0,a,i0)')'(OPNFIL): failed to open file ', trim(locfn), ' on unit ',iun,', ierr=',ioe + end if + call endrun('(OPNFIL): failed to open file '//trim(locfn)//' on unit '//int2str(iun)//', ierr='//int2str(ioe)) else - if(masterproc) write(iulog,*)'(OPNFIL): Successfully opened file ',trim(locfn), ' on unit= ',iun + if(masterproc) then + write(iulog,'(3a,i0)')'(OPNFIL): Successfully opened file ', trim(locfn), ' on unit = ', iun + end if end if - + return end subroutine opnfil - + !======================================================================= - - + + end module ioFileMod