!--------------------------------------------------------------------------------- ! 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) !------------------------------- ! Get electron temperature !------------------------------- ! call ionos_tempei(state, ptend, lchnk, ncol, ztodt, pbuf, 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) !------------------------------------------------------------------------------ ! 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,zenAngD,tN,mmrO1,mmrO2,omega,pMid,tE,tI 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 !---------------------------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) :: 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) :: mmrNO ! NO 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 !-------------------------------------------------------------------------------- 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 ndensE(:,:) = 0._r8 ndensO2p(:,:) = 0._r8 ndensNOp(:,:) = 0._r8 sourceR(:,:) = 0._r8 sourceEff(:,:) = 0._r8 sourceg4(:,:) = 0._r8 !------------------------------------------------------------------------------------------------------ ! 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 !-------------------------------------------------------------- rads2Degs = 180._r8/pi degs2Rads = pi/180._r8 !------------------------------------------------------------------------------------- ! 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(iCol,pver) - 0.5_r8 * omegai(iCol,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 !----------------------------------------------------------------------------------------------------- ! 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 ndensO2p(iCol,iVer) = mmrO2p(iCol,iVer) * mbar / rmassO2p * & pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 ndensNOp(iCol,iVer) = mmrNOp(iCol,iVer) * mbar / rmassNOp * & 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 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%ndensE(1:ncol,1:pver) = ndensE(1:ncol,1:pver) istate%ndensO2p(1:ncol,1:pver) = ndensO2p(1:ncol,1:pver) istate%ndensNOp(1:ncol,1:pver) = ndensNOp(1:ncol,1:pver) istate%ndensN2(1:ncol,1:pver) = ndensN2(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) !------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------------------- ! Calculate electron heating rate which is source in electron/ion temperature derivation !------------------------------------------------------------------------------------------------- !S!------------------------------------------------------------------------------- !S! Calculate g4 source term for electron temperature update !S!------------------------------------------------------------------------------- return end subroutine ionos_datinit ! !=============================================================================== !=============================================================================== 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 use geogrid, only : pMid,pInt use params, only : iulog,rair,gravity,avogad,pi,ztodt,mbar,kboltz 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) :: geoMagLat ! Geomagnetic latitude (radians) 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 = ( tPInt(iCol,iVer) - tPInt(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 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) ) !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,*)'2rd 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) !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,*)'lchnk, iCol, min/max ndensOp(iCol,1:oPBot) ', lchnk, iCol, ndensOp(iCol,1:oPBot) ! if (iCol < 3) write(iulog,*)'lchnk, iCol, min/max newDensOp(1:oPBot) ', lchnk, 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)) 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,*)'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