diff --git a/.release/changeVersionNumbers.sh b/.release/changeVersionNumbers.sh new file mode 100755 index 00000000..82f5410a --- /dev/null +++ b/.release/changeVersionNumbers.sh @@ -0,0 +1,108 @@ +#!/bin/bash + +#EOC +#------------------------------------------------------------------------------ +# Harmonized Emissions Component (HEMCO) ! +#------------------------------------------------------------------------------ +#BOP +# +# !MODULE: changeVersionNumbers.sh +# +# !DESCRIPTION: Bash script to change the version numbers in the appropriate +# files in the HEMCO directory structure. Run this before releasing +# a new HEMCO version. +#\\ +#\\ +# !CALLING SEQUENCE: +# $ ./changeVersionNumbers.sh X.Y.Z # X.Y.Z = HEMCO version number +#EOP +#------------------------------------------------------------------------------ +#BOC + +function replace() { + + #======================================================================== + # Function to replace text in a file via sed. + # + # 1st argument: Search pattern + # 2nd argument: Replacement text + # 3rd argument: File in which to search and replace + #======================================================================== + + sed -i -e "s/${1}/${2}/" "${3}" +} + + +function exitWithError() { + + #======================================================================== + # Display and error message and exit + #======================================================================== + + echo "Could not update version numbers in ${1}... Exiting!" + exit 1 +} + + +function main() { + + #======================================================================== + # Replaces the version number in the files listed. + # + # 1st argument: New version number to use + #======================================================================== + + # New version number + version="${1}" + + # Save this directory path and change to root directory + thisDir=$(pwd -P) + cd .. + + #======================================================================== + # Update version numbers in various files + #======================================================================== + + # Pattern to match: X.Y.Z + pattern='[0-9][0-9]*\.[0-9][0-9]*\.[0-9][0-9]*' + + # List of files to replace + files=( \ + "CMakeLists.txt" \ + "docs/source/conf.py" \ + "src/Core/hco_error_mod.F90" + ) + + # Replace version numbers in files + for file in ${files[@]}; do + replace "${pattern}" "${version}" "${file}" + [[ $? -ne 0 ]] && exitWithError "${file}" + echo "HEMCO version updated to ${version} in ${file}" + done + + #======================================================================== + # Update version number and date in CHANGELOG.md + #======================================================================== + + # Pattern to match: "[Unreleased] - TBD" + pattern='\[.*Unreleased.*\].*' + date=$(date -Idate) + replace "${pattern}" "\[${version}\] - ${date}" "CHANGELOG.md" + + # Return to the starting directory + cd "${thisDir}" +} + +# --------------------------------------------------------------------------- + +# Expect 1 argument, or exit with error +if [[ $# -ne 1 ]]; then + echo "Usage: ./changeVersionNumbers.sh VERSION" + exit 1 +fi + +# Replace version numbers +main "${1}" + +# Return status +exit $? diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f96a6b9..46930265 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,25 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] - TBD +### Changed +- Updated TOMAS_Jeagle sea salt extension + +### Fixed +- Updated IsModelLevel check for CESM and WRF-GC +- Interpolation error for 8-day MODIS LAI files (removed month loop in `GetIndex2Interp`) + +## [3.7.2] - 2023-12-01 +### Added +- Script `.release/changeVersionNumbers.sh` to change version numbers before a new HEMCO release + +### Changed +- Increased netCDF variable string length from 50 to 100 + +### Fixed +- Rename `HEMCO_Config.rc.sample` to `HEMCO_Config.rc` in `createRunDir.sh` if sample is used. +- Added fix to turn off emissions extensions when `EMISSIONS` logical is false + ## [3.7.1] - 2023-10-10 ### Changed - Updated version numbers to 3.7.1 @@ -30,7 +49,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Add GEOSIT as an allowable meteorology directory name in HEMCO_Config.rc - Added `.readthedocs.yaml` file to configure ReadTheDocs builds -# Changed +### Changed - `Verbose` is now a `true/false` variable in `run/HEMCO_sa_Config.rc` and `run/HEMCO_Config.rc.sample` - HEMCO warnings are now only generated when `Verbose: true` is found in the HEMCO configuration file (no more numerical levels) - Updated GFED4 emission factors for VOCs to Andreae et al. (2019) @@ -100,14 +119,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [3.5.2] - 2022-11-29 ### Added -- Added sanitizer option for detecting memory leaks in HEMCO - standalone during build +- Added sanitizer option for detecting memory leaks in HEMCO standalone during build ### Changed - Remove unused, commented-out code in `src/Extensions/hcox_dustdead_mod.F` -- Replaced placeholder error messages in - `src/Core/hco_config_mod.F90` with more informational messages - (often including the line of the HEMCO_Config.rc in the printout) +- Replaced placeholder error messages in `src/Core/hco_config_mod.F90` with more informational messages (often including the line of the HEMCO_Config.rc in the printout) - Added improved documentation for time cycle flag `EFYO` in ReadTheDocs ### Fixed @@ -121,8 +137,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Support for MAPL 2.16 (needed by GCHP and GEOS) - Bug fix for HEMCO standalone run directory creation -- Bug fix: If HEMCO masks are specified as `lon1/lat1/lon2/lat2`, - then don't try to read from disk +- Bug fix: If HEMCO masks are specified as `lon1/lat1/lon2/lat2`, then don't try to read from disk - Documentation from the GEOS-Chem wiki (now on ReadTheDocs) - Badges for the ReadTheDocs front page - Bug fix for masking issues in MPI environment (for WRF, CESM) @@ -265,6 +280,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Bug fix for distributing emissions in the vertical dimension - New error checks in the HEMCO standalone module - Bug fix for `ifort` compiler in soil NOx extension + ### Removed - Null string character from netCDF unit string @@ -280,8 +296,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [2.1.004] - 2017-12-30 ### Added -- Updates to remove possible issues and excessive print statements when - operating in GEOS environment +- Updates to remove possible issues and excessive print statements when operating in GEOS environment - Fixed possible tracer ID mismatch in sea salt extension - New option to normalize MEGAN LAI, HEMCO diagnostics - Now write multiple time slices into one file @@ -364,8 +379,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [1.1.009] - 2015-09-10 ### Added -- Bug fixes to allow specifying flexible diagnostics output - frequencies. +- Bug fixes to allow specifying flexible diagnostics output frequencies. ## [1.1.008] - 2015-07-06 ### Added diff --git a/CMakeLists.txt b/CMakeLists.txt index ea235e01..39638d1e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ # HEMCO/CMakeLists.txt cmake_minimum_required(VERSION 3.5) -project(HEMCO VERSION 3.7.1 LANGUAGES Fortran) +project(HEMCO VERSION 3.7.2 LANGUAGES Fortran) # Reminder: Make sure to also update version in src/Core/hco_error_mod.F90 #----------------------------------------------------------------------------- diff --git a/SUPPORT.md b/SUPPORT.md index 369a2dd0..fa4f5d40 100644 --- a/SUPPORT.md +++ b/SUPPORT.md @@ -13,7 +13,7 @@ We use GitHub issues to support user questions. To ask a question, **[open a new ## What type of support can I expect? -We will be happy to assist you in resolving bugs and technical issues that arise when compiling or running HEMCO. User support and outreach is an important part of our mission to support the [International GEOS-Chem User Community](https://geoschem.github.io/geos-chem-people-projects-map/). +We will be happy to assist you in resolving bugs and technical issues that arise when compiling or running HEMCO. User support and outreach is an important part of our mission to support the [International GEOS-Chem User Community](https://geoschem.github.io/people.html). Even though we can assist in several ways, we cannot possibly do everything. We rely on HEMCO users being resourceful and willing to try to resolve problems on their own to the greatest extent possible. diff --git a/docs/source/conf.py b/docs/source/conf.py index 325803c6..c7ccc0b0 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -23,7 +23,7 @@ author = 'GEOS-Chem Support Team' # The full version, including alpha/beta/rc tags -release = '3.7.1' +release = '3.7.2' # -- General configuration --------------------------------------------------- diff --git a/docs/source/hco-ref-guide/hemco-config.rst b/docs/source/hco-ref-guide/hemco-config.rst index 6bcc4186..876a9b02 100644 --- a/docs/source/hco-ref-guide/hemco-config.rst +++ b/docs/source/hco-ref-guide/hemco-config.rst @@ -231,7 +231,7 @@ year, month, day, or hour, and to scale emissions to a given value: If omitted, the emisison month will be set to the model simulation hour. -.. option:: EmissScale_ +.. option:: EmisScale_ Optional argument to define a uniform scale factor that will be applied across all inventories, categories, hierarchies, and diff --git a/run/createRunDir.sh b/run/createRunDir.sh index 8035a1f9..84af3edc 100755 --- a/run/createRunDir.sh +++ b/run/createRunDir.sh @@ -202,6 +202,11 @@ while [ "$valid_path" -eq 0 ]; do if [[ "$hco_config_file" == *".rc"* ]]; then hco_config_dir=$(dirname $hco_config_file) fi + + if [[ "$hco_config_file" == "./HEMCO_Config.rc.sample" ]]; then + mv "./HEMCO_Config.rc.sample" "./HEMCO_Config.rc" + hco_config_file="./HEMCO_Config.rc" + fi done diff --git a/src/Core/hco_config_mod.F90 b/src/Core/hco_config_mod.F90 index 91d4d503..1dfca34a 100644 --- a/src/Core/hco_config_mod.F90 +++ b/src/Core/hco_config_mod.F90 @@ -635,7 +635,7 @@ SUBROUTINE Config_ReadCont( HcoConfig, IU_HCO, CFDIR, & LOGICAL :: Found CHARACTER(LEN= 63) :: cName CHARACTER(LEN=255) :: srcFile - CHARACTER(LEN= 50) :: srcVar + CHARACTER(LEN=100) :: srcVar CHARACTER(LEN= 31) :: srcTime CHARACTER(LEN= 31) :: TmCycle CHARACTER(LEN= 1) :: WildCard @@ -1816,9 +1816,9 @@ SUBROUTINE ExtSwitch2Buffer( HcoConfig, IU_HCO, EOF, RC ) ! ! !USES: ! - USE HCO_CHARPAK_Mod, ONLY : STRREPL, STRSPLIT, TRANLC - USE HCO_EXTLIST_MOD, ONLY : AddExt, AddExtOpt, HCO_GetOpt - USE HCO_EXTLIST_MOD, ONLY : GetExtNr + USE HCO_CHARPAK_Mod, ONLY : STRREPL, STRSPLIT, TRANLC + USE HCO_EXTLIST_MOD, ONLY : AddExt, AddExtOpt, HCO_GetOpt + USE HCO_EXTLIST_MOD, ONLY : GetExtNr, GetExtOpt ! ! !INPUT PARAMETERS: ! @@ -1841,6 +1841,7 @@ SUBROUTINE ExtSwitch2Buffer( HcoConfig, IU_HCO, EOF, RC ) ! INTEGER :: I, N, Idx, ExtNr LOGICAL :: Enabled, NewExt + LOGICAL :: DoEmis, Found, LTMP CHARACTER(LEN=255) :: loc CHARACTER(LEN=512) :: msg CHARACTER(LEN=1023) :: OPTS @@ -1857,6 +1858,10 @@ SUBROUTINE ExtSwitch2Buffer( HcoConfig, IU_HCO, EOF, RC ) loc = 'ExtSwitch2Buffer (hco_config_mod.F90)' ExtNr = -1 + ! Initialize + DoEmis= .TRUE. + Found = .FALSE. + ! Do until exit DO @@ -1890,6 +1895,18 @@ SUBROUTINE ExtSwitch2Buffer( HcoConfig, IU_HCO, EOF, RC ) RETURN ENDIF ENDIF + + ! Check if EMISSIONS setting is found. If so, overwrite DoEmis. + IF ( .not. Found ) THEN + CALL GetExtOpt( HcoConfig, -999, 'EMISSIONS', & + OptValBool=LTMP, FOUND=Found, RC=RC ) + IF ( RC /= HCO_SUCCESS ) THEN + msg = 'Error encountered in "GetExtOpt( EMISSIONS )"!' + CALL HCO_Error( msg, RC, ThisLoc=loc ) + RETURN + ENDIF + IF ( Found ) DoEmis = LTMP + ENDIF CYCLE ENDIF @@ -1942,6 +1959,11 @@ SUBROUTINE ExtSwitch2Buffer( HcoConfig, IU_HCO, EOF, RC ) Enabled = .FALSE. ENDIF + ! Disable extension if EMISSIONS logical is false + IF ( .not. DoEmis ) THEN + Enabled = .FALSE. + ENDIF + ! Register extension name, number and species ! idx is the position of the species names idx = idx+1 diff --git a/src/Core/hco_error_mod.F90 b/src/Core/hco_error_mod.F90 index b242dc1c..7bf556e9 100644 --- a/src/Core/hco_error_mod.F90 +++ b/src/Core/hco_error_mod.F90 @@ -105,7 +105,7 @@ MODULE HCO_Error_Mod #endif ! HEMCO version number. - CHARACTER(LEN=12), PARAMETER, PUBLIC :: HCO_VERSION = '3.7.1' + CHARACTER(LEN=12), PARAMETER, PUBLIC :: HCO_VERSION = '3.7.2' INTERFACE HCO_Error MODULE PROCEDURE HCO_ErrorNoErr diff --git a/src/Core/hco_extlist_mod.F90 b/src/Core/hco_extlist_mod.F90 index 6c382e54..88f5cb29 100644 --- a/src/Core/hco_extlist_mod.F90 +++ b/src/Core/hco_extlist_mod.F90 @@ -1567,7 +1567,7 @@ SUBROUTINE HCO_SetDefaultToken( CF, RC ) ELSE IF ( TRIM(CF%MetField) == 'GEOSIT' ) THEN DEF_MET_UC = 'GEOSIT' DEF_MET_LC = 'geosit' - DEF_CN_YR = '2018' ! Constant met fld year + DEF_CN_YR = '1998' ! Constant met fld year DEF_NC_VER = 'nc' ! NetCDF extension ENDIF diff --git a/src/Core/hco_types_mod.F90 b/src/Core/hco_types_mod.F90 index 34af66cf..ab990b6d 100644 --- a/src/Core/hco_types_mod.F90 +++ b/src/Core/hco_types_mod.F90 @@ -338,7 +338,7 @@ MODULE HCO_TYPES_MOD !------------------------------------------------------------------------- TYPE :: FileData CHARACTER(LEN=255) :: ncFile ! file path+name - CHARACTER(LEN= 50) :: ncPara ! file parameter + CHARACTER(LEN=100) :: ncPara ! file parameter INTEGER :: ncYrs(2) ! year range INTEGER :: ncMts(2) ! month range INTEGER :: ncDys(2) ! day range diff --git a/src/Core/hcoio_read_std_mod.F90 b/src/Core/hcoio_read_std_mod.F90 index 02eb44b9..76f35994 100644 --- a/src/Core/hcoio_read_std_mod.F90 +++ b/src/Core/hcoio_read_std_mod.F90 @@ -695,12 +695,39 @@ SUBROUTINE HCOIO_Read( HcoState, Lct, RC ) ! going to 72 levels. Otherwise, use MESSy (nbalasus, 8/24/2023). IF ( Lct%Dct%Dta%Levels == 0 ) THEN +#if defined( MODEL_CESM ) || defined( MODEL_WRF ) + + ! In WRF/CESM, IsModelLevel has a different meaning of "GEOS-Chem levels" + ! because the models in WRF and CESM are user-defined and thus fixed input + ! files would never be on the model level. In this case, a check is added + ! in order to match the file with known GEOS-Chem levels, and if so, the + ! data will be handled later in this file accordingly to be vertically + ! regridded to the runtime model levels using MESSy. + ! This fixes a regression from the vertical regridding fixes in 3.7.1. + ! + ! The meaning of "is model levels" in WRF and CESM are different. + ! Model levels can be changed and thus data is never on the model level. + ! In this case, IsModelLevel means that the data is on standard + ! GEOS-Chem levels, and if so, the data will be handled accordingly + ! using a hard-coded set of GEOS-Chem levels to be interpolated using MESSy. + ! (hplin, 10/15/23) + IF ( TRIM(LevUnit) == "level" .or. TRIM(LevUnit) == "GEOS-Chem level" ) THEN + ! the below check will be obsolete and is unmaintainable, but would be consistent with ModelLev_Check. + ! it is more robust to check for the explicit intention of LevUnit + ! nlev == 47 .or. nlev == 48 .or. nlev == 36 .or. nlev == 72 .or. nlev == 73 ) THEN + IsModelLevel = .true. + ENDIF + +#else + CALL ModelLev_Check( HcoState, nlev, IsModelLevel, RC ) IF ( RC /= HCO_SUCCESS ) THEN CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC ) RETURN ENDIF +#endif + ! Set level indeces to be read lev1 = 1 lev2 = nlev diff --git a/src/Core/hcoio_util_mod.F90 b/src/Core/hcoio_util_mod.F90 index d14aa7f4..d84997dd 100644 --- a/src/Core/hcoio_util_mod.F90 +++ b/src/Core/hcoio_util_mod.F90 @@ -1178,27 +1178,6 @@ SUBROUTINE GetIndex2Interp ( HcoState, Lct, & IF ( tidx2 > 0 ) EXIT ENDDO - ! Repeat above but now only modify month. - IF ( tidx2 < 0 ) THEN - tmpYMDhm = availYMDhm(tidx1) - DO - ! Increase by one month - tmpYMDhm = tmpYMDhm + 1.0e6_dp - - ! Exit if we are beyond available dates - IF ( tmpYMDhm > availYMDhm(nTime) ) EXIT - - ! Check if there is a time slice with that date - DO I = tidx1,nTime - IF ( ABS( tmpYMDhm - availYMDhm(I) ) < EPSILON ) THEN - tidx2 = I - EXIT - ENDIF - ENDDO - IF ( tidx2 > 0 ) EXIT - ENDDO - ENDIF - ! Repeat above but now only modify day IF ( tidx2 < 0 ) THEN tmpYMDhm = availYMDhm(tidx1) diff --git a/src/Extensions/hcox_tomas_dustdead_mod.F b/src/Extensions/hcox_tomas_dustdead_mod.F index 1ad5c977..3c944201 100644 --- a/src/Extensions/hcox_tomas_dustdead_mod.F +++ b/src/Extensions/hcox_tomas_dustdead_mod.F @@ -517,15 +517,15 @@ SUBROUTINE HCOX_TOMAS_DustDead_Run( ExtState, HcoState, RC ) ! Get the proper species name !%%% This is a kludge, do better later %%% - IF ( N==1 ) SpcName = 'NK1' - IF ( N==2 ) SpcName = 'NK2' - IF ( N==3 ) SpcName = 'NK3' - IF ( N==4 ) SpcName = 'NK4' - IF ( N==5 ) SpcName = 'NK5' - IF ( N==6 ) SpcName = 'NK6' - IF ( N==7 ) SpcName = 'NK7' - IF ( N==8 ) SpcName = 'NK8' - IF ( N==9 ) SpcName = 'NK9' + IF ( N==1 ) SpcName = 'NK01' + IF ( N==2 ) SpcName = 'NK02' + IF ( N==3 ) SpcName = 'NK03' + IF ( N==4 ) SpcName = 'NK04' + IF ( N==5 ) SpcName = 'NK05' + IF ( N==6 ) SpcName = 'NK06' + IF ( N==7 ) SpcName = 'NK07' + IF ( N==8 ) SpcName = 'NK08' + IF ( N==9 ) SpcName = 'NK09' IF ( N==10 ) SpcName = 'NK10' IF ( N==11 ) SpcName = 'NK11' IF ( N==12 ) SpcName = 'NK12' diff --git a/src/Extensions/hcox_tomas_jeagle_mod.F90 b/src/Extensions/hcox_tomas_jeagle_mod.F90 index 8e64b5e8..a53a58ad 100644 --- a/src/Extensions/hcox_tomas_jeagle_mod.F90 +++ b/src/Extensions/hcox_tomas_jeagle_mod.F90 @@ -59,12 +59,40 @@ MODULE HCOX_TOMAS_Jeagle_Mod INTEGER :: ExtNr ! HEMCO extension # REAL(dp) :: TOMAS_COEF ! Seasalt emiss coeff. + !Tracer IDs + LOGICAL :: EmitSnowSS ! Calculate sea salt emission blowing snow + + ! Arrays INTEGER, ALLOCATABLE :: HcoIDs (: ) ! HEMCO species ID's REAL(dp), POINTER :: TOMAS_DBIN(: ) ! TOMAS bin width REAL(dp), POINTER :: DRFAC (: ) ! TOMAS area? - REAL(dp), POINTER :: TC1 (:,:,:,:) ! Aerosol mass - REAL(dp), POINTER :: TC2 (:,:,:,:) ! Aerosol number + REAL(dp), POINTER :: TC1 (:,:,:,:) ! Aerosol number + REAL(dp), POINTER :: TC2 (:,:,:,:) ! Aerosol mass + LOGICAL :: ColdSST ! Flag to correct SSA emissions over cold waters + + ! Scale factors + REAL*8 :: NSLNT_FYI ! North Hemisphere snow salinity on first year ice (FYI) (psu) + REAL*8 :: NSLNT_MYI ! North Hemisphere snow salinity on multiyear ice (MYI) (psu) + REAL*8 :: SSLNT_FYI ! South Hemisphere snow salinity on FYI (psu) + REAL*8 :: SSLNT_MYI ! South Hemisphere snow salinity on MYI (psu) + REAL*8 :: NAGE ! North Hemisphere snow age (days) + REAL*8 :: SAGE ! South Hemisphere snow age (days) + REAL*8 :: NP ! number of particle per snowflake + + !Module variables + REAL*8, POINTER :: SS_DEN(:) ! densities + REAL*8, POINTER :: F_DI_N_FYI(:,:) ! add for blowing snow for NH + REAL*8, POINTER :: F_DI_N_MYI(:,:) ! add for blowing snow for NH + REAL*8, POINTER :: F_DI_S_FYI(:,:) ! add for blowing snow for SH + REAL*8, POINTER :: F_DI_S_MYI(:,:) ! add for blowing snow for SH + REAL*8, POINTER :: F_DN_N_FYI(:,:) ! add for blowing snow for NH + REAL*8, POINTER :: F_DN_N_MYI(:,:) ! add for blowing snow for NH + REAL*8, POINTER :: F_DN_S_FYI(:,:) ! add for blowing snow for SH + REAL*8, POINTER :: F_DN_S_MYI(:,:) ! add for blowing snow for SH + + !Number densities + REAL(sp), POINTER :: MULTIICE(:,:) => NULL() ! add for blowing snow TYPE(MyInst), POINTER :: NextInst => NULL() END TYPE MyInst @@ -72,6 +100,9 @@ MODULE HCOX_TOMAS_Jeagle_Mod ! Pointer to instances TYPE(MyInst), POINTER :: AllInst => NULL() + !INTEGER, PARAMETER :: NR_MAX = 200 ! max. # of bins per mode + INTEGER, PARAMETER :: NR_MAX = 4000 ! max. # of bins per mode + CONTAINS !EOC !------------------------------------------------------------------------------ @@ -95,6 +126,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) USE HCO_GeoTools_Mod, ONLY : HCO_LandType USE HCO_FluxArr_mod, ONLY : HCO_EmisAdd USE HCO_State_Mod, ONLY : HCO_GetHcoID + USE HCO_Calc_Mod, ONLY : HCO_EvalFld ! ! !INPUT PARAMETERS: ! @@ -106,7 +138,32 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) INTEGER, INTENT(INOUT) :: RC ! Success or failure? ! ! !REMARKS: -! +! References +! ============================================================================ +! (1 ) Chin, M., P. Ginoux, S. Kinne, B. Holben, B. Duncan, R. Martin, +! J. Logan, A. Higurashi, and T. Nakajima, "Tropospheric aerosol +! optical thickness from the GOCART model and comparisons with +! satellite and sunphotometers measurements", J. Atmos Sci., 2001. +! (2 ) Gong, S., L. Barrie, and J.-P. Blanchet, "Modeling sea-salt +! aerosols in the atmosphere. 1. Model development", J. Geophys. Res., +! v. 102, 3805-3818, 1997. +! (3 ) Gong, S. L., "A parameterization of sea-salt aerosol source function +! for sub- and super-micron particles", Global Biogeochem. Cy., 17(4), +! 1097, doi:10.1029/2003GB002079, 2003. +! (4 ) Jaegle, L., P.K. Quinn, T.S. Bates, B. Alexander, J.-T. Lin, "Global +! distribution of sea salt aerosols: New constraints from in situ and +! remote sensing observations", Atmos. Chem. Phys., 11, 3137-3157, +! doi:10.5194/acp-11-3137-2011. +! (5 ) Huang, J., Jaeglé, L., "Wintertime enhancements of sea salt aerosol in +! polar regions consistent with a sea ice source from blowing snow." +! Atmos. Chem. Phys. 17, 3699–3712. https://doi.org/10.5194/acp-17-3699-2017, 2017. +! (6 ) Huang, J., Jaeglé, L., Chen, Q., Alexander, B., Sherwen, T., +! Evans, M. J., Theys, N., and Choi, S. "Evaluating the impact of +! blowing snow sea salt aerosol on springtime BrO and O3 in the Arctic, +! Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2019-1094, 2020. +! (7 ) Tschudi, M., W. N. Meier, J. S. Stewart, C. Fowler, and J. Maslanik. +! "EASE-Grid Sea Ice Age, Version 4." NASA National Snow and Ice Data Center +! Distributed Active Archive Center. doi: https://doi.org/10.5067/UTAV7490FEPB., 2019. ! ! !REVISION HISTORY: ! 01 Oct 2014 - R. Yantosca - Initial version, based on TOMAS SRCSALT30 code @@ -121,9 +178,50 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) INTEGER :: I, J, L, K, HcoID REAL(sp) :: FOCEAN, W10M, DTEMIS REAL(dp) :: F100, W, A_M2, FEMIS, NUMBER, MASS, NUMBER_TOT + REAL(dp) :: NUMBER2 REAL(dp) :: rwet, dfo, B, A, SST, SCALE CHARACTER(LEN=255):: MSG, LOC + ! New variables for blowing snow (huang, 04/09/20) + REAL*8 :: SNOWSALT + REAL*8 :: FROPEN, FRFIRST + REAL*8 :: FRICTVEL, WVMR, TEMP + REAL*8 :: PRESS, P_ICE, RH_ICE + REAL*8 :: D, FK, FD + REAL*8 :: PSI, QSPRIME, UT, APRIM + REAL*8 :: QS, QSNOWICE_FYI, QSNOWICE_MYI,QBSALT, QB0 + REAL*8 :: SLNT, SLNT_FYI, SLNT_MYI + REAL*8 :: AGE, ISFROST + + ! New parameters for blowing snow (huang, 04/09/20) + REAL*8, PARAMETER :: LS = 2839d3 ! Latent heat of sublimation @ T=-30C (J/kg). + ! Varies very little with Temperature + REAL*8, PARAMETER :: RV = 461.5d0 !J kg-1 K-1 + REAL*8, PARAMETER :: RHONACL = 2160.0d0 !kg/m3 + REAL*8, PARAMETER :: RHOICE = 900.0d0 !kg/m3 + REAL*8, PARAMETER :: K0 = 2.16d-2 !J m-1 s-1 K-1 + REAL*8, PARAMETER :: A0 = 3.78407d-1 + REAL*8, PARAMETER :: A1 = -8.64089d-2 + REAL*8, PARAMETER :: A2 = -1.60570d-2 + REAL*8, PARAMETER :: A3 = 7.25516d-4 + REAL*8, PARAMETER :: A4 = -1.25650d-1 + REAL*8, PARAMETER :: A5 = 2.48430d-2 + REAL*8, PARAMETER :: A6 = -9.56871d-4 + REAL*8, PARAMETER :: A7 = 1.24600d-2 + REAL*8, PARAMETER :: A8 = 1.56862d-3 + REAL*8, PARAMETER :: A9 = -2.93002d-4 + REAL*8, PARAMETER :: A_SALT = 2.0d0 !from Mann et al. 2000 + REAL*8, PARAMETER :: B_SALT = 37.5d0 !in um + !REAL*8, PARAMETER :: DDSNOW = 2.0d0 !in um for snow particle interval + REAL*8, PARAMETER :: DDSNOW = 0.1d0 !in um for snow particle interval + LOGICAL, SAVE :: FIRST = .TRUE. + LOGICAL, SAVE :: FIRSTSAL = .TRUE. + CHARACTER(LEN=31) :: FLDNME + INTEGER :: NDAYS!, cYYYY, cMM, cDD, K + REAL(hp), TARGET :: MULTI(HcoState%NX,HcoState%NY) + REAL(hp), TARGET :: SNOWSALA (HcoState%NX,HcoState%NY) + REAL(hp), TARGET :: SNOWSALC (HcoState%NX,HcoState%NY) + REAL*8, PARAMETER :: BETHA=2.22d0 !wet diameter (80% Rel Hum) to dry diam ! Strings @@ -134,10 +232,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) REAL(dp), POINTER :: ptr3D(:,:,:) ! For debugging - !INTEGER :: ii=50, jj=10 - - ! Error handling - LOGICAL :: ERR + INTEGER :: ii=50, jj=10 !================================================================= ! SRCSALT30 begins here! @@ -159,7 +254,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) CALL InstGet ( ExtState%TOMAS_Jeagle, Inst, RC ) IF ( RC /= HCO_SUCCESS ) THEN WRITE(MSG,*) 'Cannot find TOMAS_Jeagle instance Nr. ', ExtState%TOMAS_Jeagle - CALL HCO_ERROR(MSG,RC) + CALL HCO_ERROR(MSG, RC ) RETURN ENDIF @@ -167,11 +262,26 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) Inst%TC1 = 0.0_hp Inst%TC2 = 0.0_hp + IF ( Inst%EmitSnowSS ) THEN + ! Read in distribution of multi-year sea ice from + ! remotely sensed observations of sea ice motion and sea + ! ice extent for the Arctic (Tschudi et al., 2019). For the + ! Antarctic, the multi year sea ice extent is based on the minimum + ! MERRA-2 sea ice extent of the previous summer. + + CALL HCO_EvalFld ( HcoState, 'MULTISEAICE', MULTI, RC ) + IF ( RC /= HCO_SUCCESS ) THEN + WRITE(MSG,*) 'Cannot find MULTISEAICE data for blowing snow' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + ENDIF + ! Depending on the grid resolution. 4x5 (default) doesn't need ! adjusting coeff !### Debug - !print*, 'JACK IN HCOX TOMAS Jeagle' + print*, 'IN HCOX_TOMAS_Jeagle_Mod.F90' ! Init ptr3D => NULL() @@ -183,37 +293,131 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) DO J = 1, HcoState%NY DO I = 1, HcoState%NX + Inst%TC1(I,J,:,:) = 0d0 + Inst%TC2(I,J,:,:) = 0d0 + ! Grid box surface area [m2] A_M2 = HcoState%Grid%AREA_M2%Val(I,J) - ! Get the fraction of the box over ocean - IF ( HCO_LandType( ExtState%FRLAND%Arr%Val(I,J), & - ExtState%FRLANDIC%Arr%Val(I,J), & - ExtState%FROCEAN%Arr%Val(I,J), & - ExtState%FRSEAICE%Arr%Val(I,J), & - ExtState%FRLAKE%Arr%Val(I,J) ) == 0 ) THEN - FOCEAN = 1d0 - ExtState%FRCLND%Arr%Val(I,J) - ELSE - FOCEAN = 0.d0 + ! Advance to next grid box if it's not over water or sea ice + IF ( ExtState%FROCEAN%Arr%Val(I,J)<=0d0 .and. & + ExtState%FRSEAICE%Arr%Val(I,J)<=0d0 ) CYCLE + + ! Wind speed at 10 m altitude [m/s] + W10M = SQRT( ExtState%U10M%Arr%Val(I,J)**2 & + + ExtState%V10M%Arr%Val(I,J)**2 ) + + ! Sea surface temperature in Celcius + SST = ExtState%TSKIN%Arr%Val(I,J) - 273.15d0 + + ! Limit SST to 0-30C range + SST = MAX( SST , 0d0 ) ! limit to 0C + SST = MIN( SST , 30d0 ) ! limit to 30C + + ! Empirical SST scaling factor (jaegle 5/11/11) + SCALE = 0.329d0 + 0.0904d0*SST - & + 0.00717d0*SST**2d0 + 0.000207d0*SST**3d0 + + ! Limit the SST scaling factor to 0.25 over cold SST (below 5C) + IF ( Inst%ColdSST .and. SST<= 5.0d0 ) SCALE = 0.25d0 + + ! Apply to only the open ocean fraction of the gridbox (Huang 06/12/20) + FROPEN = ExtState%FROCEAN%Arr%Val(I,J)-ExtState%FRSEAICE%Arr%Val(I,J) + IF ( FROPEN < 0d0 ) FROPEN = 0d0 + + ! Eventually apply wind scaling factor. + !SCALE = SCALE * Inst%WindScale * FROPEN + SCALE = SCALE * FROPEN + + !---------------------------------------------------------------- + ! huang, 04/09/20: Add blowing snow emissions over sea ice + !---------------------------------------------------------------- + + IF ( Inst%EmitSnowSS ) THEN + IF ( ExtState%FRSEAICE%Arr%Val(I,J) > 0d0 )THEN + ! Friction velocity [m/s] + FRICTVEL = ExtState%USTAR%Arr%Val(I,J) + ! Convert specific humidity [g H2O/kg air] to water vapor mixing ratio [v/v] + ! QV2m is in kg H2O/kg air + WVMR = ExtState%QV2M%Arr%Val(I,J) * 28.973d0 / 18.0d0 + ! Temperature at 2M in grid box (I,J) [K] + TEMP = ExtState%T2M%Arr%Val(I,J) + ! Surface pressure at grid box (I,J). Convert from [Pa] to [hPa] + PRESS = HcoState%Grid%PSFC%Val( I, J ) /100d0 + ! Calculate saturation vapor pressure over ice [in Pa] at temperature + ! TS [K] + P_ICE = 10d0**(-2663.5d0/TEMP+12.537d0) + ! Calculate relative humidity with respect to ice [%] + RH_ICE = PRESS * WVMR / (P_ICE*0.01d0) *100.0d0 + ! Limit RH to 100% + IF (RH_ICE > 100d0) RH_ICE =100.0d0 + ! Coefficient of Diffusion of water vapor in air [m2/s] + ! Parameterization of Massman, W.J. "A review of teh molecular diffusivities of + ! H2O, CO2, CH4... in air, O2 and N2 near STP" Atmos. Env., 32, 6, 1111-1127, 1998. + D = 2.178d-5*(1000d0/PRESS)*(TEMP/273.15d0)**1.81 + ! Heat conductivity and vapor diffusion terms [m s/kg] + ! Rogers and Yau "A short course in cloud physics", 1989, Eqn 9.4, with + ! RV = 461.5 [J/kg/K] Individual gas constant for water vapor + ! LS = 2839.0*1d3 [J/kg ] Latent heat of sublimation @ T=-30C + ! K = 2.16d-2 [J/(m s K)] Coeff of thermal conductivity of Air [Table 7.1 Rogers and Yau] + FK = ( LS / (RV * TEMP ) -1d0 ) * LS / (K0 * TEMP) + FD = ( RV * TEMP ) / (D * P_ICE) + ! Variable PSI [m2/s] Equation 11 from Dery and Yau (2001) + ! RHOICE = 900 kg/m3 Density of ice + PSI = (RH_ICE/100.d0 - 1d0)/(2d0 * RHOICE * (FK + FD)) + ! Convert PSI from m2/s to units of -1x10d-12 m2/s + PSI = PSI * (-1.0d12) + ! Qs prime [mm/day snow water equivalent] Equation 11 Dery and Yau (2001) + QSPRIME = A0 + A1*PSI + A2*PSI**2d0 + A3*PSI**3d0 & + + A4* W10M + A5*PSI*W10M & + + A6*W10M*PSI**2d0 + A7*W10M**2d0 & + + A8*PSI*W10M**2d0 + A9*W10M**3d0 + IF ( QSPRIME < 0.0d0 ) QSPRIME = 0.0d0 + !APRIM + IF ( HcoState%Grid%YEDGE%Val(I,J) .lt. 0 ) AGE = Inst%SAGE*24.0d0 + IF ( HcoState%Grid%YEDGE%Val(I,J) .ge. 0 ) AGE = Inst%NAGE*24.0d0 + APRIM = (1.038d0+0.03758d0*AGE-0.00014349d0*AGE**2d0 & + + (1.911315d-7*AGE**3d0) )**(-1d0) + ! Threshold wind speed [m/s] + UT = 6.975d0 + 0.0033d0 * (TEMP - 273.15d0 + 27.27d0 )**2.0d0 + !IF (W10M > UT) THEN + ! add RH<100 too + + IF (W10M > UT .and. RH_ICE<100d0) THEN + QBSALT = 0.385d0*(1.0d0-Ut/W10M)**2.59d0/FRICTVEL + QB0 = 0.385d0*(1d0-6.975d0/W10M)**2.59d0/FRICTVEL + ! Snow sublimation rate [kg/m2/s] Equation 1 in Yang et al. (2008) + ! The constant 1.1574d-5 converts mm/day column integrated sublimation rate to kg m-2 s-1 + + QS = 1.1574d-5*APRIM*QSPRIME*QBSALT/QB0 + ELSE + QS = 0d0 + ENDIF + !set up the snow salinity + IF ( HcoState%Grid%YEDGE%Val(I,J) .lt. 0 ) SLNT_FYI = Inst%SSLNT_FYI + IF ( HcoState%Grid%YEDGE%Val(I,J) .lt. 0 ) SLNT_MYI = Inst%SSLNT_MYI + IF ( HcoState%Grid%YEDGE%Val(I,J) .ge. 0 ) SLNT_FYI = Inst%NSLNT_FYI + IF ( HcoState%Grid%YEDGE%Val(I,J) .ge. 0 ) SLNT_MYI = Inst%NSLNT_MYI + ! Sea ice fraction that is first year + FRFIRST = ExtState%FRSEAICE%Arr%Val(I,J) - MULTI(I,J) + IF ( FRFIRST < 0d0 ) FRFIRST = 0d0 + ! Apply FYI salinity to FYI seaice fraction and MYI salinity to MYI fraction + !SLNT = SLNT_FYI * FRFIRST + SLNT_MYI * MULTI(I,J) + ! Assume MYI salinity is 50% of FYI + !SLNT = SLNT * FRFIRST + SLNT * 0.5 * MULTI(I,J) + ! Convert snow sublimation rate to sea salt production rate [kg/m2/s] + ! Calculate it separately for FYI and MYI, scaled by their respective sea ice fraction + QSNOWICE_FYI = QS * SLNT_FYI * FRFIRST / 1000d0 + QSNOWICE_MYI = QS * SLNT_MYI * MULTI(I,J) / 1000d0 + !print *, 'Bettyhere is QSNOW ',QSNOWICE_FYI,QSNOWICE_MYI, I,J + ELSE + QSNOWICE_FYI = 0.0d0 + QSNOWICE_MYI = 0.0d0 + ENDIF ENDIF + ! End of added blowing snow section - ! Skip boxes that are not at least 50% water - IF ( FOCEAN > 0.5d0 ) THEN - - ! Wind speed at 10 m altitude [m/s] - W10M = SQRT( ExtState%U10M%Arr%Val(I,J)**2 & - + ExtState%V10M%Arr%Val(I,J)**2 ) - - ! Sea surface temperature in Celcius - SST = ExtState%TSKIN%Arr%Val(I,J) - 273.15d0 - - ! Limit SST to 0-30C range - SST = MAX( SST , 0d0 ) ! limit to 0C - SST = MIN( SST , 30d0 ) ! limit to 30C - - ! Empirical SST scaling factor (jaegle 5/11/11) - SCALE = 0.329d0 + 0.0904d0*SST - & - 0.00717d0*SST**2d0 + 0.000207d0*SST**3d0 + !----------------------------------------------------------------- !--------------------------------------------------------------- ! Partition TOMAS_Jeagle emissions w/in the boundary layer @@ -230,11 +434,13 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) else - dfo=0.d0 - endif + dfo=0.d0 + endif - dfo=dfo*Inst%DRFAC(k)*BETHA !hemco units???? jkodros - dfo=dfo*focean*SCALE + dfo=dfo*Inst%DRFAC(k)*BETHA !hemco units???? jkodros (?m-2 s-1) + + !dfo=dfo*focean*SCALE ! scale now includes ocean fraction - remove this line + dfo=dfo*SCALE ! note: scale now includes ocean fraction ! Loop thru the boundary layer DO L = 1, HcoState%Nz @@ -246,7 +452,44 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) IF ( FEMIS > 0d0 ) THEN ! Number - NUMBER = dfo * FEMIS + !betty NUMBER = dfo * FEMIS ! need to move after add blowing snow + + ! update seasalt from blowing snow - huang 1/4/18 + IF (( Inst%EmitSnowSS )) THEN + IF ( HcoState%Grid%YEDGE%Val(I,J) .lt. 0 ) THEN + ! Southern Hemisphere + + !if (K > 3 ) THEN + NUMBER = FEMIS* (dfo + & + (( QSNOWICE_FYI * SUM( Inst%F_DN_S_FYI(:,K) ) + & + QSNOWICE_MYI * SUM( Inst%F_DN_S_MYI(:,K) ) ) * DDSNOW)) + + !NUMBER = FEMIS* (dfo) + FEMIS* ( & + ! (( QSNOWICE_FYI * SUM( Inst%F_DI_S_FYI(:,K) ) + & + ! QSNOWICE_MYI * SUM( Inst%F_DI_S_MYI(:,K) ) ) * DDSNOW)) & + ! / (Inst%SS_DEN( K ) * ((Inst%TOMAS_DBIN(k))**3.d0) & + ! * HcoState%Phys%PI/6.0d0) + !else + ! NUMBER = FEMIS * dfo + !endif + !IF (NUMBER2 .GT. 5.e5) THEN + ! print *, 'First', NUMBER2,NUMBER, I,J,K,L,FEMIS + !ENDIF + + !IF (NUMBER .GT. 5.e5) THEN + !!print*,'Betty here are emissions',NUMBER,dfo*FEMIS,K,I,J,L, & + ! Inst%F_DN_S_FYI(:,K), DDSNOW, SUM( Inst%F_DN_S_FYI(:,K) ), & + ! QSNOWICE_FYI, QSNOWICE_MYI, FEMIS + !ENDIF + ELSE + ! Northern Hemisphere + NUMBER = FEMIS* (dfo + & + (( QSNOWICE_FYI * SUM( Inst%F_DN_N_FYI(:,K) ) + & + QSNOWICE_MYI * SUM( Inst%F_DN_N_MYI(:,K) ) ) * DDSNOW)) + ENDIF + ELSE ! only open ocean sea salt + NUMBER = dfo * FEMIS + ENDIF ! Mass MASS = NUMBER & @@ -260,15 +503,11 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) ENDIF ENDDO ENDDO - ELSE - Inst%TC1(I,J,:,:) = 0d0 - Inst%TC2(I,J,:,:) = 0d0 - ENDIF ENDDO ENDDO !### Debug - !print*, 'JACK SEASALT EMISSIONS AT 50, 10,7: ', TC2(ii,jj,1,7) + !print*, 'Aerosol mass AT 50, 10, 7: ', Inst%TC2(ii,jj,1,7) !print*, 'BINS: ', HcoState%MicroPhys%nBins ! Loop over # of microphysics bins @@ -283,7 +522,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) ! Get the proper species name IF ( K<10 ) THEN - WRITE(SpcName,'(A2,I1)') 'NK', K + WRITE(SpcName,'(A3,I1)') 'NK0', K ELSE WRITE(SpcName,'(A2,I2)') 'NK', K ENDIF @@ -292,8 +531,8 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Run( ExtState, HcoState, RC ) HcoID = HCO_GetHcoID( TRIM(SpcName), HcoState ) !### Debug - !print*, 'JACK SEASALT EMISSIONS AT 50, 10,: ', TC1(ii,jj,1,k) - !print*, 'JACK HCO ID: ', HcoID + print*, 'Aerosol number AT 50, 10,: ', Inst%TC1(ii,jj,1,k) + print*, 'HCO ID: ', K, SpcName, HcoID ! Add number to the HEMCO data structure CALL HCO_EmisAdd( HcoState, Inst%TC1(:,:,:,K), HcoID, RC) @@ -335,6 +574,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) USE HCO_State_Mod, ONLY : HCO_GetHcoID USE HCO_STATE_MOD, ONLY : HCO_GetExtHcoID USE HCO_ExtList_Mod, ONLY : GetExtNr + USE HCO_ExtList_Mod, ONLY : GetExtOpt ! ! !INPUT PARAMETERS: ! @@ -369,6 +609,16 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) ! Pointers TYPE(MyInst), POINTER :: Inst + ! Local variables for blowing snow + INTEGER :: K, ND, IH !IH for different hemisphere + REAL*8 :: D_SNOW, D_DRY + REAL*8, PARAMETER :: A_SALT = 2.0d0 !from Mann et al. 2000 + REAL*8, PARAMETER :: B_SALT = 37.5d0 !in um + !REAL*8, PARAMETER :: DDSNOW = 2.0d0 !in um for snow particle interval + REAL*8, PARAMETER :: DDSNOW = 0.1d0 !in um for snow particle interval + REAL*8, PARAMETER :: RHONACL = 2160.0d0 !kg/m3 + REAL*8, PARAMETER :: RHOICE = 900.0d0 !kg/m3 + !================================================================= ! HCOX_TOMAS_Jeagle_Init begins here! !================================================================= @@ -392,9 +642,9 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) msg = & 'Using HEMCO extension: TOMAS_Jeagle (sea salt emissions for TOMAS)' IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN - CALL HCO_Msg( HcoState%Config%Err, sep1='-' ) ! with separator + CALL HCO_Msg( HcoState%Config%Err, msg, sep1='-' ) ! with separator ELSE - CALL HCO_Msg( msg, verb=.TRUE. ) ! w/o separator + CALL HCO_Msg( msg, verb=.TRUE. ) ! w/o separator ENDIF ENDIF @@ -412,6 +662,79 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) ! Get species IDs and settings ! ---------------------------------------------------------------------- + ! betty + ! print *,'Betty start of options' + + ! fix scaling factor over cold water SST (<5 degC) + CALL GetExtOpt ( HcoState%Config, Inst%ExtNr, 'Reduce SS cold water', & + OptValBool=Inst%ColdSST, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Add a SSA source from blowing snow (by J. Huang) + CALL GetExtOpt ( HcoState%Config, Inst%ExtNr, 'Blowing Snow SS', & + OptValBool=Inst%EmitSnowSS, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + + ! Whether or not differentiate snow salinity on FYI and MYI (by J. Huang) + !CALL GetExtOpt ( HcoState%Config, Inst%ExtNr, 'Diff salinity on ice', & + ! OptValBool=Inst%FYIsnow, RC=RC ) + !IF ( RC /= HCO_SUCCESS ) RETURN + ! Add snow salinity (NH and SH), snow age and number of particles + ! per snowflake as external factor from configuration file + IF ( Inst%EmitSnowSS ) THEN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'NH FYI snow salinity', & + OptValDp=Inst%NSLNT_FYI, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'NH MYI snow salinity', & + OptValDp=Inst%NSLNT_MYI, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'SH FYI snow salinity', & + OptValDp=Inst%SSLNT_FYI, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'SH MYI snow salinity', & + OptValDp=Inst%SSLNT_MYI, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'NH snow age', & + OptValDp=Inst%NAGE, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'SH snow age', & + OptValDp=Inst%SAGE, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + CALL GetExtOpt( HcoState%Config, Inst%ExtNr, 'N per snowflake', & + OptValDp=Inst%NP, RC=RC ) + IF ( RC /= HCO_SUCCESS ) RETURN + ELSE + Inst%NSLNT_FYI = 0.1d0 ! default value 0.1 psu for NH FYI snow + Inst%NSLNT_MYI = 0.05d0 ! default value 0.05 psu for NH MYI snow + Inst%SSLNT_FYI = 0.03d0 ! default value 0.03 psu for SH FYI snow + Inst%SSLNT_FYI = 0.015d0 ! default value 0.015 psu for SH MYI snow + Inst%NAGE = 3.0d0 ! default value 3 days snow age in NH + Inst%SAGE = 1.5d0 ! default value 1.5 days snow age in SH + Inst%NP = 5.0d0 ! default value of 5 particles per snowflake + ENDIF + + ! Verbose mode + IF ( HcoState%amIRoot ) THEN + MSG = 'Use sea salt aerosol emissions (extension module)' + CALL HCO_MSG(HcoState%Config%Err,MSG, SEP1='-' ) + IF ( Inst%EmitSnowSS ) THEN + WRITE(MSG,*) ' - Arctic Snow Salinity on FYI (psu): ', Inst%NSLNT_FYI + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Arctic Snow Salinity on MYI (psu): ', Inst%NSLNT_MYI + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Antarctic Snow Salinity on FYI (psu): ', Inst%SSLNT_FYI + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Antarctic Snow Salinity on FYI (psu): ', Inst%SSLNT_MYI + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Arctic Snow age (days): ', Inst%NAGE + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Antarctic Snow age(days): ', Inst%SAGE + CALL HCO_MSG(HcoState%Config%Err,MSG) + WRITE(MSG,*) ' - Number of particle per snowflake: ', Inst%NP + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDIF + ENDIF + ! Get HEMCO species IDs CALL HCO_GetExtHcoID( HcoState, Inst%ExtNr, Inst%HcoIDs, SpcNames, & nSpc, RC ) @@ -424,6 +747,89 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) CALL HCO_ERROR(MSG, RC ) RETURN ENDIF + IF ( HcoState%amIRoot ) THEN + MSG = 'Use the following species (Name: HcoID):' + CALL HCO_MSG(HcoState%Config%Err,MSG) + DO N = 1, nSpc + WRITE(MSG,*) TRIM(SpcNames(N)), ':', Inst%HcoIDs(N) + CALL HCO_MSG(HcoState%Config%Err,MSG) + ENDDO + ENDIF + + ALLOCATE ( Inst%SS_DEN ( HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate SS_DEN' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%SS_DEN = 2200.d0 + + ! Allocate for blowing snow simulation + IF ( Inst%EmitSnowSS ) THEN + ALLOCATE ( Inst%F_DI_N_FYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DI_N_FYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DI_N_FYI = 0.0_sp + + ALLOCATE ( Inst%F_DI_N_MYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DI_N_MYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DI_N_MYI = 0.0_sp + + ALLOCATE ( Inst%F_DN_N_FYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DN_N_FYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DN_N_FYI = 0.0_sp + + ALLOCATE ( Inst%F_DN_N_MYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DN_N_MYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DN_N_MYI = 0.0_sp + + ALLOCATE ( Inst%F_DI_S_FYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DI_S_FYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DI_S_FYI = 0.0_sp + + ALLOCATE ( Inst%F_DI_S_MYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DI_S_MYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DI_S_MYI = 0.0_sp + + ALLOCATE ( Inst%F_DN_S_FYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DN_S_FYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DN_S_FYI = 0.0_sp + + ALLOCATE ( Inst%F_DN_S_MYI( NR_MAX, HcoState%MicroPhys%nBins ), STAT=AS ) + IF ( AS/=0 ) THEN + MSG = 'Cannot allocate F_DN_S_MYI' + CALL HCO_ERROR(MSG, RC ) + RETURN + ENDIF + Inst%F_DN_S_MYI = 0.0_sp + ENDIF ! Allocate TOMAS_DBIN ALLOCATE ( Inst%TOMAS_DBIN( HcoState%MicroPhys%nBins ), STAT=RC ) @@ -441,7 +847,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) RETURN ENDIF - ! JKODROS - ALLOCATE TC1 and TC2 + ! Allocate TC1 ALLOCATE ( Inst%TC1( HcoState%NX, HcoState%NY,& HcoState%NZ, HcoState%MicroPhys%nBins ), STAT=RC ) IF ( RC /= HCO_SUCCESS ) THEN @@ -452,7 +858,7 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) Inst%TC1 = 0d0 ENDIF - ! JKODROS - ALLOCATE TC1 and TC2 + ! Allocate TC2 ALLOCATE ( Inst%TC2( HcoState%NX, HcoState%NY,& HcoState%NZ, HcoState%MicroPhys%nBins ), STAT=RC ) IF ( RC /= HCO_SUCCESS ) THEN @@ -463,6 +869,107 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) Inst%TC2 = 0d0 ENDIF + !size bins for blowing snow - Huang 6/12/20 + IF ( Inst%EmitSnowSS ) THEN + + DO K = 1, HcoState%MicroPhys%nBins + ! TOMAS dry bin limits (radius) + + R0=0.5d6*((HcoState%MicroPhys%BinBound(K)*6.0d0)/ & + (Inst%SS_DEN( K )*HcoState%Phys%PI))**( 1d0 / 3d0 ) ! in um + R1=0.5d6*((HcoState%MicroPhys%BinBound(K+1)*6.0d0)/ & + (Inst%SS_DEN( K )*HcoState%Phys%PI))**( 1d0 / 3d0 ) ! in um + + !-------------- Define size distribution --------------------- + ! for northern hemisphere FYI + D_SNOW = 1.0d0 + DO ND = 1, NR_MAX + D_DRY = ( Inst%NSLNT_FYI * RHOICE / (1000.d0 & + * Inst%NP * RHONACL ) )**( 1d0 / 3d0 ) * D_SNOW + !print*,'Here is Ddry ',D_DRY, R0*2d0 , R1*2d0 + IF (D_DRY .ge. R0*2d0 .and. D_DRY .le. R1*2d0 ) THEN + !---------------------------------------------------------- + ! NOTES: + ! For size distribution + ! define the two-parameter gamma probability density funtion here + ! Yang et al 2008 eq (6) + !---------------------------------------------------------- + ! Midpoint of IRth bin + Inst%F_DI_N_FYI(ND, K) = EXP( - D_SNOW / B_SALT ) & + * D_SNOW**( A_SALT - 1.d0 ) & + / ( B_SALT**A_SALT * GAMMA( A_SALT ) ) + !print*,'Here is F_DI ', ND, K, Inst%F_DI_N_FYI(ND, K) + ELSE + Inst%F_DI_N_FYI(ND, K) = 0d0 + ENDIF + Inst%F_DN_N_FYI(ND, K) = Inst%F_DI_N_FYI(ND,K) / (4d0/3d0 * HcoState%Phys%PI & + * 1.d-18 * Inst%SS_DEN( K ) * (D_DRY/2d0)**3) + D_SNOW = D_SNOW + DDSNOW + + !print*,'Here is F_DN ', ND, K, Inst%F_DN_N_FYI(ND, K) + + ENDDO + + ! for northern hemisphere MYI + D_SNOW = 1.0d0 + DO ND = 1, NR_MAX + D_DRY = ( Inst%NSLNT_MYI * RHOICE / (1000.d0 & + * Inst%NP * RHONACL ) )**( 1d0 / 3d0 ) * D_SNOW + + IF (D_DRY .ge. R0*2d0 .and. D_DRY .le. R1*2d0 ) THEN + ! Midpoint of IRth bin + Inst%F_DI_N_MYI(ND, K) = EXP( - D_SNOW / B_SALT ) & + * D_SNOW**( A_SALT - 1.d0 ) & + / ( B_SALT**A_SALT * GAMMA( A_SALT ) ) + ELSE + Inst%F_DI_N_MYI(ND, K) = 0d0 + ENDIF + Inst%F_DN_N_MYI(ND, K) = Inst%F_DI_N_MYI(ND,K) / (4d0/3d0 * HcoState%Phys%PI & + * 1.d-18 * Inst%SS_DEN( K ) * (D_DRY/2d0)**3) + + D_SNOW = D_SNOW + DDSNOW + ENDDO + + + ! for southern hemisphere FYI + D_SNOW = 1.0d0 + DO ND = 1, NR_MAX + D_DRY = ( Inst%SSLNT_FYI * RHOICE / (1000.d0 & + * Inst%NP * RHONACL ) )**( 1d0 / 3d0 ) * D_SNOW + + IF (D_DRY .ge. R0*2d0 .and. D_DRY .le. R1*2d0 ) THEN + ! Midpoint of IRth bin + Inst%F_DI_S_FYI(ND, K) = EXP( - D_SNOW / B_SALT ) & + * D_SNOW**( A_SALT - 1.d0 ) & + / ( B_SALT**A_SALT * GAMMA( A_SALT ) ) + ELSE + Inst%F_DI_S_FYI(ND, K) = 0d0 + ENDIF + Inst%F_DN_S_FYI(ND, K) = Inst%F_DI_S_FYI(ND,K)/ (4d0/3d0 * HcoState%Phys%PI & + * 1.d-18 * Inst%SS_DEN( K ) * (D_DRY/2d0)**3) + D_SNOW = D_SNOW + DDSNOW + ENDDO + + ! for southern hemisphere MYI + D_SNOW = 1.0d0 + DO ND = 1, NR_MAX + D_DRY = ( Inst%SSLNT_MYI * RHOICE / (1000.d0 & + * Inst%NP * RHONACL ) )**( 1d0 / 3d0 ) * D_SNOW + IF (D_DRY .ge. R0*2d0 .and. D_DRY .le. R1*2d0 ) THEN + ! Midpoint of IRth bin + Inst%F_DI_S_MYI(ND, K) = EXP( - D_SNOW / B_SALT ) & + * D_SNOW**( A_SALT - 1.d0 ) & + / ( B_SALT**A_SALT * GAMMA( A_SALT ) ) + ELSE + Inst%F_DI_S_MYI(ND, K) = 0d0 + ENDIF + Inst%F_DN_S_MYI(ND, K) = Inst%F_DI_S_MYI(ND,K)/ (4d0/3d0 * HcoState%Phys%PI & + * 1.d-18 * Inst%SS_DEN( K ) * (D_DRY/2d0)**3) + D_SNOW = D_SNOW + DDSNOW + ENDDO + ENDDO !K + ENDIF + ! ----- IMPORTANT BINS ONLY CORRECTLY SET UP FOR TOMAS 15 PLEASE ADJUST OTHERS -jkodros (7/21/15) ! ---- 6/24/16 - JKodros - I have updated the DRFAC. They should (hopefully) be ! ---- correct now. DRFAC is the bin width (radius not diameter) for DRY SS @@ -575,16 +1082,21 @@ SUBROUTINE HCOX_TOMAS_Jeagle_Init( HcoState, ExtName, ExtState, RC ) ! Activate met fields ExtState%FRLAND%DoUse = .TRUE. - ExtState%FRLANDIC%DoUse = .TRUE. ExtState%FROCEAN%DoUse = .TRUE. ExtState%FRSEAICE%DoUse = .TRUE. - ExtState%FRLAKE%DoUse = .TRUE. ExtState%TSKIN%DoUse = .TRUE. ExtState%U10M%DoUse = .TRUE. ExtState%V10M%DoUse = .TRUE. ExtState%FRAC_OF_PBL%DoUse = .TRUE. ExtState%FRCLND%DoUse = .TRUE. + ! for blowing snow + IF ( Inst%EmitSnowSS ) THEN + ExtState%USTAR%DoUse = .TRUE. + ExtState%T2M%DoUse = .TRUE. + ExtState%QV2M%DoUse = .TRUE. + ENDIF + !======================================================================= ! Leave w/ success !======================================================================= @@ -726,6 +1238,9 @@ SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC ) ! Generic instance initialization ! ---------------------------------------------------------------- + !betty + !print*, 'Betty In instance' + ! Initialize Inst => NULL() @@ -742,6 +1257,17 @@ SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC ) Inst%Instance = nnInst + 1 Inst%ExtNr = ExtNr + !Init values + Inst%ColdSST = .FALSE. + Inst%EmitSnowSS = .FALSE. + Inst%NSLNT_FYI = 0.0 + Inst%NSLNT_MYI = 0.0 + Inst%SSLNT_FYI = 0.0 + Inst%SSLNT_MYI = 0.0 + Inst%NAGE = 0.0 + Inst%SAGE = 0.0 + Inst%NP = 1.0 + ! Attach to instance list Inst%NextInst => AllInst AllInst => Inst @@ -829,6 +1355,56 @@ SUBROUTINE InstRemove ( Instance ) DEALLOCATE( Inst%HcoIDs ) ENDIF + IF ( ASSOCIATED( Inst%SS_DEN ) ) THEN + DEALLOCATE( Inst%SS_DEN ) + ENDIF + Inst%SS_DEN => NULL() + + IF ( ASSOCIATED( Inst%F_DI_N_FYI ) ) THEN + DEALLOCATE( Inst%F_DI_N_FYI ) + ENDIF + Inst%F_DI_N_FYI => NULL() + + IF ( ASSOCIATED( Inst%F_DI_N_MYI ) ) THEN + DEALLOCATE( Inst%F_DI_N_MYI ) + ENDIF + Inst%F_DI_N_MYI => NULL() + + IF ( ASSOCIATED( Inst%F_DI_S_FYI ) ) THEN + DEALLOCATE( Inst%F_DI_S_FYI ) + ENDIF + Inst%F_DI_S_FYI => NULL() + + IF ( ASSOCIATED( Inst%F_DI_S_MYI ) ) THEN + DEALLOCATE( Inst%F_DI_S_MYI ) + ENDIF + Inst%F_DI_S_MYI => NULL() + + IF ( ASSOCIATED( Inst%F_DN_N_FYI ) ) THEN + DEALLOCATE( Inst%F_DN_N_FYI ) + ENDIF + Inst%F_DN_N_FYI => NULL() + + IF ( ASSOCIATED( Inst%F_DN_N_MYI ) ) THEN + DEALLOCATE( Inst%F_DN_N_MYI ) + ENDIF + Inst%F_DN_N_MYI => NULL() + + IF ( ASSOCIATED( Inst%F_DN_S_FYI ) ) THEN + DEALLOCATE( Inst%F_DN_S_FYI ) + ENDIF + Inst%F_DN_S_FYI => NULL() + + IF ( ASSOCIATED( Inst%F_DN_S_MYI ) ) THEN + DEALLOCATE( Inst%F_DN_S_MYI ) + ENDIF + Inst%F_DN_S_MYI => NULL() + + IF ( ASSOCIATED( Inst%MULTIICE ) ) THEN + DEALLOCATE( Inst%MULTIICE ) + ENDIF + Inst%MULTIICE => NULL() + !--------------------------------------------------------------------- ! Pop off instance from list !---------------------------------------------------------------------