!--------------------------------------------------------------------------------- ! this file is based on structure of gw_drag from WACCMX 3.5.48 !--------------------------------------------------------------------------------- module ionosphere !--------------------------------------------------------------------------------- ! Purpose: ! ! Module to compute the relevant physics and electrodynamics needed to improve and ! add a realistic simulation of the ionosphere. Main modules are initialization, ! ion/electron temperature calculation, ambipolar diffusion, .... ! ! Authors: Joe McInerney/Hanli Liu/Art Richmond ! !--------------------------------------------------------------------------------- use shr_kind_mod, only : r8 => shr_kind_r8 ! Real kind to declare variables use params, only : iulog implicit none save private ! Make default type private to the module ! ! PUBLIC: interfaces ! public ionos_intr ! interface to actual ionosphere simulation routines public tridag ! Generic tridiagonal solver routine ! ! PUBLIC variables ! ! ! Define public variables ! !------------------------------------------------------------------------ ! PRIVATE: Rest of the data and interfaces are private to this module !------------------------------------------------------------------------ real(r8), parameter :: kboltz_ev = 8.617E-5_r8 ! Boltzmann constant (eV/K) real(r8), parameter :: temax = 7.0E3_r8 ! maximum electron temperature (K) real(r8), parameter :: dayOPFlux = 2.0E8_r8 ! Daytime O+ flux at upper boundary ( real(r8), parameter :: nightOPFlux = -2.0E8_r8 ! Nighttime O+ flux at upper boundary ( !------------------------------------------------------------------------------ ! magnetic grid dimensions (assumed resolution of 2deg) !------------------------------------------------------------------------------ integer, parameter :: nmlat = 90 ! mlat integer, parameter :: nmlon = 180 ! mlon integer, parameter :: pcols = 1 ! number of columns integer, parameter :: pver = 97 ! number of columns integer, parameter :: pverp = 98 ! number of columns integer :: ionBot ! bottom of ionosphere calculations integer :: ionBotP ! bottom of ionosphere calculations integer :: oPBot ! bottom of O+ transport calculations integer :: oPBotP ! bottom of O+ transport calculations real(r8) :: calDay ! current calendar day real(r8) :: rads2Degs ! radians to degrees real(r8) :: degs2Rads ! degrees to radians real(r8) :: f107 ! 10.7 cm solar flux type ionos_state real(r8), dimension(pcols) :: cosZenAngR ! cosine of zenith angle (radians) real(r8), dimension(pcols) :: zenAngD ! zenith angle (degrees) real(r8), dimension(pcols,pver) :: bNorth3d ! northward component of magnetic field units? real(r8), dimension(pcols,pver) :: bEast3d ! eastward component of magnetic field real(r8), dimension(pcols,pver) :: bDown3d ! downward component of magnetic field real(r8), dimension(pcols,pver) :: sumIonPRates ! Sum of ionization rates for O+,O2+,N+,N2+,NO+ (s-2 cm-3) real(r8), dimension(pcols,pver) :: ue2d ! horizontal x(eastward) component of ExB drift with added vertical dimension real(r8), dimension(pcols,pver) :: ve2d ! horizontal y(northward) component of ExB drift with added vertical dimension real(r8), dimension(pcols,pver) :: we2d ! vertical z component of ExB drift with added vertical dimension - midpoints real(r8), dimension(pcols,pverp) :: wei2d ! vertical z component of ExB drift with added vertical dimension - interfaces real(r8), dimension(pcols,pverp) :: omegai ! vertical velocity on interface levels (Pa/s) real(r8), dimension(pcols,pver) :: dipMag ! dip angle for each column (radians) real(r8), dimension(pcols,pver) :: dipMagD ! dip angle for each column (degrees) real(r8), dimension(pcols,pverp) :: tNInt ! Interface Temperature (K) real(r8), dimension(pcols,pver) :: ndensN2 ! N2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3) real(r8), dimension(pcols,pver) :: ndensNO ! NO number density (cm-3) real(r8), dimension(pcols,pver) :: ndensN1 ! N number density (cm-3) real(r8), dimension(pcols,pver) :: ndensE ! E electron number density (cm-3) real(r8), dimension(pcols,pver) :: ndensOp ! O plus number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2p ! O2 plus ion number density (cm-3) real(r8), dimension(pcols,pver) :: ndensNOp ! NO plus ion number density (cm-3) real(r8), dimension(pcols,pver) :: sourceg4 ! g4 source term for electron/ion temperature update real(r8), dimension(pcols,pverp) :: rairvi ! Constituent dependent gas constant on interface levels real(r8), dimension(pcols,pver) :: rairv ! Constituent dependent gas constant on interface levels real(r8), dimension(pcols,pver) :: mbarv ! Constituent dependent gas constant on interface levels end type ionos_state contains !============================================================================== subroutine ionos_intr(ncol,oPBot,ionBot) !----------------------------------------------------------------------- ! Interface for improved ionosphere simulation !----------------------------------------------------------------------- ! !------------------------------Arguments-------------------------------- implicit none type(ionos_state) :: istate ! ionosphere state structure integer, intent(in) :: ionBot ! bottom of ionosphere calculations integer, intent(in) :: oPBot ! bottom of O+ transport calculations integer, intent(in) :: ncol ! number of atmospheric columns !---------------------------Local storage------------------------------- integer :: lchnk ! chunk identifier !------------------------------------------------------------ ! Initialize data needed in the ionosphere calculations !------------------------------------------------------------ write(iulog,*)'Calling ionos_datinit ' call ionos_datinit(ncol, istate, ionBot) !------------------------------- ! Get electron temperature !------------------------------- write(iulog,*)'Calling ionos_tempei ' call ionos_tempei(ncol, istate, ionBot) !------------------------------- ! Get vertical O+ transport !------------------------------- ! write(iulog,*)'Calling ionos_optransvert ' ! call ionos_optransvert(ncol, istate, oPBot) !--------------------------------------- ! Update heating of neutral atmosphere !-------------------------------------- return end subroutine ionos_intr !=============================================================================== subroutine ionos_datinit(ncol, istate, ionBot) !------------------------------------------------------------------------------ ! Time independent initialization for ionosphere simulation called in phys_init ! of physpkg module which is called in cam_comp module !------------------------------------------------------------------------------ use read_ncfile, only: bEast,bNorth,bDown,zenAngR,tN,mmrO1,mmrO2,mmrNO,omega,pMid,tE,tI,ndensE use params, only: rair,mbar,pi,kboltz,avogad implicit none type(ionos_state), intent(inout) :: istate ! ionosphere state structure integer, intent(in) :: ncol ! Number of columns in current chunk integer, intent(in) :: ionBot ! bottom of ionosphere calculations !---------------------------Local storage------------------------------- integer :: iVer ! Counter for vertical loops integer :: iCol ! Counter for column loops integer :: nstep ! current timestep number real(r8), dimension(ncol) :: geoLatR ! Latitude (radians) Make ncol because zenith aurora are ncol real(r8), dimension(ncol) :: geoLonR ! Longitude (radians) real(r8), dimension(pcols,pverp) :: pInt ! Interface pressure (Pa) real(r8), dimension(pcols,pverp) :: tNInt ! Interface Temperture (K) real(r8), dimension(ncol) :: cosZenAngR ! cosine of zenith angle (radians) real(r8), dimension(pcols,pver) :: bNorth3d ! northward component of magnetic field units? real(r8), dimension(pcols,pver) :: bEast3d ! eastward component of magnetic field real(r8), dimension(pcols,pver) :: bDown3d ! downward component of magnetic field real(r8), dimension(pcols,pver) :: ue2d ! horizontal x(eastward) component of ExB drift with added vertical dimension real(r8), dimension(pcols,pver) :: ve2d ! horizontal y(northward) component of ExB drift with added vertical dimension real(r8), dimension(pcols,pver) :: we2d ! vertical z component of ExB drift with added vertical dimension - midpoints real(r8), dimension(pcols,pver) :: sourceR ! R term of source g4 calculation real(r8), dimension(pcols,pver) :: sourceEff ! Efficiency term of source g4 calculation real(r8), dimension(pcols,pverp) :: wei2d ! vertical z component of ExB drift with added vertical dimension - interfaces real(r8), dimension(pcols,pverp) :: omegai ! vertical velocity on interface levels (Pa/s) real(r8), dimension(pcols,pverp) :: rairvi ! Constituent dependent gas constant real(r8), dimension(pcols,pver) :: dipMag ! dip angle for each column (radians) real(r8), dimension(pcols,pver) :: dipMagD ! dip angle for each column (degrees) real(r8) :: zenAngD ! Zenith angle (degrees) real(r8) :: rmassN2 ! N2 molecular weight kg/kmol real(r8) :: rmassO2 ! O2 molecular weight kg/kmol real(r8) :: rmassO1 ! O atomic weight kg/kmol real(r8) :: rmassNO ! NO molecular weight kg/kmol real(r8) :: rmassH ! H atomic weight kg/kmol real(r8) :: rmassN1 ! N atomic weight kg/kmol real(r8) :: rmassE ! Electron mass kg/kmol real(r8) :: rmassOp ! O plus ion molecular weight kg/kmol real(r8) :: rmassO2p ! O2 plus ion molecular weight kg/kmol real(r8) :: rmassNOp ! NO plus ion molecular weight kg/kmol real(r8), dimension(pcols,pver) :: mmrN2 ! N2 mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: mmrN1 ! N mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: mmrE ! E electron mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: mmrO2p ! O2 plus ion mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: mmrNOp ! NO plus ion mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: ndensN2 ! N2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3) real(r8), dimension(pcols,pver) :: ndensNO ! NO number density (cm-3) real(r8), dimension(pcols,pver) :: ndensN1 ! N number density (cm-3) ! real(r8), dimension(pcols,pver) :: ndensE ! E electron number density (cm-3) real(r8), dimension(pcols,pver) :: ndensOp ! O plus number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2p ! O2 plus ion number density (cm-3) real(r8), dimension(pcols,pver) :: ndensNOp ! NO plus ion number density (cm-3) real(r8), dimension(pcols,pver) :: sumIonPRates ! Sum of ionization rates for O+,O2+,N+,N2+,NO+ (s-2 cm-3) real(r8), dimension(:,:), pointer :: aurIPRateSum ! Auroral ion production sum for O2+,O+,N2+ (s-1 cm-3 from module mo_aurora) real(r8), dimension(pcols,pver) :: sourceg4 ! g4 source term for electron/ion temperature update logical, dimension(ncol) :: do_aurora ! Logical to test if aurora ion production rates calculated in mo_aurora (set to zero) logical :: xo2_slvd, xo1_slvd, o2p_slvd, op_slvd, nop_slvd, elec_slvd, xn_slvd, xno_slvd, xh_slvd ! Logicals for short lived species status real(r8), dimension(pver) :: aurIPRateSumVert ! Temporary WACCMX vertical auroral ion production rates real(r8), dimension(pver) :: sumIonPRatesVert ! Temporary vertical ion production rates from WACCM-X integer :: iVerWX real(r8), dimension(81) :: sumIonPRatesVertWX ! Input WACCMX vertical ion production rates data sumIonPRatesVertWX/ & 16.5765747548862, 27.4290745912036, 45.3255917144736, & 74.8891312623343, 123.760902197506, 204.571788067206, & 338.053031965589, 557.653386039374, 915.485592788164, & 1487.61186755924, 2372.30478446464, 3664.13926483388, & 5368.60033719279, 7230.88243545343, 8419.73873706100, & 6438.96635132295, 5324.29947674881, 4453.48900014561, & 3880.06555700992, 2761.23170301566, 2351.91384616588, & 2654.65367734872, 3170.60548773902, 2953.67926160213, & 1486.88451766190, 267.959117588397, 9.63693080703941, & 0.387022309460580, 0.193178223497298, 7.941694775215291E-002, & 1.522856580132134E-002, 2.471485235558893E-003, 1.265506017012632E-003, & 4.454530050343265E-004, 5.938093202320889E-005, 1.908870045659309E-006, & 1.051451121559231E-008, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000 / real(r8), dimension(81) :: aurIPRateSumVertWX ! Input vertical auroral ion production rates from WACCM-X data aurIPRateSumVertWX/ & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000, & 0.000000000000000E+000, 0.000000000000000E+000, 0.000000000000000E+000/ real(r8), dimension(81) :: pMidWX ! Temporary pressure from WACCM-X data pMidWX/ & 3.296541015681223E-007, 5.435077232800263E-007, 8.960927343535380E-007, & 1.477407135477792E-006, 2.435832543085364E-006, 4.016008881691840E-006, & 6.621279194093812E-006, 1.091664372706028E-005, 1.799850252045966E-005, & 2.967451362143404E-005, 4.892500126983817E-005, 8.066368937972138E-005, & 1.329919389978525E-004, 2.192666362574653E-004, 3.615095631957054E-004, & 5.960300000000000E-004, 9.826900000000000E-004, 1.620185000000000E-003, & 2.671225000000000E-003, 4.404100000000000E-003, 7.261275000000000E-003, & 1.197190000000000E-002, 1.973800000000000E-002, 3.254225000000000E-002, & 5.365325000000000E-002, 8.846025000000000E-002, 0.145845750000000, & 0.240457500000000, 0.397825000000000, 0.655682550000000, & 1.08138255000000, 1.78980000000000, 2.95577500000000, & 4.87307500000000, 7.99107500000000, 12.8273250000000, & 19.8120000000000, 29.2025000000000, 41.0167500000000, & 55.3470000000000, 73.0480000000000, 95.5947500000000, & 124.479500000000, 161.285000000000, 207.932500000000, & 266.742500000000, 340.487500000000, 432.457500000000, & 546.540000000000, 687.285000000000, 859.972500000000, & 1070.70500000000, 1326.47500000000, 1635.17500000000, & 2005.67500000000, 2447.90000000000, 2972.80000000000, & 3592.32500000000, 4319.37500000000, 5167.75000000000, & 6152.05000000000, 7287.45000000000, 8570.99855327839, & 10068.2524365878, 11855.3448670211, 13957.7658625785, & 16431.1566273078, 19340.9711462438, 22764.2188933347, & 26791.4946874561, 31529.3717896922, 37103.2378947354, & 43660.6021568008, 51375.0037784712, 60450.6012591619, & 70151.2448927345, 79311.1631922018, 87317.6428692450, & 93614.1191900039, 97735.9170217135, 99952.3736593074/ !-------------------------------------------------------------------------------- tNInt(:,:) = 0._r8 bNorth3d(:,:) = 0._r8 bEast3d(:,:) = 0._r8 bDown3d(:,:) = 0._r8 ue2d(:,:) = 0._r8 ve2d(:,:) = 0._r8 we2d(:,:) = 0._r8 wei2d(:,:) = 0._r8 omegai(:,:) = 0._r8 rairvi(:,:) = 0._r8 dipMag(:,:) = 0._r8 dipMagD(:,:) = 0._r8 sumIonPRates(:,:) = 0._r8 ndensN2(:,:) = 0._r8 ndensO2(:,:) = 0._r8 ndensO1(:,:) = 0._r8 sourceR(:,:) = 0._r8 sourceEff(:,:) = 0._r8 sourceg4(:,:) = 0._r8 !------------------------------------------------------------------------------------------------------- ! Convert production rates from WACCM-X pressure grid (81 levels) to TIME-GCM pressure grid (97 levels) !------------------------------------------------------------------------------------------------------- do iCol = 1, ncol do iVer = 1, pver do IVerWX = 2, 81 if ( pMid(iCol,iVer) <= pMidWX(1) ) then sumIonPRatesVert(iVer) = sumIonPRatesVertWX(1) aurIPRateSumVert(iVer) = aurIPRateSumVertWX(1) write(iulog,*) 'pMid(iCol,iVer),pMidWX(iVerWX) ',pMid(iCol,iVer),pMidWX(iVerWX) write(iulog,*) 'sumIonPRatesVert(iVer) ',sumIonPRatesVert(iVer) exit endif if ( pMid(iCol,iVer) > pMidWX(iVerWX-1) .and. pMid(iCol,iVer) < pMidWX(iVerWX) ) then sumIonPRatesVert(iVer) = (sumIonPRatesVertWX(iVerWX-1) + sumIonPRatesVertWX(iVerWX)) / 2. aurIPRateSumVert(iVer) = (aurIPRateSumVertWX(iVerWX-1) + aurIPRateSumVertWX(iVerWX)) / 2. write(iulog,*) 'pMid(iCol,iVer),pMidWX(iVerWX) ',pMid(iCol,iVer),pMidWX(iVerWX) write(iulog,*) 'sumIonPRatesVertWX(iVerWX-1),sumIonPRatesVertWX(iVerWX),sumIonPRatesVert(iVer) ',sumIonPRatesVertWX(iVerWX-1),sumIonPRatesVertWX(iVerWX),sumIonPRatesVert(iVer) exit endif enddo enddo enddo !------------------------------------------------------------------------------------------------------ ! Set the bottom of the ionosphere calculations at around 0.5 hectopascals(millibars). hypm is in ! Pascals. Since the model vertical range goes from top down this can be used as the counter of the ! number of levels in the range from the top down to the ionBot level. !------------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------------------------------ ! Set the bottom of the O+ transport calculations at around .001 Pascals or .00001 hectopascals(millibars) or ~130km. !------------------------------------------------------------------------------------------------------------------------ do iVer = 1, pver ! if (hypm(iVer) < 50._r8) oPBot = iVer !~60km ! if (hypm(iVer) < 0.01_r8) oPBot = iVer ! ~100 km !S if (hypm(iVer) < 0.001_r8) oPBot = iVer ! ~130 km ! if (hypm(iVer) < 0.0005_r8) oPBot = iVer ! ~150 km enddo !------------------------------------------------------------------- ! Set radians to degrees variable and get zenith angle in degrees !------------------------------------------------------------------- rads2Degs = 180._r8/pi degs2Rads = pi/180._r8 zenAngD = zenAngR * rads2Degs !------------------------------------------------------------------------------------- ! Calculate neutral temperature on interface levels. tN vertical dimension is pver !------------------------------------------------------------------------------------- do iVer = 2, pver do iCol = 1, ncol tNInt(iCol,iVer) = 0.5_r8 * tN(iCol,iVer) + 0.5_r8 * tN(iCol,iVer-1) enddo enddo do iCol = 1, ncol tNInt(iCol,1) = 1.5_r8 * tNInt(iCol,2) - 0.5_r8 * tNInt(iCol,3) enddo do iCol = 1, ncol tNInt(iCol,pverp) = 1.5_r8 * tNInt(iCol,pver) - 0.5_r8 * tNInt(iCol,pver-1) enddo istate%tNInt(1:ncol,1:pverp) = tNInt(1:ncol,1:pverp) cosZenAngR(:) = COS(zenAngD * degs2Rads) istate%cosZenAngR(1:ncol) = cosZenAngR(1:ncol) istate%zenAngD(1:ncol) = zenAngD istate%bNorth3d(1:ncol,1:pver) = bNorth istate%bEast3d(1:ncol,1:pver) = bEast istate%bDown3d(1:ncol,1:pver) = bDown !--------------------------------------------------------------------------------- ! Get ion ExB drift from physics buffer (they were defined by exbdrift module) ! Ion drifts are 2d output arrays, i.e., no vertical dimension but 1d here (pcols) !--------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Form 2D drifts needed for this module. Vertical level goes to ionBotP since ! needed below to get vertical drift on interface levels !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- ! Need vertical ExB drift on interface levels !------------------------------------------------------------------------------- do iVer = 2, pver do iCol = 1, ncol wei2d(iCol,iVer) = 0.5_r8 * (we2d(iCol,iVer-1) + we2d(iCol,iVer)) enddo enddo do iCol = 1, ncol wei2d(iCol,1) = 1.5_r8 * wei2d(iCol,2) - 0.5_r8 * wei2d(iCol,3) enddo do iCol = 1, ncol wei2d(iCol,pverp) = 1.5_r8 * wei2d(iCol,pver) - 0.5_r8 * wei2d(iCol,pver-1) enddo istate%wei2d(1:ncol,1:pverp) = wei2d(1:ncol,1:pverp) !-------------------------------------------------------------------------------------- ! Need to get vertical velocity on interface levels handling top level specifically !-------------------------------------------------------------------------------------- omegai(1:ncol,1:pver) = omega(1:ncol,1:pver) omegai(1:ncol,pverp) = 1.5_r8 * omegai(1:ncol,pver) - 0.5_r8 * omegai(1:ncol,pver-1) istate%omegai(1:ncol,1:pverp) = omegai(1:ncol,1:pverp) !------------------------------------------------------------------------ ! Get constituent dependent gas constant and derive on interface levels !------------------------------------------------------------------------ ! istate%rairvi(1:ncol,1:pverp) = rair !------------------------------------------------------------------------------- ! Need to get dip angle from magnetic field components !------------------------------------------------------------------------------- do iVer = 1, pver do iCol = 1, ncol dipMag(iCol,iVer) = ATAN(bDown / SQRT(bNorth**2 + bEast**2)) if (dipMag(iCol,iVer) < 0.17_r8 .and. dipMag(iCol,iVer) > 0._r8 ) dipMag(iCol,iVer) = 0.17_r8 if (dipMag(iCol,iVer) > -0.17_r8 .and. dipMag(iCol,iVer) < 0._r8 ) dipMag(iCol,iVer) = 0.17_r8 dipMagD(iCol,iVer) = dipMag(iCol,iVer) * rads2Degs enddo enddo istate%dipMag(1:ncol,1:pver) = dipMag(1:ncol,1:pver) istate%dipMagD(1:ncol,1:pver) = dipMagD(1:ncol,1:pver) !---------------------------------------------------------------------------------------- ! Set flag whether each constituent is short-lived !---------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------- ! Get molecular weights using constituent indices (kg/kmole) except N2 which is hard-wired !---------------------------------------------------------------------------------------- rmassO2 = 32._r8 rmassO1 = 16._r8 rmassNO = 30._r8 rmassH = 1._r8 rmassN1 = 14._r8 rmassE = 1.E-04_r8 rmassOP = rmassO1 rmassO2P = rmassO2 rmassNOP = rmassNO rmassN2 = 28._r8 do iVer = 1,pver do iCol = 1, ncol istate%mbarv(iCol,iVer) = 1._r8 / (mmrO1(iCol,iVer)/rmassO1 + mmrO2(iCol,iVer)/rmassO2 + (1._r8-mmrO1(iCol,iVer)-mmrO2(iCol,iVer))/rmassN2) enddo enddo istate%rairv(1:ncol,1:pver) = kboltz * avogad / istate%mbarv(1:ncol,1:pver) do iVer = 2,pver istate%rairvi(1:ncol,iVer) = 0.5_r8 * (istate%rairv(1:ncol,iVer-1)+istate%rairv(1:ncol,iVer)) enddo istate%rairvi(1:ncol,1) = 1.5_r8*istate%rairvi(1:ncol,2)-0.5_r8*istate%rairvi(1:ncol,3) istate%rairvi(1:ncol,pverp) = rair write(iulog,*)'ionos_datinit before density calcs' !----------------------------------------------------------------------------------------------------- ! Get mass mixing ratios from state structure q array (kg/kg) except N2 which is calculated from O,O2 ! or unless constituent is short-lived !----------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------------------- ! Need to get number density (cgs units) from mass mixing ratio. mbarv is kg/mole, same as rmass units ! kg/kg * (kg/mole)/(kg/mole) * (Pa or N/m*m)/((Joules/K or N*m/K) * (K)) = m-3 * 1E-06 = cm-3 !-------------------------------------------------------------------------------------------------------------- do iVer = 1, pver do iCol = 1, ncol ndensO2(iCol,iVer) = mmrO2(iCol,iVer) * mbar / rmassO2 * & pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ndensO1(iCol,iVer) = mmrO1(iCol,iVer) * mbar / rmassO1 * & pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ndensNO(iCol,iVer) = mmrNO(iCol,iVer) * mbar / rmassNO * & pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ! ndensN1(iCol,iVer) = mmrN1(iCol,iVer) * mbar / rmassN1 * & ! pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ! ndensE(iCol,iVer) = mmrE(iCol,iVer) * mbar / rmassE * & ! pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ndensN2(iCol,iVer) = (1._r8-mmrO1(iCol,iVer)-mmrO2(iCol,iVer)) * mbar / rmassN2 * & pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 enddo enddo write(iulog,*)'ionos_datinit after density calcs' istate%ndensO2(1:ncol,1:pver) = ndensO2(1:ncol,1:pver) istate%ndensO1(1:ncol,1:pver) = ndensO1(1:ncol,1:pver) istate%ndensNO(1:ncol,1:pver) = ndensNO(1:ncol,1:pver) ! istate%ndensN1(1:ncol,1:pver) = ndensN1(1:ncol,1:pver) istate%ndensN2(1:ncol,1:pver) = ndensN2(1:ncol,1:pver) ! istate%ndensE(1:ncol,1:pver) = ndensE(1:ncol,1:pver) !------------------------------------------------------------------------------- ! Write output fields to history file !------------------------------------------------------------------------------- !------------------------------------------------------------------------------------ ! Get ionization rates from physics buffer which were calculated in mo_jeuv and ! mo_jshort modules. Rates array dimensions are pcols, pver, nIonRates. Units s-1 !------------------------------------------------------------------------------------ !---------------------------------------------------------------------------------------------- ! Need to convert these ionization rates to ion production rates by multiplying number density ! of neutral species appropriate from reactions in mo_jeuv(jeuv) and mo_jshort(jshort)(for NO) !---------------------------------------------------------------------------------------------- !---------------------------------------------- ! Sum ion production rates all reactions !---------------------------------------------- !------------------------------------------------------------------------------------------- ! Get aurora ion production rate sum from physics buffer which were calculated in mo_aurora ! module. Rate array dimensions are pcols, pver. Units s-1 cm-3 !------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------- ! Check latitudes, and set aurora ion production rates for all columns and levels ! to zero if all below 32.5 deg since not calculated in mo_aurora for those latitudes ! (same as criteria in mo_aurora) !------------------------------------------------------------------------------------- WRITE(iulog,*) 'sumIonPRatesVert(1:ionBot) ', sumIonPRatesVert(1:ionBot) WRITE(iulog,*) 'aurIPRateSumVert(1:ionBot)', aurIPRateSumVert(1:ionBot) !------------------------------------------------------------------------------------------------- ! Calculate electron heating rate which is source in electron/ion temperature derivation !------------------------------------------------------------------------------------------------- do iVer = 1, ionBot do iCol = 1, ncol sourceR(iCol,iVer) = LOG( ndensE(iCol,iVer) / (ndensO2(iCol,iVer) + ndensN2(iCol,iVer) + & 0.1_r8 * ndensO1(iCol,iVer)) ) sourceEff(iCol,iVer) = EXP( -(12.75_r8 + 6.941_r8 * sourceR(iCol,iVer) + 1.166_r8 * sourceR(iCol,iVer)**2 + & 0.08043_r8 * sourceR(iCol,iVer)**3 + 0.001996_r8 * sourceR(iCol,iVer)**4) ) !S!------------------------------------------------------------------------------- !S! Calculate g4 source term for electron temperature update !S!------------------------------------------------------------------------------- ! sourceg4(iCol,iVer) = (sumIonPRates(iCol,iVer) + aurIPRateSum(iCol,iVer)) * sourceEff(iCol,iVer) WRITE(iulog,*) 'sumIonPRatesVert(iVer) ', sumIonPRatesVert(iVer) WRITE(iulog,*) 'aurIPRateSumVert(iVer) ', aurIPRateSumVert(iVer) WRITE(iulog,*) 'sourceEFf(iCol,iVer) ', sourceEff(iCol,iVer) sourceg4(iCol,iVer) = (sumIonPRatesVert(iVer) + aurIPRateSumVert(iVer)) * sourceEff(iCol,iVer) enddo enddo WRITE(iulog,*) 'sourceR(1,1:ionBot) ',sourceR(1,1:ionBot) WRITE(iulog,*) 'sourceg4(1,1:ionBot) ',sourceg4(1,1:ionBot) istate%sourceg4(1:ncol,1:pver) = sourceg4(1:ncol,1:pver) return end subroutine ionos_datinit ! !=============================================================================== subroutine ionos_tempei(ncol, istate, ionBot) !----------------------------------------------------------------------- ! Routine to compute the electron and ion temperature needed for ! ambipolar diffusion calculations. !----------------------------------------------------------------------- !T use cam_history, only : addfld, add_default, phys_decomp !T use phys_grid, only : get_lat_p, get_lon_p, get_rlat_p,get_rlon_p !T use physconst, only : gravit ! Gravity (m/s2) !T use cam_control_mod, only : nsrest ! Variable to determine if this is an initial run or a restart/branch !T use time_manager, only : get_nstep ! Routine to get current time step !T use physconst, only : rairv, mbarv ! Constituent dependent rair and mbar use read_ncfile, only : tE,tI,omega,tN,uN,vN,ed1_geo,ed2_geo,bMag,bEast,bNorth,ndensOp,ndensO2p,ndensNOp,ndensE,geoMagLat,geoMagLon use geogrid, only : pMid,pInt use params, only : iulog,rair,gravity,avogad,pi,ztodt,mbar,kboltz implicit none !------------------------------Arguments-------------------------------- !T type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer !T type(physics_state), intent(in) :: state ! physics state structure !T type(physics_ptend), intent(inout) :: ptend ! parameterization tendency structure type(ionos_state), intent(inout) :: istate ! ionosphere state structure !T real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) !T integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: ionBot ! bottom of ionosphere calculations !---------------------------Local storage------------------------------- integer :: ionBotP ! bottom of ionosphere calculations plus one more level integer :: i,k ! loop indexes integer :: iVer ! Counter for vertical loops integer :: iCol ! Counter for column loops integer :: iter ! Counter for iteration loop integer :: nstep ! current timestep number integer :: indxTe ! pbuf index for electron temperature integer :: indxTi ! pbuf index for ion temperature integer, parameter :: maxIter = 1 ! maximum number of iterations to solve for electron/ion temperature real(r8), parameter :: Kec1 = 7.5E5_r8 ! c1 constant for calculation of electron conductivity(Ke) real(r8), parameter :: Kec2 = 3.22E4_r8 ! c2 constant for calculation of electron conductivity(Ke) real(r8), parameter :: stepweight = 1.0_r8 ! weight of previous and current times step for diagonals real(r8) :: wrk1 ! 2/3/kboltz_ev real(r8) :: sqrt2 real(r8) :: lossc5 ! c5 prime of Lc(eN2) component of loss term real(r8) :: lossc7 ! c7 of Lc(eO2) component of loss term equation real(r8) :: lossc9 ! c9 of Lc(eO) component of loss term equation real(r8) :: lossc13 ! c13 of Lc(eO)f component of loss term equation real(r8) :: FeDB ! B term of electron heat flux of UB real(r8) :: FeD ! Day time flux real(r8) :: FeN ! Night time flux real(r8) :: f1Ted1 ! d1 of f1(Te) calculation used to get electron conductivity real(r8) :: f1Ted2 ! d2 of f1(Te) calculation used to get electron conductivity real(r8) :: f1Ted3 ! d3 of f1(Te) calculation used to get electron conductivity real(r8) :: f1Te !T real(r8), dimension(:,:), pointer :: tE ! Pointer to electron temperature in pbuf (K) !T real(r8), dimension(:,:), pointer :: ti ! Pointer to ion temperature in pbuf (K) !T !T real(r8), dimension(pcols,pver) :: pMid ! Midpoint pressure (Pa) !T real(r8), dimension(pcols,pver) :: tN ! Neutral temperature (K) real(r8), dimension(pcols,pver) :: tE0 ! Electron temperature from last time step (K) real(r8), dimension(pcols,pver) :: tEPrevI ! Electron temperature from previous iteration (K) !T real(r8), dimension(pcols,pverp) :: pInt ! Interface pressure (Pa) real(r8), dimension(pcols,pverp) :: tNInt ! Interface Temperture (K) real(r8), dimension(pcols,pverp) :: rairvi ! Constituent dependent gas constant on interface levels real(r8), dimension(pcols,pver) :: rairv ! Constituent dependent gas constant on midpoint levels real(r8), dimension(pcols,pver) :: ndensN2 ! N2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3) real(r8), dimension(pcols,pver) :: ndensNO ! NO number density (cm-3) real(r8), dimension(pcols,pver) :: ndensN1 ! N number density (cm-3) !T real(r8), dimension(pcols,pver) :: ndensE ! E electron number density (cm-3) !T real(r8), dimension(pcols,pver) :: ndensOp ! O plus number density (cm-3) !T real(r8), dimension(pcols,pver) :: ndensO2p ! O2 plus ion number density (cm-3) !T real(r8), dimension(pcols,pver) :: ndensNOp ! NO plus ion number density (cm-3) real(r8), dimension(pcols,pver) :: dipMag ! dip angle for each column (radians) real(r8), dimension(pcols,pver) :: dipMagD ! dip angle for each column (degrees) real(r8), dimension(ncol) :: zenAngD ! zenith angle (degrees) real(r8), dimension(pcols) :: Feub ! electron heat flux at upper boundary real(r8), dimension(pcols,pver) :: sqrtTE ! Square root of electron temperature real(r8), dimension(pcols,pver) :: sqrtTN ! Square root of electron temperature real(r8), dimension(pver) :: Ke ! electron conductivity real(r8), dimension(pverp) :: Kei ! electron conductivity interface levels real(r8), dimension(pcols,pver) :: lossc4p ! c4 prime of Lc(eN2) component of loss term real(r8), dimension(pcols,pver) :: lossceN2 ! Lc(eN2) component of loss term equation real(r8), dimension(pcols,pver) :: lossc6p ! c6 prime of Lc(eO2) component of loss term equation real(r8), dimension(pcols,pver) :: lossceO2 ! Lc(eO2) component of loss term equation real(r8), dimension(pcols,pver) :: lossc8p ! c8 prime of Lc(eO) component of loss term equation real(r8), dimension(pcols,pver) :: lossceO1 ! Lc(eO) component of loss term equation real(r8), dimension(pcols,pver) :: lossc10p ! c10 prime of Lc(eN2) component of loss term equation real(r8), dimension(pcols,pver) :: losscA ! A of Lc(eN2)v component of loss term equation real(r8), dimension(pcols,pver) :: tENDiff ! Difference between electron and neutral temperatures real(r8), dimension(pcols,pver) :: lossceN2v ! Lc(eN2)v component of loss term equation real(r8), dimension(pcols,pver) :: lossc11p ! c11 prime of Lc(eO2)v component of loss term equation real(r8), dimension(pcols,pver) :: lossceO2v ! Lc(eO2)v component of loss term equation real(r8), dimension(pcols,pver) :: lossc12p ! c12 prime of Lc(eO)f component of loss term equation real(r8), dimension(pcols,pver) :: lossceOf ! Lc(eO)f component of loss term equation real(r8), dimension(pcols,pver) :: lossc14p ! c14 prime of Lc(eO)1D component of loss term equation real(r8), dimension(pcols,pver) :: losscf2d ! d of f2 of Lc(eO)1D component of loss term equation real(r8), dimension(pcols,pver) :: losscf2 ! f2 of Lc(eO)1D component of loss term equation real(r8), dimension(pcols,pver) :: losscf3 ! f3 of Lc(eO)1D component of loss term equation real(r8), dimension(pcols,pver) :: lossceO1D ! Lc(eO)1D component of loss term equation real(r8), dimension(pcols,pver) :: lossc15p ! c15 prime of Lc(eN2)Rot component of loss term equation real(r8), dimension(pcols,pver) :: lossceN2Rot ! Lc(eN2)Rot component of loss term equation real(r8), dimension(pcols,pver) :: lossc16p ! c16 prime of Lc(eO2)Rot component of loss term equation real(r8), dimension(pcols,pver) :: lossceO2Rot ! Lc(eO2)Rot component of loss term equation real(r8), dimension(pcols,pver) :: lossc3p ! c3 prime of Lc(ei) component of loss term equation real(r8), dimension(pcols,pver) :: losscei ! Lc(ei) component of loss term equation real(r8), dimension(pcols,pver) :: losscin ! ion-neutral heating coeff. real(r8), dimension(pcols,pver) :: lossg3 ! g3 loss term for Te tendency real(r8), dimension(pcols,pver) :: delTEN ! Difference between electron and neutral temperatures from production/loss real(r8), dimension(pcols,pverp) :: delZi ! Delta z: interfaces real(r8), dimension(pcols,pver) :: delZ ! Delta z: midpoints real(r8), dimension(pcols,pver) :: sourceg4 ! g4 source term for electron/ion temperature update real(r8), dimension(pcols,pver) :: qjoule ! joule heating real(r8), dimension(pcols,pver) :: qen ! electron-neutral heating real(r8), dimension(pcols,pver) :: qei ! electron-ion Coulomb heating real(r8), dimension(pcols,pver) :: rho ! mass density real(r8), dimension(pcols,pver) :: wrk2 real(r8), dimension(ionBot) :: subdiag ! subdiagonal values for Te tendency solving real(r8), dimension(ionBot) :: superdiag ! superdiagonal values for Te tendency solving real(r8), dimension(ionBot) :: diag ! diagonal values for Te tendency solving real(r8), dimension(ionBot) :: rHSInit ! initial RHS of electron temperature update real(r8), dimension(ionBot) :: rHSH ! h for RHS of electron temperature update real(r8), dimension(ionBot) :: rHS ! RHS of electron temperature update real(r8), dimension(ionBot) :: tETemp ! temporary electron temperature array for input to tridag logical, dimension(pcols) :: colConv ! flag for column converging = 1 if converged otherwise = 0 logical :: converged !Flag for convergence in electron temperature calculation iteration loop !----------------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------------- ! Initialize arrays to zero !--------------------------------------------------------------------------------------------------------- ! pMid(:,:) = 0._r8 ! tN(:,:) = 0._r8 ! pInt(:,:) = 0._r8 tNInt(:,:) = 0._r8 rairvi(:,:) = 0._r8 ndensO2(:,:) = 0._r8 ndensO1(:,:) = 0._r8 ndensNO(:,:) = 0._r8 ndensN1(:,:) = 0._r8 ! ndensE(:,:) = 0._r8 ! ndensOp(:,:) = 0._r8 ! ndensO2p(:,:) = 0._r8 ! ndensNOp(:,:) = 0._r8 ndensN2(:,:) = 0._r8 zenAngD(:) = 0._r8 sqrtTE(:,:) = 0._r8 sqrtTN(:,:) = 0._r8 Ke(:) = 0._r8 Kei(:) = 0._r8 lossc4p(:,:) = 0._r8 lossceN2(:,:) = 0._r8 lossc6p(:,:) = 0._r8 lossc6p(:,:) = 0._r8 lossceO2(:,:) = 0._r8 lossc8p(:,:) = 0._r8 lossceO1(:,:) = 0._r8 lossc10p(:,:) = 0._r8 losscA(:,:) = 0._r8 tENDiff(:,:) = 0._r8 lossceN2v(:,:) = 0._r8 lossc11p(:,:) = 0._r8 lossceO2v(:,:) = 0._r8 lossc12p(:,:) = 0._r8 lossceOf(:,:) = 0._r8 lossc14p(:,:) = 0._r8 losscf2d(:,:) = 0._r8 losscf2(:,:) = 0._r8 losscf3(:,:) = 0._r8 lossceO1D(:,:) = 0._r8 lossc15p(:,:) = 0._r8 lossceN2Rot(:,:) = 0._r8 lossc16p(:,:) = 0._r8 lossceO2Rot(:,:) = 0._r8 lossc3p(:,:) = 0._r8 losscei(:,:) = 0._r8 losscin(:,:) = 0._r8 lossg3(:,:) = 0._r8 delTEN(:,:) = 0._r8 delZi(:,:) = 0._r8 delZ(:,:) = 0._r8 subDiag(:) = 0._r8 superDiag(:) = 0._r8 diag(:) = 0._r8 sourceg4(:,:) = 0._r8 dipMag(:,:) = 0._r8 dipMagD(:,:) = 0._r8 rHSInit(:) = 0._r8 rHSH(:) = 0._r8 rHS(:) = 0._r8 teTemp(:) = 0._r8 qjoule(:,:) = 0._r8 qei(:,:) = 0._r8 qen(:,:) = 0._r8 rho(:,:) = 0._r8 colConv(:) = .false. !------------------------------------------- ! Calculate some commonly used variables !------------------------------------------- wrk1 = 2._r8 / 3._r8/ kboltz_ev sqrt2 = SQRT(2._r8) ionBotP = ionBot + 1 ndensO1(1:ncol,1:pver) = istate%ndensO1(1:ncol,1:pver) ndensN2(1:ncol,1:pver) = istate%ndensN2(1:ncol,1:pver) ndensO2(1:ncol,1:pver) = istate%ndensO2(1:ncol,1:pver) ndensNO(1:ncol,1:pver) = istate%ndensNO(1:ncol,1:pver) tNInt(1:ncol,1:pverp) = istate%tNInt(1:ncol,1:pverp) rairvi(1:ncol,1:pverp) = istate%rairvi(1:ncol,1:pverp) rairv(1:ncol,1:pver) = istate%rairv(1:ncol,1:pver) !------------------------------------------------------------------------------------------------------- ! Need to get midpoint and interface pressure and neutral temperature from state structure (ncol,ionBot) !------------------------------------------------------------------------------------------------------- !T pMid(1:ncol,1:pver) = state%pmid(1:ncol,1:pver) !T tN(1:ncol,1:pver) = state%t(1:ncol,1:pver) rho(1:ncol,1:pver) = pMid(1:ncol,1:pver)/rairv(1:ncol,1:pver)/tN(1:ncol,1:pver) * 1.E-3_r8 ! convert to g/cm3 !T qjoule(1:ncol,1:ionBot) = ptend%s(1:ncol,1:ionBot) * 6.24E15_r8 ! convert from J/kg/s to ev/g/s !T pInt(1:ncol,1:pverp) = state%pint(1:ncol,1:pverp) !T tNInt(1:ncol,1:pverp) = istate%tNInt(1:ncol,1:pverp) !T rairvi(1:ncol,1:pverp) = istate%rairvi(1:ncol,1:pverp) sqrtTN(1:ncol,1:pver) = SQRT(tN(1:ncol,1:pver)) !---------------------------------------------------------------- ! Get variables needed from the ionosphere state structure !---------------------------------------------------------------- ndensO2(1:ncol,1:pver) = istate%ndensO2(1:ncol,1:pver) ndensO1(1:ncol,1:pver) = istate%ndensO1(1:ncol,1:pver) ndensNO(1:ncol,1:pver) = istate%ndensNO(1:ncol,1:pver) !???? ndensN1(1:ncol,1:pver) = istate%ndensN1(1:ncol,1:pver) ndensN1(1:ncol,1:pver) = 1.E6_r8 !T ndensE(1:ncol,1:pver) = istate%ndensE(1:ncol,1:pver) !T ndensOp(1:ncol,1:pver) = istate%ndensOp(1:ncol,1:pver) !T ndensO2p(1:ncol,1:pver) = istate%ndensO2p(1:ncol,1:pver) !T ndensNOp(1:ncol,1:pver) = istate%ndensNOp(1:ncol,1:pver) ndensN2(1:ncol,1:pver) = istate%ndensN2(1:ncol,1:pver) sourceg4(1:ncol,1:pver) = istate%sourceg4(1:ncol,1:pver) dipMag(1:ncol,1:pver) = istate%dipMag(1:ncol,1:pver) dipMagD(1:ncol,1:pver) = istate%dipMagD(1:ncol,1:pver) zenAngD(1:ncol) = istate%zenAngD(1:ncol) !------------------------------------------------------------------------------------------------------------------- ! Get electron temperature from physics buffer and if this is the first time calculated then initialize to neutral ! temperature. nsrest=0 means this is an initial run and nstep=0 means this is the first time step. !------------------------------------------------------------------------------------------------------------------- !T indxTe = pbuf_get_index( 'ElecTemp' ) !T call pbuf_get_field(pbuf, indxTe, tE) !T !T indxTi = pbuf_get_index( 'IonTemp' ) !T call pbuf_get_field(pbuf, indxTi, ti) !T !T nstep = get_nstep() nstep = 1 !T if(nsrest == 0 .and. nstep == 0) then !T !T ti(1:ncol,1:pver) = tN(1:ncol,1:pver) !T tE(1:ncol,1:pver) = tN(1:ncol,1:pver) !T !T else WRITE(iulog,*) 'tE(1,1:ionBot) 1', tE(1,1:ionBot) !T tE(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),tE(1:ncol,1:pver)) !T tE(1:ncol,1:pver) = MIN(temax,tE(1:ncol,1:pver)) !T !T ti(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),ti(1:ncol,1:pver)) !T ti(1:ncol,1:pver) = MIN(ti(1:ncol,1:pver),tE(1:ncol,1:pver)) !T !T tE(1:ncol,ionBotP:pver) = tN(1:ncol,ionBotP:pver) !T ti(1:ncol,ionBotP:pver) = tN(1:ncol,ionBotP:pver) WRITE(iulog,*) 'tE(1,1:ionBot) 2', tE(1,1:ionBot) !T endif tE0(1:ncol,1:pver) = tE(1:ncol,1:pver) wrk2(1:ncol,1:ionBot) = ndensE(1:ncol,1:ionBot)/wrk1/(SIN(dipMag(1:ncol,1:ionBot)))**2._r8 !----------------------------------------------------------------------------- ! Get constant terms needed for loss term g3 for electron temperature update !----------------------------------------------------------------------------- lossc5 = 1.21E-4_r8 lossc7 = 3.6E-2_r8 lossc9 = 5.7E-4_r8 lossc13 = 7.E-5_r8 WRITE(iulog,*) 'tE(1,1:ionBot) 3', tE(1,1:ionBot) !----------------------------------------------------------------------------- ! Get terms needed for loss term g3 for electron temperature update which do ! not need to be updated in iteration loop. !----------------------------------------------------------------------------- do iCol = 1, ncol if (.not. colConv(iCol)) then do iVer = 1, ionBot lossc4p(iCol,iVer) = 1.77E-19_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) lossc6p(iCol,iVer) = 1.21E-18_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) lossc8p(iCol,iVer) = 7.9E-19_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) lossc10p(iCol,iVer) = 1.3E-4_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) lossc11p(iCol,iVer) = 3.125E-21_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) lossc12p(iCol,iVer) = 3.4E-12_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) lossc14p(iCol,iVer) = 1.57E-12_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) lossc15p(iCol,iVer) = 2.9E-14_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) lossc16p(iCol,iVer) = 6.9E-14_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) lossc3p(iCol,iVer) = 3.2E-8_r8 * 15._r8 * (ndensOP(iCol,iVer) + & 0.5_r8 * ndensO2P(iCol,iVer) + 0.53_r8 * ndensNOP(iCol,iVer)) * ndensE(iCol,iVer) losscin(iCol,iVer) = (6.6e-14_r8*ndensN2(iCol,iVer) + 5.8e-14_r8*ndensO2(iCol,iVer) & + 0.21e-14_r8*ndensO1(iCol,iVer)*sqrt2*sqrtTN(iCol,iVer))*ndensOP(iCol,iVer) & +(5.9e-14_r8*ndensN2(iCol,iVer) + 5.45e-14_r8*ndensO2(iCol,iVer) & + 4.5e-14_r8*ndensO1(iCol,iVer))*ndensOP(iCol,iVer) & +(5.8e-14_r8*ndensN2(iCol,iVer) + 0.14e-14_r8*ndensO2(iCol,iVer)*sqrtTN(iCol,iVer) & + 4.4e-14_r8*ndensO1(iCol,iVer)) * ndensO2P(iCol,iVer) enddo !iVer loop !---------------------------------------------------------------------------------- ! Calculate upper boundary heat flux !---------------------------------------------------------------------------------- if (ABS(dipMagD(iCol,1)) < 60.0_r8) FeDB = 0.5_r8 * & (1._r8 + SIN(pi * (ABS(dipMagD(iCol,1)) - 30.0_r8) /60.0_r8)) if (ABS(dipMagD(iCol,1)) >= 60.0_r8) FeDB = 1._r8 FeD = -4.5E-7_r8 * f107 * 0.5_r8 * FeDB FeN = .2_r8 * FeD !--------------------------------------------------- ! Set upper boundary condition for right hand side !--------------------------------------------------- if (zenAngD(iCol) <= 80.0_r8) Feub(iCol) = FeD if (zenAngD(iCol) > 80.0_r8 .AND. zenAngD(iCol) < 100.0_r8) Feub(iCol) = 0.5_r8 * (FeD + FeN) & + 0.5_r8 * (FeD - FeN) * & COS(pi * ((zenAngD(iCol) - 80.0_r8) / 20.0_r8)) if (zenAngD(iCol) >= 100.0_r8) Feub(iCol) = FeN !------------------------------------------------------------------------------------------ ! Calculate thickness terms for vertical derivative !------------------------------------------------------------------------------------------ do iVer = 1, ionBot delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rairv(iCol,iVer) * tN(iCol,iVer) / pMid(iCol,iVer) / gravity enddo do iVer = 2, ionBotP ! Assuming ionBotP < pverp delZi(iCol,iVer) = (pMid(iCol,iVer) - pMid(iCol,iVer-1)) * rairvi(iCol,iVer) * & tNInt(iCol,iVer) / pInt(iCol,iVer) / gravity enddo delZi(iCol,1) = 1.5_r8*delZi(iCol,2) - .5_r8*delZi(iCol,3) !---------------------------------------------------------- ! Convert delZ variables from meters to centimeters !---------------------------------------------------------- delZi(iCol,1:ionBotP) = delZi(iCol,1:ionBotP)*100._r8 delZ(iCol,1:ionBot) = delZ(iCol,1:ionBot)*100._r8 endif ! Column not converged enddo !iCol loop !------------------------------------------------------------------------------------------------------- ! Iterate to calculate new electron temperature. ! Time splitting is used: first solve the heating/cooling equation, then solve the diffusion equations. ! Also, set convergence flag to false and iterate until true or 6 iterations, whichever comes first !------------------------------------------------------------------------------------------------------- converged = .false. iter = 0 do while (.not. converged .and. iter < maxIter) !-------------------------------------------------------------------------------------------------------- ! Increment iteration loop counter and save electron temperature from previous iteration for convergence ! test at end of interation loop. Also, take square root of electron temperature to be used later !-------------------------------------------------------------------------------------------------------- iter = iter + 1 tEPrevI(1:ncol,1:ionBot) = tE(1:ncol,1:ionBot) sqrtTE(1:ncol,1:ionBot) = SQRT(tE(1:ncol,1:ionBot)) WRITE(iulog,*) 'tE(1,1:ionBot) 4', tE(1,1:ionBot) !-------------------------------------------------------------------------------------------------------- ! Loop over columns then vertical levels and call tridiagonal solver for each column to get electron ! temperature !-------------------------------------------------------------------------------------------------------- do iCol = 1, ncol if (.not. colConv(iCol)) then do iVer = 1, ionBot !----------------------------------------------------------------------------- ! Get loss term g3 for electron temperature update. Need to calculate ! constituent dependent loss terms which make up g3 !----------------------------------------------------------------------------- lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer) lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iCol,iVer)) * sqrtTE(iCol,iVer) lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iCol,iVer) if (tE(iCol,iVer) < 1000.0_r8) losscA(iCol,iVer) = 5.71E-8_r8 * EXP(-3352.6_r8 / tE(iCol,iVer)) if (tE(iCol,iVer) >= 1000.0_r8 .AND. tE(iCol,iVer) <= 2000.0_r8) & losscA(iCol,iVer) = 2.0E-7_r8 * EXP(-4605.2_r8 / tE(iCol,iVer)) if (tE(iCol,iVer) > 2000.0_r8) losscA(iCol,iVer) = 2.53E-6_r8 * sqrtTE(iCol,iVer) * & EXP(-17620._r8 / tE(iCol,iVer)) tENDiff(iCol,iVer) = tE(iCol,iVer) - tN(iCol,iVer) if (ABS(tENDiff(iCol,iVer)) < 0.1_r8) tENDiff(iCol,iVer) = 0.1_r8 lossceN2v(iCol,iVer) = lossc10p(iCol,iVer) * losscA(iCol,iVer) * & (1._r8 - EXP(3200._r8 * (1._r8 / tE(iCol,iVer) - 1._r8 / tN(iCol,iVer)))) / tENDiff(iCol,iVer) lossceO2v(iCol,iVer) = lossc11p(iCol,iVer) * tE(iCol,iVer)**2 lossceOf(iCol,iVer) = lossc12p(iCol,iVer) * (1._r8 - lossc13 * tE(iCol,iVer)) * & (0.4_r8 + 150._r8 / tE(iCol,iVer)) / tN(iCol,iVer) losscf2d(iCol,iVer) = 2.4E+4_r8 + 0.3_r8 * (tE(iCol,iVer) - 1500._r8) + & 1.947E-5_r8 * (tE(iCol,iVer) - 1500._r8) * (tE(iCol,iVer) - 4000._r8) losscf2(iCol,iVer) = losscf2d(iCol,iVer) * (1._r8 / 3000._r8 - 1._r8 / tE(iCol,iVer)) losscf3(iCol,iVer) = -22713._r8 * (1._r8 / tN(iCol,iVer) - 1._r8 / tE(iCol,iVer)) lossceO1D(iCol,iVer) = lossc14p(iCol,iVer) * EXP(losscf2(iCol,iVer)) * & (1._r8 - EXP(losscf3(iCol,iVer))) / tENDiff(iCol,iVer) lossceN2Rot(iCol,iVer) = lossc15p(iCol,iVer) / sqrtTE(iCol,iVer) lossceO2Rot(iCol,iVer) = lossc16p(iCol,iVer) / sqrtTE(iCol,iVer) losscei(iCol,iVer) = lossc3p(iCol,iVer) / tE(iCol,iVer)**1.5_r8 lossg3(iCol,iVer) = lossceN2(iCol,iVer) + lossceO2(iCol,iVer) + lossceO1(iCol,iVer) + lossceN2v(iCol,iVer) & + lossceO2v(iCol,iVer) + lossceOf(iCol,iVer) + lossceO1D(iCol,iVer) & + lossceN2Rot(iCol,iVer) + lossceO2Rot(iCol,iVer) tE(iCol,iVer) = (wrk2(iCol,iVer)/ztodt * tE0(iCol,iVer) + (lossg3(iCol,iVer) * tN(iCol,iVer) & + losscei(iCol,iVer) * ti(iCol,iVer) + sourceg4(iCol,iVer)) & /(SIN(dipMag(iCol,iVer)))**2) / & (wrk2(iCol,iVer)/ztodt + (lossg3(iCol,iVer) + losscei(iCol,iVer))/(SIN(dipMag(iCol,iVer)))**2._r8) if (iVer == 1) WRITE(iulog,*) 'wrk2(iCol,1)', iCol, wrk2(iCol,1) if (iVer == 1) WRITE(iulog,*) 'ndensE(iCol,1)', iCol, ndensE(iCol,1) if (iVer == 1) WRITE(iulog,*) 'tE0(iCol,1)', iCol, tE0(iCol,1) if (iVer == 1) WRITE(iulog,*) 'lossg3(iCol,1)', iCol, lossg3(iCol,1) if (iVer == 1) WRITE(iulog,*) 'tN(iCol,1)', iCol, tN(iCol,1) if (iVer == 1) WRITE(iulog,*) 'losscei(iCol,1)', iCol, losscei(iCol,1) if (iVer == 1) WRITE(iulog,*) 'tI(iCol,1)', iCol, tI(iCol,1) if (iVer == 1) WRITE(iulog,*) 'sourceg4(iCol,1)', iCol, sourceg4(iCol,1) if (iVer == 1) WRITE(iulog,*) 'dipMag(iCol,1)', iCol, dipMag(iCol,1) !------------------------------------------------------------------------------------------------------------------- ! Limit maximum value of electron temperature to be a locally set maximum (temax) and minimum value of electron ! temperature to be at least the neutral temperature for a given column and level !------------------------------------------------------------------------------------------------------------------- tE(iCol,iVer) = min(temax,tE(iCol,iVer)) tE(iCol,iVer) = max(tN(iCol,iVer),tE(iCol,iVer)) enddo !iVer loop endif ! Column not converged enddo ! End of column loop WRITE(iulog,*) 'tE(1,1:ionBot) 5', tE(1,1:ionBot) WRITE(iulog,*) 'wrk2(1,1:ionBot) ', wrk2(1,1:ionBot) WRITE(iulog,*) 'ndensE(1,1:ionBot) ', ndensE(1,1:ionBot) WRITE(iulog,*) 'tE0(1,1:ionBot) ', tE0(1,1:ionBot) WRITE(iulog,*) 'lossg3(1,1:ionBot) ', lossg3(1,1:ionBot) WRITE(iulog,*) 'tN(1,1:ionBot) ', tN(1,1:ionBot) WRITE(iulog,*) 'losscei(1,1:ionBot) ', losscei(1,1:ionBot) WRITE(iulog,*) 'tI(1,1:ionBot) ', tI(1,1:ionBot) WRITE(iulog,*) 'sourceg4(1,1:ionBot) ', sourceg4(1,1:ionBot) ! fvitt -- make sqrtTE consistent with new tE sqrtTE(1:ncol,1:ionBot) = SQRT(tE(1:ncol,1:ionBot)) !----------------------------------------------------- ! Calculate thermal conductivity of electron gas !----------------------------------------------------- do iCol = 1, ncol if (.not. colConv(iCol)) then WRITE(iulog,*) 'Kec1 ', Kec1 WRITE(iulog,*) 'Kec2 ', Kec2 do iVer = 1, ionBot f1Ted1 = 2.82E-17_r8 * sqrtTE(iCol,iVer) - 3.41E-21_r8 * tE(iCol,iVer)**1.5_r8 f1Ted2 = 2.2E-16_r8 + 7.92E-18_r8 * sqrtTE(iCol,iVer) f1Ted3 = 1.1E-16_r8 * (1._r8 + 5.7E-4_r8 * tE(iCol,iVer)) f1Te = ndensN2(iCol,iVer) / ndensE(iCol,iVer) * f1Ted1 + ndensO2(iCol,iVer) / & ndensE(iCol,iVer) * f1Ted2 + ndensO1(iCol,iVer) / ndensE(iCol,iVer) * f1Ted3 !----------------------------------------------------------------------------- ! Calculate electron conductivity using parameters set in module and f1(Te) !----------------------------------------------------------------------------- Ke(iVer) = Kec1 * tE(iCol,iVer)**2.5_r8 / (1._r8 + Kec2 * tE(iCol,iVer)**2._r8 * f1Te) enddo !iVer loop !---------------------------------------------------------------------- ! Get electron conductivity at interface levels to be used later !---------------------------------------------------------------------- do iVer = 2,ionBot Kei(iVer) = SQRT(Ke(iVer-1)*Ke(iVer)) enddo Kei(1) = 1.5_r8*Kei(2)-.5_r8*Kei(3) Kei(ionBotP) = 1.5_r8*Kei(ionBot)-.5_r8*Kei(ionBot-1) WRITE(iulog,*) 'Kei(1:ionBot) ', Kei(1:ionBot) WRITE(iulog,*) 'delZi(iCol,1:ionBot) ', delZi(iCol,1:ionBot) WRITE(iulog,*) 'delZ(iCol,1:ionBot) ', delZ(iCol,1:ionBot) !------------------------------------------------------------------------------------------------------ ! Derive subdiagonal, superdiagonal, and diagonal as input to solver for electron temperature tendency !------------------------------------------------------------------------------------------------------ do iVer = 2, ionBot-1 subDiag(iVer) = -Kei(iVer) / delZi(iCol,iVer) / delZ(iCol,iVer) superDiag(iVer) = -Kei(iVer+1) / delZi(iCol,iVer+1) / delZ(iCol,iVer) diag(iVer) = wrk2(iCol,iVer)/ztodt-subDiag(iVer)-superDiag(iVer) rHS(iVer) = tE(iCol,iVer) * wrk2(iCol,iVer)/ztodt enddo !iVer loop !------------------------------------------------------------------------------------- ! Calculate diagonal, superdiagonal, and right hand side upper boundary values !------------------------------------------------------------------------------------- superDiag(1) = -Kei(2) / delZi(iCol,2) / delZ(iCol,1) diag(1) = wrk2(iCol,1)/ztodt - superDiag(1) rHS(1) = tE(iCol,1) * wrk2(iCol,1) / ztodt - Feub(iCol) / delZ(iCol,1) !--------------------------------------------------------------------------------------------- ! Calculate subdiagonal, diagonal, superdiagonal, and right hand side lower boundary values !--------------------------------------------------------------------------------------------- subDiag(ionBot) = -Kei(ionBot) / delZi(iCol,ionBot) / delZ(iCol,ionBot) superDiag(ionBot) = -Kei(ionBotP) / delZi(iCol,ionBotP) / delZ(iCol,ionBot) diag(ionBot) = wrk2(iCol,ionBot)/ztodt-subDiag(ionBot)-superDiag(ionBot) rHS(ionBot) = tE(iCol,ionBot) * wrk2(iCol,ionBot)/ztodt -superDiag(ionBot)*tN(iCol,ionBotP) WRITE(iulog,*) 'subdiag(1:ionBot) ', subdiag(1:ionBot) WRITE(iulog,*) 'superdiag(1:ionBot) ', superdiag(1:ionBot) WRITE(iulog,*) 'diag(1:ionBot) ', diag(1:ionBot) WRITE(iulog,*) 'rHS(1:ionBot) ', rHS(1:ionBot) !------------------------------------------------- ! Call solver to get electron temperature update !------------------------------------------------- call tridag(subDiag,diag,superDiag,rHS,tETemp,ionBot) WRITE(iulog,*) 'tETemp(1:ionBot) ', tETemp(1:ionBot) tE(iCol,1:ionBot) = tETemp(1:ionBot) do iVer = 1,ionBot tE(iCol,iVer) = min(temax,tE(iCol,iVer)) tE(iCol,iVer) = max(tN(iCol,iVer),tE(iCol,iVer)) enddo !--------------------------------------------------------------------------------------------------------- ! Calculate ion temperature from electron temperature, ion-neutral and electron-ion loss terms, neutral ! temperature, mass density and joule heating. Set minimum value to neutral temperature and maximum ! value to electron temperature for each column and vertical level !--------------------------------------------------------------------------------------------------------- !T do iVer = 1,ionBot !T ti(iCol,iVer) = (losscei(iCol,iVer) * tE(iCol,iVer) + losscin(iCol,iVer) * tN(iCol,iVer) + & !T rho(iCol,iVer) * qjoule(iCol,iVer))/(losscei(iCol,iVer) + losscin(iCol,iVer)) !T ti(iCol,iVer) = max(tN(iCol,iVer),ti(iCol,iVer)) !T ti(iCol,iVer) = min(tE(iCol,iVer),ti(iCol,iVer)) !T enddo WRITE(iulog,*) 'tE(iCol,1:ionBot) ',tE(iCol,1:ionBot) WRITE(iulog,*) 'pMid(iCol,1:ionBot) ',pMid(iCol,1:ionBot) WRITE(iulog,*) 'tEPrevI(iCol,1:ionBot) ',tE(iCol,1:ionBot), tEPrevI(iCol,1:ionBot) WRITE(iulog,*) 'tE(iCol,1:ionBot)/tEPrevI(iCol,1:ionBot) ',tE(iCol,1:ionBot)/tEPrevI(iCol,1:ionBot) !-------------------------------------------------------------------------------------------------------- ! Check for convergence which is a change of electron temperature ratio to previous loop for all levels ! and columns of less than 0.05K. Had to modify this to do convergence check on each column since ! checking all columns in a chunk gives different answers depending on number of tasks and tasks per node. !-------------------------------------------------------------------------------------------------------- if (ALL(ABS(tE(iCol,1:ionBot) / tEPrevI(iCol,1:ionBot) - 1._r8) < 0.05_r8)) then colConv(iCol) = .true. endif endif ! Column not converged enddo ! iCol loop !-------------------------------------------------------------- ! Check to see if all columns have converged and set flag !-------------------------------------------------------------- if (ALL(colConv(1:ncol))) converged = .true. enddo ! End of iteration loop WRITE(iulog,*) 'MIN/MAX tE(1,1:ionBot)) ',MINVAL(tE(1,1:ionBot)), MAXVAL(tE(1,1:ionBot)) if (MAXVAL(tE(1,1:ionBot)) > 5000.) stop write(iulog,*)'lchnk, Number of tempei iterations, converged flag ', iter, converged !-------------------------------------------------------------------------------------------------------- ! Calculate electron-neutral heating and electron-ion Coulomb heating. Then update dry static energy. !-------------------------------------------------------------------------------------------------------- !T do iVer = 1, ionBot !T do iCol = 1, ncol !T sqrtTE(iCol,iVer) = SQRT(tE(iCol,iVer)) !T lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer) !T lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iCol,iVer)) * sqrtTE(iCol,iVer) !T lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iCol,iVer) !T enddo !T enddo !T qen(1:ncol,1:ionBot) = (lossceN2(1:ncol,1:ionBot)+lossceO2(1:ncol,1:ionBot)+lossceO1(1:ncol,1:ionBot)) * & !T (tE(1:ncol,1:ionBot)-tN(1:ncol,1:ionBot)) / rho(1:ncol,1:ionBot) !T qei(1:ncol,1:ionBot) = losscei(1:ncol,1:ionBot) * (tE(1:ncol,1:ionBot)-ti(1:ncol,1:ionBot)) / rho(1:ncol,1:ionBot) !T ptend%s(1:ncol,1:ionBot) = ptend%s(1:ncol,1:ionBot) + (qei(1:ncol,1:ionBot)+qen(1:ncol,1:ionBot))/6.24E15_r8 !T !T call outfld ('TE' , tE , pcols, lchnk) !T call outfld ('TI' , ti , pcols, lchnk) !T call outfld ('LOSS_g3' , lossg3 , pcols, lchnk) !T call outfld ('LOSS_EI' , losscei , pcols, lchnk) !T call outfld ('LOSS_IN' , losscin , pcols, lchnk) !T call outfld ('QIN' , qei , pcols, lchnk) return end subroutine ionos_tempei !=============================================================================== !=============================================================================== subroutine ionos_optransvert(ncol, istate, oPBot) !----------------------------------------------------------------------- ! Routine to compute O+ ion transport in the vertical direction through ! ambipolar diffusion. !----------------------------------------------------------------------- use read_ncfile, only : tE,tI,omega,tN,uN,vN,ed1_geo,ed2_geo,bMag,bEast,bNorth,ndensOp,geoMagLat,geoMagLon use geogrid, only : pMid,pInt use params, only : iulog,rair,gravity,avogad,pi,ztodt,mbar,kboltz use namelist, only : inLat,inLon implicit none !------------------------------Arguments-------------------------------- type(ionos_state), intent(inout) :: istate ! ionosphere state structure integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: oPBot ! bottom of O+ transport calculations !---------------------------Local storage------------------------------- integer :: oPBotP ! bottom of O+ transport calculations plus one more level integer :: iVer ! Counter for vertical loops integer :: iCol ! Counter for column loops integer :: nstep ! current timestep number integer :: indxTe ! pbuf index for electron temperature integer :: indxTi ! pbuf index for ion temperature integer :: indxO1 ! pbuf index atomic O integer :: indxOp,sindxOp,pindxOp ! pbuf or state%q index for O+ mixing ratio integer :: indxE,sindxE,pindxE ! pbuf or state%q index for E mixing ratio integer :: indxO2p ! cnst index for O2+ integer :: indxNOp ! cnst index for NO+ integer :: indxNp ! cnst index for N+ integer :: indxN2p ! cnst index for N2+ real(r8) :: cBF ! Burnside factor real(r8) :: rMassOp ! Molecular mass of O+ (kg/kmole) real(r8) :: tR ! Reduced temperature (K) real(r8) :: dAPrime ! Factor needed for dA (cm*cm/s/K real(r8) :: delTpDelZ ! Partial derivative of Tp with z real(r8) :: delTpDelZ1 ! Top level of partial derivative of Tp with z real(r8) :: pScaleHeight ! Pressure scale height real(r8) :: bMagT ! Magnetic field magnitude in Teslas (tesla) real(r8) :: c31,c41,c43r,c5,c6 ! Coefficients for upper boundary O+ calculations real(r8) :: aOpFluxCoeff ! Coefficient for O+ flux calculation real(r8) :: oPFlux ! O+ flux (?) real(r8) :: ndensOpTop ! Top (zeroth) level O+ flux (?) real(r8), dimension(pcols) :: epot_geo ! electric potential on geographic grid real(r8), dimension(pcols) :: mlt ! mag.local time of WACCM geographic grid points real(r8), dimension(pcols) :: decMag ! magnetic declination angle (radians) real(r8), dimension(pcols) :: sinDec ! sine of the magnetic declination angle real(r8), dimension(pcols) :: cosDec ! cosine of the magnetic declination angle real(r8), dimension(pcols) :: eHPB ! Horizontal electric field perpendicular to magnetic field [V/m] real(r8), dimension(:,:), pointer :: mmrPE ! Pointer to access electrons in pbuf real(r8), dimension(:,:), pointer :: mmrPOp ! Pointer to access O+ in pbuf real(r8), dimension(pcols,pver) :: mmrE ! E electron mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: tP ! Plasma temperature (K) real(r8), dimension(pcols,pverp) :: tNInt ! Interface Neutral Temperture (K) real(r8), dimension(pcols,pverp) :: tIInt ! Interface Ion Temperture (K) real(r8), dimension(pcols,pverp) :: tEInt ! Interface Electron Temperture (K) real(r8), dimension(pcols,pverp) :: tPInt ! Interface Plasma Temperture (K) real(r8), dimension(pcols,pverp) :: rairvi ! Constituent dependent gas constant on interface levels real(r8), dimension(pcols,pver) :: rairv ! Constituent dependent gas constant on interface levels real(r8), dimension(pcols,pver) :: wN ! Neutral vertical wind velocity (m/s) real(r8), dimension(pcols,pver) :: ndensN2 ! N2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3) real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3) real(r8), dimension(pcols,pver) :: dipMag ! dip angle for each column (radians) real(r8), dimension(pcols,pver) :: dipMagD ! dip angle for each column (degrees) real(r8), dimension(ncol) :: zenAngD ! zenith angle (degrees) real(r8), dimension(pcols) :: Feub ! electron heat flux at upper boundary real(r8), dimension(pcols,pver) :: delZ ! Delta z: midpoints (cm) real(r8), dimension(pcols,pverp) :: delZi ! Delta z: interfaces (cm) real(r8), dimension(pcols,pver) :: dA ! Ambipolar diffusion coefficient for c1 and c2 coefficients (cm*cm/s) real(r8), dimension(pcols,pver) :: c1 ! c1 coefficient for matrix diagonals: midpoints real(r8), dimension(pcols,pver) :: c2 ! c1 coefficient for matrix diagonals: midpoints real(r8), dimension(pcols,pverp) :: c1i ! c2 coefficient for matrix diagonals: interfaces real(r8), dimension(pcols,pverp) :: c2i ! c2 coefficient for matrix diagonals: interfaces real(r8), dimension(oPBot) :: subdiag ! subdiagonal values for Te tendency solving real(r8), dimension(oPBot) :: superdiag ! superdiagonal values for Te tendency solving real(r8), dimension(oPBot) :: diag ! diagonal values for Te tendency solving real(r8), dimension(oPBot) :: rHS ! RHS of electron temperature update real(r8), dimension(oPBot) :: newDensOp ! Op result from tridag solver real(r8), dimension(oPBot) :: newMMROp ! O+ mass mixing ratio converted from O+ number density from solver real(r8), dimension(oPBot) :: newMMRE ! Electron mass mixing ratio from new O+ summed with O2+ and NO+ real(r8), dimension(pcols) :: vdrift !----------------------------------------------------------------------------------------------------------------- write(iulog,*)'lchnk, Inside optransvert ' !--------------------------------------------------------------------------------------------------------- ! Initialize arrays to zero !--------------------------------------------------------------------------------------------------------- epot_geo(:) = 0._r8 mlt(:) = 0._r8 decMag(:) = 0._r8 sinDec(:) = 0._r8 cosDec(:) = 0._r8 ! geoMagLat(:) = 0._r8 eHPB(:) = 0._r8 tP(:,:) = 0._r8 tNInt(:,:) = 0._r8 tIInt(:,:) = 0._r8 tEInt(:,:) = 0._r8 rairvi(:,:) = 0._r8 rairv(:,:) = 0._r8 wN(:,:) = 0._r8 ndensN2(:,:) = 0._r8 ndensO2(:,:) = 0._r8 ndensO1(:,:) = 0._r8 dipMag(:,:) = 0._r8 dipMagD(:,:) = 0._r8 zenAngD(:) = 0._r8 Feub(:) = 0._r8 delZ(:,:) = 0._r8 delZi(:,:) = 0._r8 dA(:,:) = 0._r8 c1(:,:) = 0._r8 c2(:,:) = 0._r8 c1i(:,:) = 0._r8 c2i(:,:) = 0._r8 subDiag(:) = 0._r8 superDiag(:) = 0._r8 diag(:) = 0._r8 rHS(:) = 0._r8 newDensOp(:) = 0._r8 newMMROp(:) = 0._r8 newMMRE(:) = 0._r8 !------------------------------------------------------------------------------- ! Set Burnside Factor and bottom level plus one to calculate O+ transport !------------------------------------------------------------------------------- cBF = 0.8_r8 oPBotP = oPBot + 1 !---------------------------------------------------------------- ! Get variables needed from the ionosphere state structure !---------------------------------------------------------------- dipMag(1:ncol,1:pver) = istate%dipMag(1:ncol,1:pver) ndensO1(1:ncol,1:pver) = istate%ndensO1(1:ncol,1:pver) ndensN2(1:ncol,1:pver) = istate%ndensN2(1:ncol,1:pver) ndensO2(1:ncol,1:pver) = istate%ndensO2(1:ncol,1:pver) tNInt(1:ncol,1:pverp) = istate%tNInt(1:ncol,1:pverp) rairvi(1:ncol,1:pverp) = istate%rairvi(1:ncol,1:pverp) rairv(1:ncol,1:pver) = istate%rairv(1:ncol,1:pver) zenAngD(1:ncol) = istate%zenAngD(1:ncol) !---------------------------------------------------------------- ! Get O+ mass !---------------------------------------------------------------- rmassOP = 16._r8 !---------------------------------------------------------------- ! Get variables needed from the physics state structure !---------------------------------------------------------------- !------------------------------------------------------------------------------------------------------------------- ! Get electron temperature from physics buffer and if this is the first time calculated then initialize to neutral ! temperature. nsrest=0 means this is an initial run and nstep=0 means this is the first time step. !------------------------------------------------------------------------------------------------------------------- nstep = 1 !-------------------------------------------------- ! Calculate cosine/sine of declination angle !-------------------------------------------------- decMag(1:ncol) = -atan2( beast,bnorth ) cosDec(1:ncol) = cos( decMag(1:ncol) ) sinDec(1:ncol) = sin( decMag(1:ncol) ) !------------------------------------------------------------------------------------------------------------ ! Calculate the magnetic local time of WACCM geographic grid points needed for mag to geo conversion !------------------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------------------------ ! Convert e_field module zonal and meridional electric field components to geographic coords using exbdrift module ! map_mag2geo method !------------------------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------- ! Calculate ion temperature on interface levels !------------------------------------------------------------------------------------- do iVer = 2, oPBotP do iCol = 1, ncol tIInt(iCol,iVer) = 0.5_r8 * tI(iCol,iVer) + 0.5_r8 * tI(iCol,iVer-1) enddo enddo !-------------------------- ! Get top interface level !-------------------------- do iCol = 1, ncol tIInt(iCol,1) = 1.5_r8 * tIInt(iCol,2) - 0.5_r8 * tIInt(iCol,3) enddo !-------------------------- ! Get bottom interface level !-------------------------- do iCol = 1, ncol tIInt(iCol,pverp) = 1.5_r8 * tIInt(iCol,pver) - 0.5_r8 * tIInt(iCol,pver-1) enddo !------------------------------------------------------------------------------------- ! Calculate electron temperature on interface levels !------------------------------------------------------------------------------------- do iVer = 2, oPBotP do iCol = 1, ncol tEInt(iCol,iVer) = 0.5_r8 * tE(iCol,iVer) + 0.5_r8 * tE(iCol,iVer-1) enddo enddo !-------------------------- ! Get top interface level !-------------------------- do iCol = 1, ncol tEInt(iCol,1) = 1.5_r8 * tEInt(iCol,2) - 0.5_r8 * tEInt(iCol,3) enddo !-------------------------- ! Get bottom interface level !-------------------------- ! do iCol = 1, ncol ! ! tEInt(iCol,pverp) = 1.5_r8 * tEInt(iCol,pver) - 0.5_r8 * tEInt(iCol,pver-1) ! ! enddo !-------------------------------------------------------------------------------------------------------- ! Get plasma temperature needed for delTp/delz !-------------------------------------------------------------------------------------------------------- do iVer = 1, oPBot do iCol = 1, ncol tP(iCol,iVer) = ( tI(iCol,iVer) + tE(iCol,iVer) ) / 2._r8 enddo enddo !-------------------------- ! Get bottom mid-point level !-------------------------- do iCol = 1, ncol tP(iCol,oPBotP) = 1.5_r8 * tP(iCol,oPBot) - 0.5_r8 * tP(iCol,oPBot-1) enddo do iVer = 1, oPBotP do iCol = 1, ncol tPInt(iCol,iVer) = 0.5_r8 * (tEInt(iCol,iVer) + tIInt(iCol,iVer)) enddo enddo !------------------------------------------------------------------------------------------ ! Calculate thickness terms for vertical derivative !------------------------------------------------------------------------------------------ do iVer = 1, oPBot do iCol = 1, ncol delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rairv(iCol,iVer) * tN(iCol,iVer) / pMid(iCol,iVer) / gravity enddo enddo do iVer = 2, oPBotP ! Assuming oPBotP < pverp do iCol = 1, ncol delZi(iCol,iVer) = (pMid(iCol,iVer) - pMid(iCol,iVer-1)) * rairvi(iCol,iVer) * & tNInt(iCol,iVer) / pInt(iCol,iVer) / gravity enddo enddo do iCol = 1, ncol delZi(iCol,1) = 1.5_r8 * delZi(iCol,2) - .5_r8 * delZi(iCol,3) !---------------------------------------------------------- ! Convert delZ variables from meters to centimeters !---------------------------------------------------------- delZi(iCol,1:oPBotP) = delZi(iCol,1:oPBotP) * 100._r8 delZ(iCol,1:oPBot) = delZ(iCol,1:oPBot) * 100._r8 enddo write(iulog,*)'Tn', tN(1,1:oPBot) !----------------------------------------------------------------------------- ! Compute the horizontal electric field perpendicular to the magnetic field !----------------------------------------------------------------------------- eHPB(1:ncol) = ed1_geo(1:ncol) * cosDec(1:ncol) - ed2_geo(1:ncol) * sinDec(1:ncol) !------------------------------------------------------------------ ! Calculate pressure scale height and vertical velocity !------------------------------------------------------------------ do iVer = 2,oPBot do iCol = 1,ncol pScaleHeight = .5_r8*(rairv(iCol,iVer)*tN(iCol,iVer)+rairv(iCol,iVer)*tN(iCol,iVer-1))/gravity wN(iCol,iVer) = -omega(iCol,iVer) * pScaleHeight ! write(iulog,*)'iVer,iCol,wN(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight,rairv,tN(iCol,iVer),gravity ', iVer,iCol,wN(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight,rairv(iCol,iVer),tN(iCol,iVer),gravity enddo enddo !-------------------------- ! Get top midpoint level !-------------------------- do iCol = 1, ncol wN(iCol,1) = 1.5_r8 * wN(iCol,2) - 0.5_r8 * wN(iCol,3) enddo write(iulog,*)'vN(1,1:oPBot) ', vN(1,1:oPBot) ! write(iulog,*)'lchnk, min/max ndensO1(1:ncol,1:oPBot) ', lchnk, MINVAL(ndensO1(1:ncol,1:oPBot)), MAXVAL(ndensO1(1:ncol,1:oPBot)) ! write(iulog,*)'lchnk, min/max ndensN2(1:ncol,1:oPBot) ', lchnk, MINVAL(ndensN2(1:ncol,1:oPBot)), MAXVAL(ndensN2(1:ncol,1:oPBot)) ! write(iulog,*)'lchnk, min/max ndensO2(1:ncol,1:oPBot) ', lchnk, MINVAL(ndensO2(1:ncol,1:oPBot)), MAXVAL(ndensO2(1:ncol,1:oPBot)) ! write(iulog,*)'lchnk, Before vertical loop 2 optransvert ' !-------------------------------------------------------------------------------------------------------- ! Loop over vertical levels then columns and call tridiagonal solver for each column to get O+ !-------------------------------------------------------------------------------------------------------- do iVer = 1, oPBot do iCol = 1, ncol !------------------------------------------------------------------------------- ! Calculate Da ambipolar diffusion coefficient for c1 and c2 coefficients (cgs) !------------------------------------------------------------------------------- tR = ( tI(iCol,iVer) + tN(iCol,iVer) ) *.5_r8 dAPrime = 1.42E+17_r8 / (SQRT(tR) * (1._r8 - 0.064_r8 * LOG10(tR))**2 * ndensO1(iCol,iVer) * cBF + 18.6_r8 * & ndensN2(iCol,iVer) + 18.1_r8 * ndensO2(iCol,iVer)) ! if (iVer == 1) write(iulog,*)'lchnk, iCol, tR ', lchnk, iCol, tR ! if (iVer == 1) write(iulog,*)'lchnk, iCol, ndensO1,(iCol,iVer), ndensN2(iCol,iVer), ndensO2(iCol,iVer) ', lchnk, iCol, ndensO1(iCol,iVer), ndensN2(iCol,iVer), ndensO2(iCol,iVer) ! if (iVer == 1) write(iulog,*)'lchnk, iCol, dAPrime ', lchnk, iCol, dAPrime ! if (iVer == 1) write(iulog,*)'lchnk, iCol, tP ', lchnk, iCol, tP(iCol,iVer) dA(iCol,iVer) = dAPrime * 2._r8 * tP(iCol,iVer) !------------------------------------------------------------------------------- ! Calculate derivative for c1 and c2 coefficients (cgs) !------------------------------------------------------------------------------- delTpDelZ = ( tP(iCol,iVer) - tP(iCol,iVer+1) ) / delZ(iCol,iVer) !------------------------------------------------------------- ! Need to save top value for upper boundary calculations !------------------------------------------------------------- if (iVer == 1) delTPDelZ1 = delTpDelZ !----------------------------------------------------------------------------- ! Convert magnetic field magnitude units from nanoTesla to Tesla !----------------------------------------------------------------------------- ! bMagT = bmag*1.e-9_r8 bMagT = bmag*1.e-4_r8 !--------------------------------------------------------------------------------- ! Compute C1 coefficient for matrix formation and RHS ! eHPB(V/m) bMagT(nT=V*s/m2) so eHPB/bMagT=(V/m/V/s*m2=m/s=.01*cm/s) ! Factors .01 and 100. is to convert from SI to cgs units (meters to centimeters) !--------------------------------------------------------------------------------- vdrift(iCol) = eHPB(iCol) / bMagT * cos(dipMag(iCol,iVer)) * 100._r8 ! c1(iCol,iVer) = -vdrift(iCol) + & ! ( ( vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * sin(decMag(iCol)) ) * & ! cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ) * & ! sin(dipMag(iCol,iVer)) * 100._r8 + & ! dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + & ! 0.5_r8 * rmassOp / avogad * gravity / kboltz / tP(iCol,iVer) * .01_r8 ) * & ! sin(dipMag(iCol,iVer))**2 c1(iCol,iVer) = 2000.+( ( vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * sin(decMag(iCol)) ) * & cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ) * & sin(dipMag(iCol,iVer)) * 100._r8 + & dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + & 0.5_r8 * rmassOp / avogad * gravity / kboltz / tP(iCol,iVer) * .01_r8 ) * & sin(dipMag(iCol,iVer))**2 ! c1(iCol,iVer) = dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + & ! 0.5_r8 * rmassOp / avogad * gravity / kboltz / tP(iCol,iVer) * .01_r8 ) * & ! sin(dipMag(iCol,iVer))**2 ! c1(iCol,iVer) = -eHPB(iCol) / bMagT * cos(dipMag(iCol,iVer)) * 100._r8 + & ! ( ( vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * sin(decMag(iCol)) ) * & ! cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ) * & ! sin(dipMag(iCol,iVer)) * 100._r8 + & ! dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + & ! 0.5_r8 * rmassOp / avogad * gravity / kboltz / tP(iCol,iVer) * .01_r8 ) * & ! sin(dipMag(iCol,iVer))**2 if(c1(iCol,iVer) <= 0._r8) write(iulog,*)'neg c1, iVer, delTpDelZ, tP, c1',iVer, delTpDelZ, & tP(iCol,iVer),c1(iCol,iVer) !-------------------------------------------------------------- ! Compute C2 coefficient for matrix formation ! dA(cm2/s) !-------------------------------------------------------------- c2(iCol,iVer) = dA(iCol,iVer) * sin(dipMag(iCol,iVer))**2 enddo enddo write(iulog,*)'lchnk, After vertical loop optransvert ' write(iulog,*)'eHPB(1), bMagT, cos(decMag(1) ', eHPB(1), bMagT , cos(dipMag(1,1)) write(iulog,*)'vdrift',vdrift(1) write(iulog,*)'dipMag(1,1:oPBot) ', dipMag(1,1:oPBot) write(iulog,*)'vN(1,1:oPBot) ', vN(1,1:oPBot) write(iulog,*)'uN(1,1:oPBot) ', uN(1,1:oPBot) write(iulog,*)'wN(1,1:oPBot) ', wN(1,1:oPBot) write(iulog,*)'dA(1,1:oPBot) ', dA(1,1:oPBot) write(iulog,*)'delTpDelZ ', delTpDelZ write(iulog,*)'tP(1,1:oPBot) ', tP(1,1:oPBot) write(iulog,*)'rairv(1,1:oPBot) ', rairv(1,1:oPBot) write(iulog,*)'rmassOp, avogad, gravity, kboltz ', rmassOp, avogad, gravity, kboltz !---------------------------------------------------- ! Get c1 and c2 coefficients on interface levels !---------------------------------------------------- do iVer = 2, oPBot ! Assuming oPBot < pverp do iCol = 1, ncol c1i(iCol,iVer) = 0.5_r8 * c1(iCol,iVer) + 0.5_r8 * c1(iCol,iVer-1) c2i(iCol,iVer) = 0.5_r8 * c2(iCol,iVer) + 0.5_r8 * c2(iCol,iVer-1) enddo enddo !------------------------------------------------------- ! Get top interface level for c1i and c2i coefficients !------------------------------------------------------- do iCol = 1, ncol ! write(iulog,*)'lchnk,iCol,ztodt,c1(iCol,1),c1(iCol,2) ',lchnk,iCol,ztodt,c1(iCol,1),c1(iCol,2) ! write(iulog,*)'lchnk, iCol, min/max dA(iCol,1:oPBot) ', lchnk, iCol, MINVAL(dA(iCol,1:oPBot)), MAXVAL(dA(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max c1(iCol,1:oPBot) ', lchnk, iCol, MINVAL(c1(iCol,1:oPBot)), MAXVAL(c1(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max c2(iCol,1:oPBot) ', lchnk, iCol, MINVAL(c2(iCol,1:oPBot)), MAXVAL(c2(iCol,1:oPBot)) c1i(iCol,1) = 1.5_r8 * c1i(iCol,2) - 0.5_r8 * c1i(iCol,3) c2i(iCol,1) = 1.5_r8 * c2i(iCol,2) - 0.5_r8 * c2i(iCol,3) ! write(iulog,*)'lchnk,iCol,ztodt,c1i(iCol,1),c1i(iCol,2) ',lchnk,iCol,ztodt,c1i(iCol,1),c1i(iCol,2) enddo !--------------------------------------------------------- ! Get bottom interface level for c1i and c2i coefficients !--------------------------------------------------------- do iCol = 1, ncol c1i(iCol,oPBotP) = 1.5_r8 * c1i(iCol,oPBot) - 0.5_r8 * c1i(iCol,oPBot-1) c2i(iCol,oPBotP) = 1.5_r8 * c2i(iCol,OPBot) - 0.5_r8 * c2i(iCol,oPBot-1) enddo ! write(iulog,*)'lchnk, Before column loop optransvert ' ! write(iulog,*)'lchnk,iCol,ztodt,c1i(nCol,oPBot),c1i(nCol,oPBot-1),c1i(nCol,oPBotP)',lchnk,nCol,ztodt,c1i(nCol,oPBot),c1i(nCol,oPBot-1),c1i(nCol,oPBotP) !-------------------------------------------------------------------------------------------------------- ! Loop over columns then vertical levels and call tridiagonal solver for each column to get O+ !-------------------------------------------------------------------------------------------------------- do iCol = 1, ncol ! write(iulog,*)'lchnk,iCol,ztodt,c1i(iCol,1),c1i(iCol,2) ',lchnk,iCol,ztodt,c1i(iCol,1),c1i(iCol,2) ! write(iulog,*)'lchnk, iCol, ztodt, min/max c1i(iCol,1:oPBot) ', lchnk, iCol, ztodt, MINVAL(c1i(iCol,1:oPBot)), MAXVAL(c1i(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max delZ(iCol,1:oPBot) ', lchnk, iCol, MINVAL(delZ(iCol,1:oPBot)), MAXVAL(delZ(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max c2i(iCol,1:oPBot) ', lchnk, iCol, MINVAL(c2i(iCol,1:oPBot)), MAXVAL(c2i(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max delZi(iCol,1:oPBot) ', lchnk, iCol, MINVAL(delZi(iCol,1:oPBot)), MAXVAL(delZi(iCol,1:oPBot)) do iVer = 2, oPBot-1 !----------------------------------------------------------------------------------------------------- ! Calculate sub-diagonal, diagonal, superdiagonal, and right hand side !----------------------------------------------------------------------------------------------------- subDiag(iVer) = -( 0.25_r8 * c1i(iCol,iVer) / delZ(iCol,iVer) + 0.5_r8 *c2i(iCol,iVer) / & delZi(iCol,iVer) / delZ(iCol,iVer) ) diag(iVer) = 1._r8 / ztodt - 0.25_r8 * ( c1i(iCol,iVer) - c1i(iCol,iVer+1) ) / delZ(iCol,iVer) + 0.5_r8 * & ( c2i(iCol,iVer) / delZi(iCol,iVer) /delZ(iCol,iVer) + c2i(iCol,iVer+1) / delZ(iCol,iVer) / & delZi(iCol,iVer+1) ) superDiag(iVer) = 0.25_r8 * c1i(iCol,iVer+1) / delZ(iCol,iVer) - 0.5_r8 * c2i(iCol,iVer+1) / & delZ(iCol,iVer) / delZi(iCol,iVer+1) rHS(iVer) = -subDiag(iVer) * ndensOp(iCol,iVer-1) + & ( 1._r8 / ztodt + 0.25_r8 * ( c1i(iCol,iVer) - c1i(iCol,iVer+1) ) / delZ(iCol,iVer) - & 0.5_r8 * ( c2i(iCol,iVer) / delZi(iCol,iVer) /delZ(iCol,iVer) + & c2i(iCol,iVer+1) / delZ(iCol,iVer) / delZi(iCol,iVer+1) ) ) * ndensOp(iCol,iVer) - & superDiag(iVer) * ndensOp(iCol,iVer+1) enddo !----------------------------------------------------------------- ! Get upper boundary conditions for diagonal and RHS in cgs units !----------------------------------------------------------------- ! Can we make the flux a constant array, so to save some computing each time we go through the iCol loop? ?-HL !----------------------------------------------------------------- ! Calculate O+ flux !----------------------------------------------------------------- ! write(iulog,*)'geoMagLat(iCol)',geoMagLat(iCol) ! if (ABS(geoMagLat(iCol)) > pi/12._r8) aOPFluxCoeff = 1._r8 ! if (ABS(geoMagLat(iCol)) <= pi/12._r8) aOPFluxCoeff = 0.5_r8 * ( 1._r8 + sin((ABS(geoMagLat(iCol)) - pi / 24._r8) * 12._r8) ) write(iulog,*)'geoMagLat',geoMagLat if (ABS(geoMagLat) > pi/12._r8) aOPFluxCoeff = 1._r8 if (ABS(geoMagLat) <= pi/12._r8) aOPFluxCoeff = 0.5_r8 * ( 1._r8 + sin((ABS(geoMagLat) - pi / 24._r8) * 12._r8) ) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,lat,state%lat(iCol),lon,state%lon(iCol)',lchnk,iCol,get_lat_p(lchnk,iCol),state%lat(iCol),get_lon_p(lchnk,iCol),state%lon(iCol) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,geoMagLat(iCol),aOPFluxCoeff',lchnk,iCol,geoMagLat(iCol),aOPFluxCoeff if (zenAngD(iCol) <= 80._r8) oPFlux = dayOPFlux * aOPFluxCoeff if (zenAngD(iCol) > 80._r8 .and. zenAngD(iCol) < 100._r8) & oPFlux = ( 0.5_r8 * (dayOPFlux + nightOPFlux) + 0.5_r8 *(dayOPFlux - nightOPFlux) * cos(pi * (zenAngD(iCol) - 80._r8) / 20._r8) ) * aOPFluxCoeff if (zenAngD(iCol) >= 100._r8) oPFlux = nightOPFlux * aOPFluxCoeff write(iulog,*)'zenAngD,oPFlux,dayOPflux,aOPFluxCoeff,nightOPFlux',zenAngD(iCol),oPFlux,dayOPFlux,aOPFluxCoeff,nightOPFlux !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,zenAngD(iCol),OPFlux,dayOPFlux,nightOpFlux',lchnk,iCol,zenAngD(iCol),OPFlux,dayOPFlux,nightOpFlux !----------------------------------------------------------------- ! Calculate coefficients needed for diagonal and RHS top values !----------------------------------------------------------------- c31 = -dA(iCol,1) * sin(dipMag(iCol,1))**2 c43r = ( 2._r8 * delTpDelZ1 + rmassOp / avogad * gravity / kboltz * .01_r8 )/2._r8/tP(iCol,1) c41 = c31 * c43r !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c31,c41',lchnk,iCol,c31,c41 !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c31,dA(iCol,1),dipmag(iCol,1)',lchnk,iCol,c31,dA(iCol,1),dipMag(iCol,1) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c41,tP(iCol,1),delTPDelZ1,rmassOp,avogad,gravit,kboltz',lchnk,iCol,c41,tP(iCol,1),delTPDelZ1,rmassOp,avogad,gravit,kboltz if (dipMag(iCol,1) == 0._r8) c5 = 0._r8 if (dipMag(iCol,1) /= 0._r8) c5 = 2._r8 * delZ(iCol,1) / c31 * oPFlux if (dipMag(iCol,1) == 0._r8) c6 = 1._r8 if (dipMag(iCol,1) /= 0._r8) c6 = - 2._r8 * delZ(iCol,1) * c43r !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c5,c6',lchnk,iCol,c5,c6 !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c5,c6,dA(iCol,1),dipmag(iCol,1),c31,c41,delZ(iCol,1)',lchnk,iCol,c5,c6,dA(iCol,1),dipmag(iCol,1),c31,c41,delZ(iCol,1) !----------------------------------------------------------------- ! Calculate diagonals and RHS top level values !----------------------------------------------------------------- subDiag(1) = -( 0.25_r8 * c1i(iCol,1) / delZ(iCol,1) + 0.5_r8 *c2i(iCol,1) / & delZi(iCol,1) / delZ(iCol,1) ) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,subDiag(1),c1i(iCol,1),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1),c2i(iCol,2),delZi(iCol,2) ',lchnk,iCol,subDiag(1),c1i(iCol,1),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1),c2i(iCol,2),delZi(iCol,2) diag(1) = 1._r8 / ztodt - 0.25_r8 * ( c1i(iCol,1) - c1i(iCol,2) ) / delZ(iCol,1) + 0.5_r8 * & ( c2i(iCol,1) / delZi(iCol,1) /delZ(iCol,1) + c2i(iCol,2) / delZ(iCol,1) / & delZi(iCol,2) ) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,diag(1),ztodt,c1i(iCol,1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1),c2i(iCol,2),delZi(iCol,2) ',lchnk,iCol,diag(1),ztodt,c1i(iCol,1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1),c2i(iCol,2),delZi(iCol,2) write(iulog,*)'ztodt,c1i(iCol,1:oPBot) ',ztodt,c1i(iCol,1:oPBot) !Swrite(iulog,*)'delZ(iCol,1:oPBot) ', delZ(iCol,1:oPBot) !Swrite(iulog,*)'c2i(iCol,1:oPBot) ',c2i(iCol,1:oPBot) !Swrite(iulog,*)'delZi(iCol,1:oPBot) ',delZi(iCol,1:oPBot) superDiag(1) = 0.25_r8 * c1i(iCol,2) / delZ(iCol,1) - 0.5_r8 * c2i(iCol,2) / & delZ(iCol,1) / delZi(iCol,2) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,superDiag(1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,2),delZi(iCol,2) ',lchnk,iCol,superDiag(1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,2),delZi(iCol,2) ndensOpTop = c5 + c6 * ndensOP(iCol,1) + ndensOP(iCol,2) rHS(1) = - subDiag(1) * ndensOpTop + & ( 1._r8 / ztodt + 0.25_r8 * ( c1i(iCol,1) - c1i(iCol,2) ) / delZ(iCol,1) - & 0.5_r8 * ( c2i(iCol,1) / delZi(iCol,1) /delZ(iCol,1) + & c2i(iCol,2) / delZ(iCol,1) / delZi(iCol,2) ) ) * ndensOp(iCol,1) - & superDiag(1) * ndensOp(iCol,2) write(iulog,*)'1st rHS',rHS(1),ndensOpTop,c5,c6,c43r,ndensOP(iCol,1),ndensOP(iCol,2),delZ(iCol,1) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,rHS(1),subDiag(1),ndensOpTop,ztodt,c1i(iCol,1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1)',lchnk,iCol,rHS(1),subDiag(1),ndensOpTop,ztodt,c1i(iCol,1),c1i(iCol,2),delZ(iCol,1),c2i(iCol,1),delZi(iCol,1) !if(ABS(state%lat(iCol)) < 0.2 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,c2i(iCol,2),delZi(iCol,2),ndensOp(iCol,1),superDiag(1),ndensOp(iCol,2) ',lchnk,iCol,c2i(iCol,2),delZi(iCol,2),ndensOp(iCol,1),superDiag(1),ndensOp(iCol,2) !-------------------------------------------------------------- ! Set upper boundary conditions for diagonal and RHS !-------------------------------------------------------------- diag(1) = subDiag(1) * c6 + diag(1) superDiag(1) = superDiag(1) + subDiag(1) rHS(1) = rHS(1) - subDiag(1) * c5 write(iulog,*)'2nd rHS',rHS(1) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,diag(1),subdiag(1),superDiag(1),rHS(1) ',lchnk,iCol,diag(1),subdiag(1),superDiag(1),rHS(1) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,diag(1),subdiag(1),c6 ',lchnk,iCol,diag(1),subdiag(1),c6,superDiag(1) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,rHS(1),subdiag(1),c5 ',lchnk,iCol,rHS(1),subdiag(1),c5 !---------------------------------------------- ! Calculate diagonal and RHS bottom values !---------------------------------------------- subDiag(oPBot) = -( 0.25_r8 * c1i(iCol,oPBot) / delZ(iCol,oPBot) + 0.5_r8 *c2i(iCol,oPBot) / & delZi(iCol,oPBot) / delZ(iCol,oPBot) ) diag(oPBot) = 1._r8 / ztodt - 0.25_r8 * ( c1i(iCol,oPBot) - c1i(iCol,oPBotP) ) / delZ(iCol,oPBot) + 0.5_r8 * & ( c2i(iCol,oPBot) / delZi(iCol,oPBot) /delZ(iCol,oPBot) + c2i(iCol,oPBotP) / delZ(iCol,oPBot) / & delZi(iCol,oPBotP) ) superDiag(oPBot) = 0.25_r8 * c1i(iCol,oPBotP) / delZ(iCol,oPBot) - 0.5_r8 * c2i(iCol,oPBotP) / & delZ(iCol,oPBot) / delZi(iCol,oPBotP) rHS(oPBot) = -subDiag(oPBot) * ndensOp(iCol,oPBot-1) + & ( 1._r8 / ztodt + 0.25_r8 * ( c1i(iCol,oPBot) - c1i(iCol,oPBotP) ) / delZ(iCol,oPBot) - & 0.5_r8 * ( c2i(iCol,oPBot) / delZi(iCol,oPBot) /delZ(iCol,oPBot) + & c2i(iCol,oPBotP) / delZ(iCol,oPBot) / delZi(iCol,oPBotP) ) ) * ndensOp(iCol,oPBot) - & superDiag(oPBot) * ndensOp(iCol,oPBotP) !-------------------------------------------------------------- ! Set lower boundary condition for RHS !-------------------------------------------------------------- rHS(oPBot) = rHS(oPBot) - superDiag(oPBot) * ndensOp(iCol,oPBotP) ! write(iulog,*)'lchnk, iCol, min/max subdiag(1:oPBot) ', lchnk, iCol, MINVAL(subdiag(1:oPBot)), MAXVAL(subdiag(1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max diag(1:oPBot) ', lchnk, iCol, MINVAL(diag(1:oPBot)), MAXVAL(diag(1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max superdiag(1:oPBot) ', lchnk, iCol, MINVAL(superdiag(1:oPBot)), MAXVAL(superdiag(1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max rHS(1:oPBot) ', lchnk, iCol, MINVAL(rHS(1:oPBot)), MAXVAL(rHS(1:oPBot)) write(iulog,*)' subdiag(1:oPBot) ', subdiag(1:oPBot) write(iulog,*)' diag(1:oPBot) ', diag(1:oPBot) write(iulog,*)' superdiag(1:oPBot) ', superdiag(1:oPBot) write(iulog,*)' rHS(1:oPBot) ', rHS(1:oPBot) !------------------------------------------------- ! Call solver to get O+ number density update !------------------------------------------------- call tridag(subDiag,diag,superDiag,rHS,newDensOp,oPBot) write(iulog,*)'iCol, lat, lon, zenAngD(iCol), ndensOp(iCol,1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), ndensOp(iCol,1:oPBot) write(iulog,*)'iCol, lat, lon, zenAngD(iCol), newDensOp(1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), newDensOp(1:oPBot) !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,newDensOp(1) ',lchnk,iCol,newDensOp(1) !if(newDensOp(1) < 0.) write(iulog,*)'lchnk,iCol,newDensOp(1) ',lchnk,iCol,newDensOp(1) ! if (newDensOp(1) < 0._r8) write(iulog,*)'lchnk, iCol, !newDensOp(1:OPBot) = ABS(newDensOp(1:oPBot)) ! if (newDensOp(1) < 0._r8) write(iulog,*)'lchnk, iCol, lat, lon, newDensOp(1:oPBot) ', lchnk, iCol, get_lat_p(lchnk,iCol),get_lon_p(lchnk,iCol), newDensOp(1:oPBot) ! if (iCol < 3) write(iulog,*)'iCol, lat, lon, zenAngD(iCol),ndensOp(iCol,1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), ndensOp(iCol,1:oPBot) ! if (iCol < 3) write(iulog,*)'iCol, lat, lon, zenAngD(iCol), newDensOp(1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), newDensOp(1:oPBot) ! write(iulog,*)'lchnk, iCol, min/max ndensOp ', lchnk, iCol, MINVAL(ndensOp(iCol,1:oPBot)), MAXVAL(ndensOp(iCol,1:oPBot)) ! write(iulog,*)'lchnk, iCol, min/max newDensOp ', lchnk, iCol, MINVAL(newDensOp(1:oPBot)), MAXVAL(newDensOp(1:oPBot)) ! write(iulog,*)'lchnk, min/max mbarv(iCol,1:oPBot,lchnk) ', lchnk, MINVAL(mbarv(iCol,1:oPBot,lchnk)), MAXVAL(mbarv(iCol,1:oPBot,lchnk)) ! write(iulog,*)'lchnk, min/max pMid(iCol,1:oPBot) ', lchnk, MINVAL(pMid(iCol,1:oPBot)), MAXVAL(pMid(iCol,1:oPBot)) ! write(iulog,*)'lchnk, min/max tN(iCol,1:oPBot) ', lchnk, MINVAL(tN(iCol,1:oPBot)), MAXVAL(tN(iCol,1:oPBot)) ! write(iulog,*)'lchnk, min/max mbarv(iCol,1:oPBot,lchnk) ', lchnk, MINVAL(mbarv(iCol,1:oPBot,lchnk)), MAXVAL(mbarv(iCol,1:oPBot,lchnk)) !------------------------------------------------------------------------------------------------------------------------------- ! Convert new O+ number density to mixing ratio and output to pbuf or state%q depending on whether it is short lived or not !------------------------------------------------------------------------------------------------------------------------------- ! write(iulog,*)'lchnk, min/max before update newMMROp(1:oPBot) ', lchnk, MINVAL(newMMROp(1:oPBot)), MAXVAL(newMMROp(1:oPBot)) newMMROp(1:oPBot) = newDensOp(1:oPBot) / mbar * rmassOp / pMid(iCol,1:oPBot) * (kboltz * tN(iCol,1:oPBot)) / 1.E-06_r8 ! write(iulog,*)'lchnk, min/max after update newMMROp(1:oPBot) ', lchnk, MINVAL(newMMROp(1:oPBot)), MAXVAL(newMMROp(1:oPBot)) if (ANY(newMMROp(1:oPBot) < 0._r8)) then write(iulog,*)'lat,lon,zenAngD(iCol) negativemmr newMMROp(1:oPBot) ', inLat,inLon, zenAngD(iCol), newMMROp(1:oPBot) stop endif do iVer = 1,3 if (newMMROp(iVer) < 1.E-09_r8) newMMROp(iVer) = 1.E-09_r8 enddo !if(ABS(state%lat(iCol)) < 0.5 .and. state%lon(iCol) == 0.) write(iulog,*)'lchnk,iCol,newMMROp(1),mmrOp(iCol,1) ',lchnk,iCol,newMMROp(1),mmrOp(iCol,1) write(iulog,*)'iCol, lat, lon, zenAngD(iCol), after update newMMROp(1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), newMMROp(1:oPBot) ! write(iulog,*)'iCol, lat, lon, zenAngD(iCol), after update rmassOp, mbar ', iCol, inLat, inLon, zenAngD(iCol), rmassOp, mbar ! write(iulog,*)'iCol, lat, lon, zenAngD(iCol), after update pMid(iCol,1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), pMid(iCol,1:oPBot) ! write(iulog,*)'iCol, lat, lon, zenAngD(iCol), after update kboltz, tN(iCol,1:oPBot) ', iCol, inLat, inLon, zenAngD(iCol), kboltz, tN(iCol,1:oPBot) ! write(iulog,*)'lchnk, min/max newMMROp(1:oPBot) ', lchnk, MINVAL(newMMROp(1:oPBot)), MAXVAL(newMMROp(1:oPBot)) !---------------------------------------------------------------- ! Get O2+ and NO+ indices !---------------------------------------------------------------- !S call cnst_get_ind( 'O2p', indxO2p ) !S call cnst_get_ind( 'NOp', indxNOp ) !S call cnst_get_ind( 'Np' , indxNp ) !S call cnst_get_ind( 'N2p', indxN2p ) enddo ! iCol loop write(iulog,*) 'O+ output ndens ', newDensOp(1:oPBot) write(iulog,*) 'O+ output mmr', newMMROp(1:oPBot) !------------------------------------------------------------------------------------------------------------------------------- ! Call routine to calculate new E mass mixing ratio from ion sum and output to pbuf or state%q depending on whether it is short ! lived or not !------------------------------------------------------------------------------------------------------------------------------- !S call charge_fix( ncol, state%q(:,:,:), pbuf, lchnk ) ! call outfld ('TE' , tE , pcols, lchnk) ! call outfld ('TI' , ti , pcols, lchnk) ! call outfld ('LOSS_g3' , lossg3 , pcols, lchnk) ! call outfld ('LOSS_EI' , losscei , pcols, lchnk) ! call outfld ('LOSS_IN' , losscin , pcols, lchnk) ! call outfld ('QIN' , qei , pcols, lchnk) !S call outfld ('OpI' , mmrPOp , pcols, lchnk) !S call outfld ('eI' , mmrPE , pcols, lchnk) return end subroutine ionos_optransvert !=============================================================================== !----------------------------------------------------------------------- ! Simple tridiagonal solver routine !----------------------------------------------------------------------- SUBROUTINE tridag(a,b,c,r,u,n) !S use abortutils, only: endrun INTEGER n REAL(r8) :: a(n),b(n),c(n),r(n),u(n) INTEGER j REAL(r8) :: bet,gam(n) !S if(b(1).eq.0._r8) call endrun('ionosphere: bt(1)=0 in tridag') if(b(1).eq.0._r8) write(iulog,*) 'ionosphere: bt(1)=0 in tridag' bet=b(1) u(1)=r(1)/bet do j=2,n gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) !S if(bet.eq.0._r8) call endrun('ionosphere: bet=0 in tridag') if(bet.eq.0._r8) write(iulog,*) 'ionosphere: bet=0 in tridag' u(j)=(r(j)-a(j)*u(j-1))/bet end do do j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) end do return END SUBROUTINE tridag end module ionosphere