!--------------------------------------------------------------------------------- ! 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 !S use ppgrid, only : pcols, pver, pverp ! Dimensions and chunk bounds !S use cam_history, only : outfld ! Routine to output fields to history files !S use physics_types, only : physics_state, physics_ptend !Structures containing physics state and tendency variables !S use physics_buffer, only : pbuf_add_field, pbuf_get_index,dtype_r8,physics_buffer_desc, pbuf_get_field, pbuf_set_field, pbuf_get_chunk !! Needed to get variables from physics buffer !S use cam_logfile, only : iulog ! Output unit for run.out file !S use mo_jeuv, only : nIonRates ! Number of ionization rates in mo_photo !S use shr_const_mod, only : kboltz => shr_const_boltz, pi => shr_const_pi ! Boltzmann constant and pi !S use geogrid, only : pver,pverp use params, only : iulog implicit none save private ! Make default type private to the module ! ! PUBLIC: interfaces ! !S public ionos_init ! Initialization !S public ionos_register ! Registration of ionosphere variables in pbuf physics buffer 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 !S real(r8), dimension(pcols,pver,nIonRates) :: ionPRates ! ionization rates temporary array (s-1 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(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 end type ionos_state contains !============================================================================== !S subroutine ionos_init() !S !S!----------------------------------------------------------------------- !S! Time independent initialization for ionosphere simulation. !S!----------------------------------------------------------------------- !S use cam_history, only : addfld, add_default, phys_decomp ! Routines and variables for adding fields to history output !S use dycore, only : get_resolution !S use interpolate_data, only : lininterp !S !S implicit none !S !S!------------------------------------------------------------------------------- !S! Add history variables for ionosphere !S!------------------------------------------------------------------------------- !S call addfld ('IONOS_NDENSN2', '1/m3',pver, 'I','N2 Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSO2', '1/m3',pver, 'I','O2 Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSO1', '1/m3',pver, 'I','O1 Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSNO', '1/m3',pver, 'I','NO Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSN1', '1/m3',pver, 'I','NO Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSE' , '1/m3',pver, 'I','E Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSOP' ,'1/m3',pver, 'I','OP Number Density-Ionos', phys_decomp) !S call addfld ('IONOS_NDENSO2P','1/m3',pver, 'I','O2P Number Density-Ionos',phys_decomp) !S call addfld ('IONOS_NDENSNOP','1/m3',pver, 'I','NOP Number Density-Ionos',phys_decomp) !S !S call addfld ('TE' , 'K',pver, 'I','Electron Temperature', phys_decomp) !S call addfld ('TI' , 'K',pver, 'I','Ion Temperature', phys_decomp) !S call addfld ('QIN' , 'J/kg/s',pver, 'I','JOULE+IN Heating', phys_decomp) !S call addfld ('LOSS_g3' , ' ',pver , 'I','Loss Term g3', phys_decomp) !S call addfld ('LOSS_EI' , ' ',pver , 'I','Loss Term EI', phys_decomp) !S call addfld ('LOSS_IN' , ' ',pver , 'I','Loss Term IN', phys_decomp) !S call addfld ('IONI_RATE' , ' ',pver , 'I','Total ionization', phys_decomp) !S call addfld ('IONI_Eff' , ' ',pver , 'I','ionization efficiency', phys_decomp) !S call addfld ('SOURCE_g4' , ' ',pver , 'I','SOURCE g4', phys_decomp) !S call addfld ('AUR_IRATESUM' , ' ',pver , 'I','Auroral ionization', phys_decomp) !S !S call addfld ('OpI' , ' ',pver , 'I','O+ Ionosphere', phys_decomp) !S call addfld ('eI' , ' ',pver , 'I','e Ionosphere', phys_decomp) !S !S!------------------------------------------------------------------------------- !S! Set default values for ionosphere history variables !S!------------------------------------------------------------------------------- !S call add_default ('IONOS_NDENSN2' , 1, ' ') !S call add_default ('IONOS_NDENSO2' , 1, ' ') !S call add_default ('IONOS_NDENSO1' , 1, ' ') !S call add_default ('IONOS_NDENSNO' , 1, ' ') !S call add_default ('IONOS_NDENSN1' , 1, ' ') !S call add_default ('IONOS_NDENSE' , 1, ' ') !S call add_default ('IONOS_NDENSOP' , 1, ' ') !S call add_default ('IONOS_NDENSO2P', 1, ' ') !S call add_default ('IONOS_NDENSNOP', 1, ' ') !S call add_default ('TE' , 1, ' ') !S call add_default ('TI' , 1, ' ') !S call add_default ('QIN' , 1, ' ') !S call add_default ('LOSS_g3' , 1, ' ') !S call add_default ('LOSS_EI' , 1, ' ') !S call add_default ('LOSS_IN' , 1, ' ') !S call add_default ('IONI_RATE' , 1, ' ') !S call add_default ('IONI_Eff' , 1, ' ') !S call add_default ('SOURCE_g4' , 1, ' ') !S call add_default ('AUR_IRATESUM' , 1, ' ') !S !S call add_default ('OpI' , 1, ' ') !S call add_default ('eI' , 1, ' ') !S !S return !S !S end subroutine ionos_init !============================================================================== !S subroutine ionos_register !S !S!----------------------------------------------------------------------- !S! Register ionosphere variables with physics buffer: !S! !S! Ion production rates pcols,pver,nIonRates, !S! so firstdim = 1 middledim = pver lastdim = nIonRates. !S! !S! pcols dimension and lchnk assumed here !S! !S!----------------------------------------------------------------------- !S !S implicit none !S !S integer :: idx !S !S call pbuf_add_field('ElecTemp','global',dtype_r8,(/pcols,pver/),idx) ! Electron temperature (global so can write to history files) !S !S call pbuf_add_field('IonTemp', 'global',dtype_r8,(/pcols,pver/),idx) ! Ion temperature (global so can write to history files) !S !S end subroutine ionos_register ! !============================================================================== !S! subroutine ionos_intr(state, ptend, pbuf, ztodt) subroutine ionos_intr(ncol,oPBot,ionBot) !----------------------------------------------------------------------- ! Interface for improved ionosphere simulation !----------------------------------------------------------------------- ! !------------------------------Arguments-------------------------------- !S use time_manager, only : get_nstep,get_step_size implicit none !S type(physics_state), intent(inout) :: state ! physics state structure !S type(physics_ptend), intent(inout) :: ptend ! parameterization tendency structure !S type(physics_buffer_desc),pointer :: pbuf(:) ! physics buffer type(ionos_state) :: istate ! ionosphere state structure !S real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) 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 !S integer :: ncol ! number of atmospheric columns !S integer :: ionBot ! bottom of ionosphere calculations !S integer :: oPBot ! bottom of O+ transport calculations !---------------------------------------------------------------- ! Get current chunk number and number of columns in this chunk !---------------------------------------------------------------- !S lchnk = state%lchnk !S ncol = state%ncol !------------------------------------------------------------ ! 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 !=============================================================================== !S subroutine ionos_datinit(state, lchnk, ncol, pbuf, istate, ionBot, oPBot) 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 !------------------------------------------------------------------------------ !S use cam_history, only : addfld, add_default, phys_decomp ! Routines/variables needed for outputting fields !S use addfld_mod, only : addfld ! Routines/variables needed for outputting fields !S use dycore, only : get_resolution !S use interpolate_data, only : lininterp !S use constituents, only : cnst_get_ind, cnst_mw ! Routines to get molecular weights for constituents !S use hycoef, only : hypm ! Model vertical pressure grid in Pascals !S use mo_apex, only : bnorth, beast, bdown ! Magnetic field components !S use time_manager, only : get_curr_calday ! Routine to get current calendar day !S use mo_solar_parms, only : get_solar_parms ! Routine to get solar parameters, i.e. f107 !S use cam_control_mod, only : nsrest ! Variable to determine if this is an initial run or a restart/branch !S use time_manager, only : get_nstep ! Routine to get current time step !S use physconst, only : rairv, mbarv ! Constituent dependent rair and mbar !S use short_lived_species, only : slvd_index,slvd_pbf_ndx => pbf_idx ! Routines to access short lived species in storage and pbuf use read_ncfile, only: bEast,bNorth,bDown,zenAngD,tN,mmrO1,mmrO2,omega,pMid,tE,tI use params, only: rair,mbar,pi,kboltz implicit none !S type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer !S type(physics_state), intent(in) :: state ! physics state structure type(ionos_state), intent(inout) :: istate ! ionosphere state structure !S integer, intent(in) :: lchnk ! Chunk number integer, intent(in) :: ncol ! Number of columns in current chunk !S integer, intent(out) :: ionBot ! bottom of ionosphere calculations !S integer, intent(out) :: oPBot ! bottom of ionosphere calculations !---------------------------Local storage------------------------------- !S integer :: indxUe ! pbuf index for eastward ExB drift !S integer :: indxVe ! pbuf index for northward ExB drift !S integer :: indxWe ! pbuf index for upward ExB drift !S integer :: indxIR ! pbuf index for ionization rates !S integer :: indxAIPRS ! pbuf index for aurora ion production rate sum !S integer :: indxTe ! pbuf index for electron temperature !S integer :: indxTi ! pbuf index for ion temperature !S integer :: indxN2 ! cnst index for N2 !S integer :: indxO1 ! cnst index for O !S integer :: indxO2 ! cnst index for O2 !S integer :: indxNO ! cnst index for NO !S integer :: indxH ! cnst index for H !S integer :: indxN1 ! cnst index for N !S integer :: indxE ! cnst index for electrons !S integer :: indxOp ! cnst index for Op !S integer :: indxO2p ! cnst index for O2p !S integer :: indxNOp ! cnst index for NOp integer :: iVer ! Counter for vertical loops integer :: iCol ! Counter for column loops !S integer :: iIonR ! Counter for ionization rates loops !S integer :: indxSP ! pbuf index for Pedersen Conductivity !S integer :: indxSH ! pbuf index for Hall Conductivity integer :: nstep ! current timestep number !S real(r8), dimension(:), pointer :: ue ! pointer to eastward ExB drift in pbuf (from module iondrag) !S real(r8), dimension(:), pointer :: ve ! pointer to northward ExB drift in pbuf !S real(r8), dimension(:), pointer :: we ! pointer to upward ExB drift in pbuf !S real(r8), dimension(pcols,pver) :: tE ! Pointer to electron temperature in pbuf (K) !S real(r8), dimension(pcols,pver) :: ti ! Pointer to ion temperature in pbuf (K) !S real(r8), dimension(:,:,:), pointer :: ionRates ! Pointer to ionization rates for O+,O2+,N+,N2+,NO+ in pbuf (s-1 from modules mo_jeuv and mo_jshort) !S real(r8), dimension(:,:), pointer :: sigma_ped ! Pointer to Pedersen Conductivity in pbuf (siemens/m) from module iondrag !S real(r8), dimension(:,:), pointer :: sigma_hall ! Pointer to Hall Conductivity in pbuf (siemens/m) !S real(r8), dimension(:,:), pointer :: mmrP ! Pointer to access short lived species in pbuf real(r8), dimension(ncol) :: geoLatR ! Latitude (radians) Make ncol because zenith aurora are ncol real(r8), dimension(ncol) :: geoLonR ! Longitude (radians) !S real(r8), dimension(pcols,pver) :: pMid ! Midpoint pressure (Pa) !S real(r8), dimension(pcols,pver) :: tN ! Neutral temperature (K) 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) !S real(r8), dimension(ncol) :: 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) :: 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 !S real(r8), dimension(pcols,pver) :: omega ! vertical velocity at midpoint levels (Pa/s) 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 !S real(r8), dimension(pcols,pver) :: mmrO2 ! O2 mass mixing ratio kg/kg !S real(r8), dimension(pcols,pver) :: mmrO1 ! O 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 !S real(r8), dimension(pcols,pver) :: mmrOp ! O plus ion 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) !S real(r8), dimension(pcols,pver,nIonRates) :: ionPRates ! ionization rates temporary array (s-1 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 !S zenAngD(:) = 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 !S ionPRates(:,:,:) = 0._r8 sumIonPRates(:,:) = 0._r8 ndensN2(:,:) = 0._r8 ndensO2(:,:) = 0._r8 ndensO1(:,:) = 0._r8 ndensE(:,:) = 0._r8 !S ndensOp(:,:) = 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. !------------------------------------------------------------------------------------------------------ !S do iVer = 1, pver !S !S if (hypm(iVer) < 50._r8) ionBot = iVer !S !S enddo !------------------------------------------------------------------------------------------------------------------------ ! 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 !---------------------------------------------------------------- ! Get latitude and longitude of each column in this chunk !---------------------------------------------------------------- !S geoLatR(1:ncol) = state%lat(1:ncol) !S geoLonR(1:ncol) = state%lon(1:ncol) !------------------------------------------------------------------------------------------------------- ! Need to get midpoint and interface pressure and neutral temperature from state structure (pcols,pver) !------------------------------------------------------------------------------------------------------- !S pMid(1:ncol,1:pver) = state%pmid(1:ncol,1:pver) !S tN(1:ncol,1:pver) = state%t(1:ncol,1:pver) !------------------------------------------------------------------------------------- ! 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) !-------------------------------------------------------------- ! Get zenith angle !-------------------------------------------------------------- !S calDay = get_curr_calday() !S call zenith(calDay,geoLatR,geoLonR,cosZenAngR,ncol) !S do iCol = 1, ncol !S !S zenAngD(iCol) = ACOS(cosZenAngR(iCol)) * rads2Degs !S zenAngD(iCol) = ACOS(cosZenAngR(iCol)) * rads2Degs !S !S enddo cosZenAngR(:) = COS(zenAngD * degs2Rads) istate%cosZenAngR(1:ncol) = cosZenAngR(1:ncol) !S istate%zenAngD(1:ncol) = zenAngD(1:ncol) istate%zenAngD(1:ncol) = zenAngD !-------------------------------------------------------------- ! Get F10.7 solar flux !-------------------------------------------------------------- !S call get_solar_parms( f107_s = f107 ) !--------------------------------------------------------------------------------------- ! Expand magnetic field components in vertical to make 3D, pcols,pver,begchunk:endchunk ! These are used in calculation of magnetic dip angle and magnetic declination angle so ! store in local ionosphere module structure. !--------------------------------------------------------------------------------------- !S do iVer = 1, pver !S !S do iCol = 1, ncol !S !S bNorth3d(iCol,iVer) = bnorth(iCol,lchnk) !S bEast3d(iCol,iVer) = beast(iCol,lchnk) !S bDown3d(iCol,iVer) = bdown(iCol,lchnk) !S !S enddo !S !S enddo !S istate%bNorth3d(1:ncol,1:pver) = bNorth3d(1:ncol,1:pver) !S istate%bEast3d(1:ncol,1:pver) = bEast3d(1:ncol,1:pver) !S istate%bDown3d(1:ncol,1:pver) = bDown3d(1:ncol,1:pver) 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) !--------------------------------------------------------------------------------- !S indxUe = pbuf_get_index( 'UE' ) !S indxVe = pbuf_get_index( 'VE' ) !S indxWe = pbuf_get_index( 'WE' ) !S call pbuf_get_field(pbuf, indxUe, ue) !S call pbuf_get_field(pbuf, indxVe, ve) !S call pbuf_get_field(pbuf, indxWe, we) !------------------------------------------------------------------------------- ! Form 2D drifts needed for this module. Vertical level goes to ionBotP since ! needed below to get vertical drift on interface levels !------------------------------------------------------------------------------- !S do iVer = 1, pver !S !S do iCol = 1, ncol !S !S ue2d(iCol,iVer) = ue(iCol) !S ve2d(iCol,iVer) = ve(iCol) !S we2d(iCol,iVer) = we(iCol) !S !S enddo !S !S enddo !S !S istate%ue2d(1:ncol,1:pver) = ue2d(1:ncol,1:pver) !S istate%ve2d(1:ncol,1:pver) = ve2d(1:ncol,1:pver) !S istate%we2d(1:ncol,1:pver) = we2d(1:ncol,1:pver) !------------------------------------------------------------------------------- ! 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 !-------------------------------------------------------------------------------------- !S omega(1:ncol,1:pver) = state%omega(1:ncol,1:pver) !S !S do iVer = 2, pver !S do iCol = 1, ncol !S omegai(iCol,iVer) = 0.5_r8 * (omega(iCol,iVer-1) + omega(iCol,iVer)) !S enddo !S enddo !S do iCol = 1, ncol !S omegai(iCol,1) = 1.5_r8 * omegai(iCol,2) - 0.5_r8 * omegai(iCol,3) !S enddo !S do iCol = 1, ncol !S omegai(iCol,pverp) = 1.5_r8 * omegai(iCol,pver) - 0.5_r8 * omegai(iCol,pver-1) !S enddo 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 !------------------------------------------------------------------------ !S do iVer = 2, pver !S do iCol = 1, ncol !S rairvi(iCol,iVer) = 0.5_r8 * rairv(iCol,iVer-1,lchnk) + 0.5_r8 * rairv(iCol,iVer,lchnk) !S enddo !S enddo !S !S do iCol = 1, ncol !S rairvi(iCol,1) = 1.5_r8 * rairvi(iCol,2) - 0.5_r8 * rairvi(iCol,3) !S enddo !S do iCol = 1, ncol !S rairvi(iCol,pverp) = 1.5_r8 * rairvi(iCol,pver) - 0.5_r8 * rairvi(iCol,pver-1) !S enddo !S istate%rairvi(1:ncol,1:pverp) = rairvi(1:ncol,1:pverp) istate%rairvi(1:ncol,1:pverp) = rair !------------------------------------------------------------------------------- ! Need to get dip angle from magnetic field components !------------------------------------------------------------------------------- do iVer = 1, pver do iCol = 1, ncol !S dipMag(iCol,iVer) = ATAN(bDown3d(iCol,iVer) / SQRT(bNorth3d(iCol,iVer)**2 + bEast3d(iCol,iVer)**2)) 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 !---------------------------------------------------------------------------------------- !S call cnst_get_ind( 'O', indxO1, abort=.false. ) !S if (indxO1 < 0) then !S indxO1 = slvd_index( 'O' ) !S if (indxO1 > 0) then !S xo1_slvd = .true. !S endif !S else !S xo1_slvd = .false. !S endif !S call cnst_get_ind( 'O2', indxO2, abort=.false. ) !S if (indxO2 < 0) then !S indxO2 = slvd_index( 'O2' ) !S if (indxO2 > 0) then !S xo2_slvd = .true. !S endif !S else !S xo2_slvd = .false. !S endif !S call cnst_get_ind( 'NO', indxNO, abort=.false. ) !S if (indxNO < 0) then !S indxNO = slvd_index( 'NO' ) !S if (indxNO > 0) then !S xno_slvd = .true. !S endif !S else !S xno_slvd = .false. !S endif !S call cnst_get_ind( 'H', indxH, abort=.false. ) !S if (indxH < 0) then !S indxH = slvd_index( 'H' ) !S if (indxH > 0) then !S xh_slvd = .true. !S endif !S else !S xh_slvd = .false. !S endif !S call cnst_get_ind( 'N', indxN1, abort=.false. ) !S if (indxN1 < 0) then !S indxN1 = slvd_index( 'N' ) !S if (indxN1 > 0) then !S xn_slvd = .true. !S endif !S else !S xn_slvd = .false. !S endif !S call cnst_get_ind( 'e', indxE, abort=.false. ) !S if (indxE < 0) then !S indxE = slvd_index( 'e' ) !S if (indxE > 0) then !S elec_slvd = .true. !S endif !S else !S elec_slvd = .false. !S endif !S call cnst_get_ind( 'Op', indxOp, abort=.false. ) !S if (indxOp < 0) then !S indxOp = slvd_index( 'Op' ) !S if (indxOp > 0) then !S op_slvd = .true. !S endif !S else !S op_slvd = .false. !S endif !S call cnst_get_ind( 'O2p', indxO2p, abort=.false. ) !S if (indxO2p < 0) then !S indxO2p = slvd_index( 'O2p' ) !S if (indxO2p > 0) then !S o2p_slvd = .true. !S endif !S else !S o2p_slvd = .false. !S endif !S call cnst_get_ind( 'NOp', indxNOp, abort=.false. ) !S if (indxNOp < 0) then !S indxNOp = slvd_index( 'NOp' ) !S if (indxNOp > 0) then !S nop_slvd = .true. !S endif !S else !S nop_slvd = .false. !S endif !---------------------------------------------------------------------------------------- ! Get molecular weights using constituent indices (kg/kmole) except N2 which is hard-wired !---------------------------------------------------------------------------------------- !S rmassO2 = cnst_mw(indxO2) !S rmassO1 = cnst_mw(indxO1) !S rmassNO = cnst_mw(indxNO) !S rmassH = cnst_mw(indxH) !S rmassN1 = cnst_mw(indxN1) !S rmassE = cnst_mw(indxE) !S rmassOP = cnst_mw(indxOp) !S rmassO2P = cnst_mw(indxO2p) !S rmassNOP = cnst_mw(indxNOp) 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 !----------------------------------------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------------------------------------- !S if ( xo2_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxO2/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxO2) !S mmrO2(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrO2(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxO2) !S endif !S if ( xo1_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxO1/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxO1) !S mmrO1(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrO1(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxO1) !S endif !S if ( xno_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxNO/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxNO) !S mmrNO(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrNO(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxNO) !S endif !S if ( xn_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxN1/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxN1) !S mmrN1(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrN1(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxN1) !S endif !S if ( elec_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxE/), kount=(/pcols,pver,1/) ) !S mmrE(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrE(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxE) !S endif !S if ( op_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxOp/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxOp) !S mmrOp(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrOp(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxOp) !S endif !S if ( o2p_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxO2p/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxO2p) !S mmrO2p(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrO2p(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxO2p) !S endif !S if ( nop_slvd ) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrP, start=(/1,1,indxNOp/), kount=(/pcols,pver,1/) ) !S! mmrP => pbuf(slvd_pbf_ndx)%fld_ptr(1,:,:,lchnk,indxNOp) !S mmrNOp(1:ncol,1:pver) = mmrP(1:ncol,1:pver) !S else !S mmrNOp(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxNOp) !S endif !S !S mmrN2(1:ncol,1:pver) = 1._r8 - (mmrO2(1:ncol,1:pver) + mmrO1(1:ncol,1:pver)) !S mmrN2(1:ncol,1:pver) = MAX(1.e-20_r8,mmrN2(1:ncol,1:pver)) !S !-------------------------------------------------------------------------------------------------------------- ! 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 !S ndensO2(iCol,iVer) = mmrO2(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassO2 * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensO1(iCol,iVer) = mmrO1(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassO1 * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensNO(iCol,iVer) = mmrNO(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassNO * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensN1(iCol,iVer) = mmrN1(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassN1 * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensE(iCol,iVer) = mmrE(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassE * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensOp(iCol,iVer) = mmrOp(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassOp * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensO2p(iCol,iVer) = mmrO2p(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassO2p * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensNOp(iCol,iVer) = mmrNOp(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassNOp * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 !S ndensN2(iCol,iVer) = mmrN2(iCol,iVer) * mbarv(iCol,iVer,lchnk) / rmassN2 * & !S pMid(iCol,iVer) / (kboltz * tN(iCol,iVer)) * 1.E-06_r8 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 !S ndensOp(iCol,iVer) = mmrOp(iCol,iVer) * mbar / rmassOp * & !S 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) = mmrN2(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) !S istate%ndensOp(1:ncol,1:pver) = ndensOp(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 !------------------------------------------------------------------------------- !S call outfld ('IONOS_NDENSN2' , ndensN2 , pcols, lchnk) !S call outfld ('IONOS_NDENSO1' , ndensO1 , pcols, lchnk) !S call outfld ('IONOS_NDENSO2' , ndensO2 , pcols, lchnk) !S call outfld ('IONOS_NDENSNO' , ndensNO , pcols, lchnk) !S call outfld ('IONOS_NDENSN1' , ndensN1 , pcols, lchnk) !S call outfld ('IONOS_NDENSE' , ndensE , pcols, lchnk) !S call outfld ('IONOS_NDENSOP' , ndensOp , pcols, lchnk) !S call outfld ('IONOS_NDENSO2P', ndensO2p , pcols, lchnk) !S call outfld ('IONOS_NDENSNOP', ndensNOp , pcols, lchnk) !------------------------------------------------------------------------------------ ! 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 !------------------------------------------------------------------------------------ !S indxIR = pbuf_get_index( 'IonRates' ) !S call pbuf_get_field(pbuf, indxIR, ionRates, start=(/1,1,1/), kount=(/pcols,pver,nIonRates/)) !---------------------------------------------------------------------------------------------- ! 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) !---------------------------------------------------------------------------------------------- !S do iVer = 1, pver !S do iCol = 1, ncol !S do iIonR = 1, nIonRates !S IF (iIonR <= 3) ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensO1(iCol,iVer) !S IF (iIonR == 4) ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensN1(iCol,iVer) !S IF ((iIonR == 5) .OR. (iIonR >= 7 .AND. iIonR <= 9)) & !S ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensO2(iCol,iVer) !S IF (iIonR == 6 .OR. iIonR == 10 .OR. iIonR == 11) & !S ionPRates(iCol,iVer,iIonR) = ionRates(iCol,iVer,iIonR) * ndensN2(iCol,iVer) !S enddo !---------------------------------------------- ! Sum ion production rates all reactions !---------------------------------------------- !S sumIonPRates(iCol,iVer) = SUM(ionPRates(iCol,iVer,1:11)) !S enddo !S enddo !S !S istate%ionPRates(1:ncol,1:pver,1:nIonRates) = ionPRates(1:ncol,1:pver,1:nIonRates) !S istate%sumIonPRates(1:ncol,1:pver) = sumIonPRates(1:ncol,1:pver) !------------------------------------------------------------------------------------------- ! 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 !------------------------------------------------------------------------------------------- !S indxAIPRS = pbuf_get_index( 'AurIPRateSum' ) !S call pbuf_get_field(pbuf, indxAIPRS, aurIPRateSum) !------------------------------------------------------------------------------------- ! 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) !------------------------------------------------------------------------------------- !S do_aurora(:) = abs( geoLatR(:) ) > pi/6._r8 !S if( all( .not. do_aurora(:) ) ) then !S aurIPRateSum(:,:) = 0._r8 !S end if !------------------------------------------------------------------------------------------------- ! Calculate electron heating rate which is source in electron/ion temperature derivation !------------------------------------------------------------------------------------------------- !S do iVer = 1, ionBot !S do iCol = 1, ncol !S sourceR(iCol,iVer) = LOG( ndensE(iCol,iVer) / (ndensO2(iCol,iVer) + ndensN2(iCol,iVer) + & !S 0.1_r8 * ndensO1(iCol,iVer)) ) !S sourceEff(iCol,iVer) = EXP( -(12.75_r8 + 6.941_r8 * sourceR(iCol,iVer) + 1.166_r8 * sourceR(iCol,iVer)**2 + & !S 0.08043_r8 * sourceR(iCol,iVer)**3 + 0.001996_r8 * sourceR(iCol,iVer)**4) ) !S !S!------------------------------------------------------------------------------- !S! Calculate g4 source term for electron temperature update !S!------------------------------------------------------------------------------- !S sourceg4(iCol,iVer) = (sumIonPRates(iCol,iVer) + aurIPRateSum(iCol,iVer)) * sourceEff(iCol,iVer) !S !S enddo !S !S enddo !S !S istate%sourceg4(1:ncol,1:pver) = sourceg4(1:ncol,1:pver) !S call outfld ('IONI_RATE' , sumIonPRates , pcols, lchnk) !S call outfld ('IONI_Eff' , sourceEff , pcols, lchnk) !S call outfld ('SOURCE_g4' , sourceg4 , pcols, lchnk) !S call outfld ('AUR_IRATESUM' , aurIPRateSum , pcols, lchnk) !---------------------------------------------------------------------------------------------- ! Get Pedersen and Hall Conductivities from physics buffer which were calculated in iondrag ! module. Conductivity array dimensions are pcols, pver !------------------------------------------------------------------------------- !S indxSP = pbuf_get_index( 'PedConduct' ) !S indxSH = pbuf_get_index( 'HallConduct' ) !S call pbuf_get_field(pbuf, indxSP, sigma_ped) !S call pbuf_get_field(pbuf, indxSH, sigma_hall) return end subroutine ionos_datinit ! !=============================================================================== !S subroutine ionos_tempei(state, ptend, lchnk, ncol, ztodt, pbuf, istate, ionBot) !S !S!----------------------------------------------------------------------- !S! Routine to compute the electron and ion temperature needed for !S! ambipolar diffusion calculations. !S!----------------------------------------------------------------------- !S !S use cam_history, only : addfld, add_default, phys_decomp !S use phys_grid, only : get_lat_p, get_lon_p, get_rlat_p,get_rlon_p !S use physconst, only : gravit ! Gravity (m/s2) !S use cam_control_mod, only : nsrest ! Variable to determine if this is an initial run or a restart/branch !S use time_manager, only : get_nstep ! Routine to get current time step !S use physconst, only : rairv, mbarv ! Constituent dependent rair and mbar !S !S implicit none !S !S!------------------------------Arguments-------------------------------- !S !S type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer !S type(physics_state), intent(in) :: state ! physics state structure !S type(physics_ptend), intent(inout) :: ptend ! parameterization tendency structure !S type(ionos_state), intent(inout) :: istate ! ionosphere state structure !S !S real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) !S !S !S integer, intent(in) :: lchnk ! chunk identifier !S integer, intent(in) :: ncol ! number of atmospheric columns !S integer, intent(in) :: ionBot ! bottom of ionosphere calculations !S !S!---------------------------Local storage------------------------------- !S integer :: ionBotP ! bottom of ionosphere calculations plus one more level !S !S integer :: i,k ! loop indexes !S integer :: iVer ! Counter for vertical loops !S integer :: iCol ! Counter for column loops !S integer :: iter ! Counter for iteration loop !S !S integer :: nstep ! current timestep number !S !S integer :: indxTe ! pbuf index for electron temperature !S integer :: indxTi ! pbuf index for ion temperature !S !S integer, parameter :: maxIter = 6 ! maximum number of iterations to solve for electron/ion temperature !S !S real(r8), parameter :: Kec1 = 7.5E5_r8 ! c1 constant for calculation of electron conductivity(Ke) !S real(r8), parameter :: Kec2 = 3.22E4_r8 ! c2 constant for calculation of electron conductivity(Ke) !S real(r8), parameter :: stepweight = 1.0_r8 ! weight of previous and current times step for diagonals !S !S real(r8) :: wrk1 ! 2/3/kboltz_ev !S real(r8) :: sqrt2 !S !S real(r8) :: lossc5 ! c5 prime of Lc(eN2) component of loss term !S real(r8) :: lossc7 ! c7 of Lc(eO2) component of loss term equation !S !S real(r8) :: lossc9 ! c9 of Lc(eO) component of loss term equation !S real(r8) :: lossc13 ! c13 of Lc(eO)f component of loss term equation !S real(r8) :: FeDB ! B term of electron heat flux of UB !S real(r8) :: FeD ! Day time flux !S real(r8) :: FeN ! Night time flux !S real(r8) :: f1Ted1 ! d1 of f1(Te) calculation used to get electron conductivity !S real(r8) :: f1Ted2 ! d2 of f1(Te) calculation used to get electron conductivity !S real(r8) :: f1Ted3 ! d3 of f1(Te) calculation used to get electron conductivity !S real(r8) :: f1Te !S !S real(r8), dimension(:,:), pointer :: tE ! Pointer to electron temperature in pbuf (K) !S real(r8), dimension(:,:), pointer :: ti ! Pointer to ion temperature in pbuf (K) !S !S real(r8), dimension(pcols,pver) :: pMid ! Midpoint pressure (Pa) !S real(r8), dimension(pcols,pver) :: tN ! Neutral temperature (K) !S real(r8), dimension(pcols,pver) :: tE0 ! Electron temperature from last time step (K) !S real(r8), dimension(pcols,pver) :: tEPrevI ! Electron temperature from previous iteration (K) !S !S real(r8), dimension(pcols,pverp) :: pInt ! Interface pressure (Pa) !S real(r8), dimension(pcols,pverp) :: tNInt ! Interface Temperture (K) !S real(r8), dimension(pcols,pverp) :: rairvi ! Constituent dependent gas constant on interface levels !S !S real(r8), dimension(pcols,pver) :: ndensN2 ! N2 number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensO2 ! O2 number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensO1 ! O number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensNO ! NO number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensN1 ! N number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensE ! E electron number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensOp ! O plus number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensO2p ! O2 plus ion number density (cm-3) !S real(r8), dimension(pcols,pver) :: ndensNOp ! NO plus ion number density (cm-3) !S !S real(r8), dimension(pcols,pver) :: dipMag ! dip angle for each column (radians) !S real(r8), dimension(pcols,pver) :: dipMagD ! dip angle for each column (degrees) !S !S real(r8), dimension(ncol) :: zenAngD ! zenith angle (degrees) !S !S real(r8), dimension(pcols) :: Feub ! electron heat flux at upper boundary !S !S real(r8), dimension(pcols,pver) :: sqrtTE ! Square root of electron temperature !S !S real(r8), dimension(pcols,pver) :: sqrtTN ! Square root of electron temperature !S !S real(r8), dimension(pver) :: Ke ! electron conductivity !S !S real(r8), dimension(pverp) :: Kei ! electron conductivity interface levels !S !S real(r8), dimension(pcols,pver) :: lossc4p ! c4 prime of Lc(eN2) component of loss term !S real(r8), dimension(pcols,pver) :: lossceN2 ! Lc(eN2) component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc6p ! c6 prime of Lc(eO2) component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceO2 ! Lc(eO2) component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc8p ! c8 prime of Lc(eO) component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceO1 ! Lc(eO) component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc10p ! c10 prime of Lc(eN2) component of loss term equation !S real(r8), dimension(pcols,pver) :: losscA ! A of Lc(eN2)v component of loss term equation !S real(r8), dimension(pcols,pver) :: tENDiff ! Difference between electron and neutral temperatures !S real(r8), dimension(pcols,pver) :: lossceN2v ! Lc(eN2)v component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc11p ! c11 prime of Lc(eO2)v component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceO2v ! Lc(eO2)v component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc12p ! c12 prime of Lc(eO)f component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceOf ! Lc(eO)f component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc14p ! c14 prime of Lc(eO)1D component of loss term equation !S real(r8), dimension(pcols,pver) :: losscf2d ! d of f2 of Lc(eO)1D component of loss term equation !S real(r8), dimension(pcols,pver) :: losscf2 ! f2 of Lc(eO)1D component of loss term equation !S real(r8), dimension(pcols,pver) :: losscf3 ! f3 of Lc(eO)1D component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceO1D ! Lc(eO)1D component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc15p ! c15 prime of Lc(eN2)Rot component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceN2Rot ! Lc(eN2)Rot component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc16p ! c16 prime of Lc(eO2)Rot component of loss term equation !S real(r8), dimension(pcols,pver) :: lossceO2Rot ! Lc(eO2)Rot component of loss term equation !S !S real(r8), dimension(pcols,pver) :: lossc3p ! c3 prime of Lc(ei) component of loss term equation !S real(r8), dimension(pcols,pver) :: losscei ! Lc(ei) component of loss term equation !S real(r8), dimension(pcols,pver) :: losscin ! ion-neutral heating coeff. !S !S real(r8), dimension(pcols,pver) :: lossg3 ! g3 loss term for Te tendency !S !S real(r8), dimension(pcols,pver) :: delTEN ! Difference between electron and neutral temperatures from production/loss !S !S real(r8), dimension(pcols,pverp) :: delZi ! Delta z: interfaces !S real(r8), dimension(pcols,pver) :: delZ ! Delta z: midpoints !S !S real(r8), dimension(pcols,pver) :: sourceg4 ! g4 source term for electron/ion temperature update !S !S !S real(r8), dimension(pcols,pver) :: qjoule ! joule heating !S real(r8), dimension(pcols,pver) :: qen ! electron-neutral heating !S real(r8), dimension(pcols,pver) :: qei ! electron-ion Coulomb heating !S real(r8), dimension(pcols,pver) :: rho ! mass density !S !S real(r8), dimension(pcols,pver) :: wrk2 !S !S real(r8), dimension(ionBot) :: subdiag ! subdiagonal values for Te tendency solving !S real(r8), dimension(ionBot) :: superdiag ! superdiagonal values for Te tendency solving !S real(r8), dimension(ionBot) :: diag ! diagonal values for Te tendency solving !S real(r8), dimension(ionBot) :: rHSInit ! initial RHS of electron temperature update !S real(r8), dimension(ionBot) :: rHSH ! h for RHS of electron temperature update !S real(r8), dimension(ionBot) :: rHS ! RHS of electron temperature update !S real(r8), dimension(ionBot) :: tETemp ! temporary electron temperature array for input to tridag !S !S logical, dimension(pcols) :: colConv ! flag for column converging = 1 if converged otherwise = 0 !S !S logical :: converged !Flag for convergence in electron temperature calculation iteration loop !S !S!----------------------------------------------------------------------------------------------------------------- !S !S!--------------------------------------------------------------------------------------------------------- !S! Initialize arrays to zero !S!--------------------------------------------------------------------------------------------------------- !S pMid(:,:) = 0._r8 !S tN(:,:) = 0._r8 !S pInt(:,:) = 0._r8 !S tNInt(:,:) = 0._r8 !S rairvi(:,:) = 0._r8 !S ndensO2(:,:) = 0._r8 !S ndensO1(:,:) = 0._r8 !S ndensNO(:,:) = 0._r8 !S ndensN1(:,:) = 0._r8 !S ndensE(:,:) = 0._r8 !S ndensOp(:,:) = 0._r8 !S ndensO2p(:,:) = 0._r8 !S ndensNOp(:,:) = 0._r8 !S ndensN2(:,:) = 0._r8 !S zenAngD(:) = 0._r8 !S sqrtTE(:,:) = 0._r8 !S sqrtTN(:,:) = 0._r8 !S Ke(:) = 0._r8 !S Kei(:) = 0._r8 !S lossc4p(:,:) = 0._r8 !S lossceN2(:,:) = 0._r8 !S lossc6p(:,:) = 0._r8 !S lossc6p(:,:) = 0._r8 !S lossceO2(:,:) = 0._r8 !S lossc8p(:,:) = 0._r8 !S lossceO1(:,:) = 0._r8 !S lossc10p(:,:) = 0._r8 !S losscA(:,:) = 0._r8 !S tENDiff(:,:) = 0._r8 !S lossceN2v(:,:) = 0._r8 !S lossc11p(:,:) = 0._r8 !S lossceO2v(:,:) = 0._r8 !S lossc12p(:,:) = 0._r8 !S lossceOf(:,:) = 0._r8 !S lossc14p(:,:) = 0._r8 !S losscf2d(:,:) = 0._r8 !S losscf2(:,:) = 0._r8 !S losscf3(:,:) = 0._r8 !S lossceO1D(:,:) = 0._r8 !S lossc15p(:,:) = 0._r8 !S lossceN2Rot(:,:) = 0._r8 !S lossc16p(:,:) = 0._r8 !S lossceO2Rot(:,:) = 0._r8 !S lossc3p(:,:) = 0._r8 !S losscei(:,:) = 0._r8 !S losscin(:,:) = 0._r8 !S lossg3(:,:) = 0._r8 !S delTEN(:,:) = 0._r8 !S delZi(:,:) = 0._r8 !S delZ(:,:) = 0._r8 !S subDiag(:) = 0._r8 !S superDiag(:) = 0._r8 !S diag(:) = 0._r8 !S sourceg4(:,:) = 0._r8 !S dipMag(:,:) = 0._r8 !S dipMagD(:,:) = 0._r8 !S rHSInit(:) = 0._r8 !S rHSH(:) = 0._r8 !S rHS(:) = 0._r8 !S teTemp(:) = 0._r8 !S qjoule(:,:) = 0._r8 !S qei(:,:) = 0._r8 !S qen(:,:) = 0._r8 !S rho(:,:) = 0._r8 !S colConv(:) = .false. !S !S!------------------------------------------- !S! Calculate some commonly used variables !S!------------------------------------------- !S wrk1 = 2._r8 / 3._r8/ kboltz_ev !S sqrt2 = SQRT(2._r8) !S ionBotP = ionBot + 1 !S !S!------------------------------------------------------------------------------------------------------- !S! Need to get midpoint and interface pressure and neutral temperature from state structure (ncol,ionBot) !S!------------------------------------------------------------------------------------------------------- !S pMid(1:ncol,1:pver) = state%pmid(1:ncol,1:pver) !S tN(1:ncol,1:pver) = state%t(1:ncol,1:pver) !S rho(1:ncol,1:pver) = pMid(1:ncol,1:pver)/rairv(1:ncol,1:pver,lchnk)/tN(1:ncol,1:pver) * 1.E-3_r8 ! convert to g/cm3 !S !S qjoule(1:ncol,1:ionBot) = ptend%s(1:ncol,1:ionBot) * 6.24E15_r8 ! convert from J/kg/s to ev/g/s !S !S pInt(1:ncol,1:pverp) = state%pint(1:ncol,1:pverp) !S tNInt(1:ncol,1:pverp) = istate%tNInt(1:ncol,1:pverp) !S rairvi(1:ncol,1:pverp) = istate%rairvi(1:ncol,1:pverp) !S !S sqrtTN(1:ncol,1:pver) = SQRT(tN(1:ncol,1:pver)) !S !S!---------------------------------------------------------------- !S! Get variables needed from the ionosphere state structure !S!---------------------------------------------------------------- !S ndensO2(1:ncol,1:pver) = istate%ndensO2(1:ncol,1:pver) !S ndensO1(1:ncol,1:pver) = istate%ndensO1(1:ncol,1:pver) !S ndensNO(1:ncol,1:pver) = istate%ndensNO(1:ncol,1:pver) !S ndensN1(1:ncol,1:pver) = istate%ndensN1(1:ncol,1:pver) !S ndensE(1:ncol,1:pver) = istate%ndensE(1:ncol,1:pver) !S ndensOp(1:ncol,1:pver) = istate%ndensOp(1:ncol,1:pver) !S ndensO2p(1:ncol,1:pver) = istate%ndensO2p(1:ncol,1:pver) !S ndensNOp(1:ncol,1:pver) = istate%ndensNOp(1:ncol,1:pver) !S ndensN2(1:ncol,1:pver) = istate%ndensN2(1:ncol,1:pver) !S !S sourceg4(1:ncol,1:pver) = istate%sourceg4(1:ncol,1:pver) !S !S dipMag(1:ncol,1:pver) = istate%dipMag(1:ncol,1:pver) !S dipMagD(1:ncol,1:pver) = istate%dipMagD(1:ncol,1:pver) !S !S zenAngD(1:ncol) = istate%zenAngD(1:ncol) !S !S!------------------------------------------------------------------------------------------------------------------- !S! Get electron temperature from physics buffer and if this is the first time calculated then initialize to neutral !S! temperature. nsrest=0 means this is an initial run and nstep=0 means this is the first time step. !S!------------------------------------------------------------------------------------------------------------------- !S indxTe = pbuf_get_index( 'ElecTemp' ) !S call pbuf_get_field(pbuf, indxTe, tE) !S !S indxTi = pbuf_get_index( 'IonTemp' ) !S call pbuf_get_field(pbuf, indxTi, ti) !S !S nstep = get_nstep() !S !S if(nsrest == 0 .and. nstep == 0) then !S !S ti(1:ncol,1:pver) = tN(1:ncol,1:pver) !S tE(1:ncol,1:pver) = tN(1:ncol,1:pver) !S !S else !S !S tE(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),tE(1:ncol,1:pver)) !S tE(1:ncol,1:pver) = MIN(temax,tE(1:ncol,1:pver)) !S !S ti(1:ncol,1:pver) = MAX(tN(1:ncol,1:pver),ti(1:ncol,1:pver)) !S ti(1:ncol,1:pver) = MIN(ti(1:ncol,1:pver),tE(1:ncol,1:pver)) !S !S tE(1:ncol,ionBotP:pver) = tN(1:ncol,ionBotP:pver) !S ti(1:ncol,ionBotP:pver) = tN(1:ncol,ionBotP:pver) !S !S endif !S !S tE0(1:ncol,1:pver) = tE(1:ncol,1:pver) !S !S wrk2(1:ncol,1:ionBot) = ndensE(1:ncol,1:ionBot)/wrk1/(SIN(dipMag(1:ncol,1:ionBot)))**2._r8 !S !S!----------------------------------------------------------------------------- !S! Get constant terms needed for loss term g3 for electron temperature update !S!----------------------------------------------------------------------------- !S lossc5 = 1.21E-4_r8 !S lossc7 = 3.6E-2_r8 !S lossc9 = 5.7E-4_r8 !S lossc13 = 7.E-5_r8 !S !S!----------------------------------------------------------------------------- !S! Get terms needed for loss term g3 for electron temperature update which do !S! not need to be updated in iteration loop. !S!----------------------------------------------------------------------------- !S do iCol = 1, ncol !S !S if (.not. colConv(iCol)) then !S !S do iVer = 1, ionBot !S !S lossc4p(iCol,iVer) = 1.77E-19_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) !S lossc6p(iCol,iVer) = 1.21E-18_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) !S lossc8p(iCol,iVer) = 7.9E-19_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) !S lossc10p(iCol,iVer) = 1.3E-4_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) !S lossc11p(iCol,iVer) = 3.125E-21_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) !S lossc12p(iCol,iVer) = 3.4E-12_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) !S lossc14p(iCol,iVer) = 1.57E-12_r8 * ndensO1(iCol,iVer) * ndensE(iCol,iVer) !S lossc15p(iCol,iVer) = 2.9E-14_r8 * ndensN2(iCol,iVer) * ndensE(iCol,iVer) !S lossc16p(iCol,iVer) = 6.9E-14_r8 * ndensO2(iCol,iVer) * ndensE(iCol,iVer) !S lossc3p(iCol,iVer) = 3.2E-8_r8 * 15._r8 * (ndensOP(iCol,iVer) + & !S 0.5_r8 * ndensO2P(iCol,iVer) + 0.53_r8 * ndensNOP(iCol,iVer)) * ndensE(iCol,iVer) !S !S losscin(iCol,iVer) = (6.6e-14_r8*ndensN2(iCol,iVer) + 5.8e-14_r8*ndensO2(iCol,iVer) & !S + 0.21e-14_r8*ndensO1(iCol,iVer)*sqrt2*sqrtTN(iCol,iVer))*ndensOP(iCol,iVer) & !S +(5.9e-14_r8*ndensN2(iCol,iVer) + 5.45e-14_r8*ndensO2(iCol,iVer) & !S + 4.5e-14_r8*ndensO1(iCol,iVer))*ndensOP(iCol,iVer) & !S +(5.8e-14_r8*ndensN2(iCol,iVer) + 0.14e-14_r8*ndensO2(iCol,iVer)*sqrtTN(iCol,iVer) & !S + 4.4e-14_r8*ndensO1(iCol,iVer)) * ndensO2P(iCol,iVer) !S !S enddo !iVer loop !S !S!---------------------------------------------------------------------------------- !S! Calculate upper boundary heat flux !S!---------------------------------------------------------------------------------- !S if (ABS(dipMagD(iCol,1)) < 60.0_r8) FeDB = 0.5_r8 * & !S (1._r8 + SIN(pi * (ABS(dipMagD(iCol,1)) - 30.0_r8) /60.0_r8)) !S !S if (ABS(dipMagD(iCol,1)) >= 60.0_r8) FeDB = 1._r8 !S !S FeD = -4.5E-7_r8 * f107 * 0.5_r8 * FeDB !S FeN = .2_r8 * FeD !S !S!--------------------------------------------------- !S! Set upper boundary condition for right hand side !S!--------------------------------------------------- !S if (zenAngD(iCol) <= 80.0_r8) Feub(iCol) = FeD !S if (zenAngD(iCol) > 80.0_r8 .AND. zenAngD(iCol) < 100.0_r8) Feub(iCol) = 0.5_r8 * (FeD + FeN) & !S + 0.5_r8 * (FeD - FeN) * & !S COS(pi * ((zenAngD(iCol) - 80.0_r8) / 20.0_r8)) !S if (zenAngD(iCol) >= 100.0_r8) Feub(iCol) = FeN !S !S!------------------------------------------------------------------------------------------ !S! Calculate thickness terms for vertical derivative !S!------------------------------------------------------------------------------------------ !S do iVer = 1, ionBot !S !S delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rairv(iCol,iVer,lchnk) * tN(iCol,iVer) / pMid(iCol,iVer) / gravit !S !S enddo !S !S do iVer = 2, ionBotP ! Assuming ionBotP < pverp !S delZi(iCol,iVer) = (pMid(iCol,iVer) - pMid(iCol,iVer-1)) * rairvi(iCol,iVer) * & !S tNInt(iCol,iVer) / pInt(iCol,iVer) / gravit !S enddo !S delZi(iCol,1) = 1.5_r8*delZi(iCol,2) - .5_r8*delZi(iCol,3) !S !S!---------------------------------------------------------- !S! Convert delZ variables from meters to centimeters !S!---------------------------------------------------------- !S delZi(iCol,1:ionBotP) = delZi(iCol,1:ionBotP)*100._r8 !S delZ(iCol,1:ionBot) = delZ(iCol,1:ionBot)*100._r8 !S !S endif ! Column not converged !S !S enddo !iCol loop !S !S!------------------------------------------------------------------------------------------------------- !S! Iterate to calculate new electron temperature. !S! Time splitting is used: first solve the heating/cooling equation, then solve the diffusion equations. !S! Also, set convergence flag to false and iterate until true or 6 iterations, whichever comes first !S!------------------------------------------------------------------------------------------------------- !S converged = .false. !S iter = 0 !S do while (.not. converged .and. iter < maxIter) !S !S!-------------------------------------------------------------------------------------------------------- !S! Increment iteration loop counter and save electron temperature from previous iteration for convergence !S! test at end of interation loop. Also, take square root of electron temperature to be used later !S!-------------------------------------------------------------------------------------------------------- !S iter = iter + 1 !S !S tEPrevI(1:ncol,1:ionBot) = tE(1:ncol,1:ionBot) !S !S sqrtTE(1:ncol,1:ionBot) = SQRT(tE(1:ncol,1:ionBot)) !S !S!-------------------------------------------------------------------------------------------------------- !S! Loop over columns then vertical levels and call tridiagonal solver for each column to get electron !S! temperature !S!-------------------------------------------------------------------------------------------------------- !S do iCol = 1, ncol !S !S if (.not. colConv(iCol)) then !S !S do iVer = 1, ionBot !S !S!----------------------------------------------------------------------------- !S! Get loss term g3 for electron temperature update. Need to calculate !S! constituent dependent loss terms which make up g3 !S!----------------------------------------------------------------------------- !S lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer) !S lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iCol,iVer)) * sqrtTE(iCol,iVer) !S lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iCol,iVer) !S !S if (tE(iCol,iVer) < 1000.0_r8) losscA(iCol,iVer) = 5.71E-8_r8 * EXP(-3352.6_r8 / tE(iCol,iVer)) !S if (tE(iCol,iVer) >= 1000.0_r8 .AND. tE(iCol,iVer) <= 2000.0_r8) & !S losscA(iCol,iVer) = 2.0E-7_r8 * EXP(-4605.2_r8 / tE(iCol,iVer)) !S if (tE(iCol,iVer) > 2000.0_r8) losscA(iCol,iVer) = 2.53E-6_r8 * sqrtTE(iCol,iVer) * & !S EXP(-17620._r8 / tE(iCol,iVer)) !S !S tENDiff(iCol,iVer) = tE(iCol,iVer) - tN(iCol,iVer) !S !S if (ABS(tENDiff(iCol,iVer)) < 0.1_r8) tENDiff(iCol,iVer) = 0.1_r8 !S !S lossceN2v(iCol,iVer) = lossc10p(iCol,iVer) * losscA(iCol,iVer) * & !S (1._r8 - EXP(3200._r8 * (1._r8 / tE(iCol,iVer) - 1._r8 / tN(iCol,iVer)))) / tENDiff(iCol,iVer) !S lossceO2v(iCol,iVer) = lossc11p(iCol,iVer) * tE(iCol,iVer)**2 !S lossceOf(iCol,iVer) = lossc12p(iCol,iVer) * (1._r8 - lossc13 * tE(iCol,iVer)) * & !S (0.4_r8 + 150._r8 / tE(iCol,iVer)) / tN(iCol,iVer) !S losscf2d(iCol,iVer) = 2.4E+4_r8 + 0.3_r8 * (tE(iCol,iVer) - 1500._r8) + & !S 1.947E-5_r8 * (tE(iCol,iVer) - 1500._r8) * (tE(iCol,iVer) - 4000._r8) !S losscf2(iCol,iVer) = losscf2d(iCol,iVer) * (1._r8 / 3000._r8 - 1._r8 / tE(iCol,iVer)) !S losscf3(iCol,iVer) = -22713._r8 * (1._r8 / tN(iCol,iVer) - 1._r8 / tE(iCol,iVer)) !S lossceO1D(iCol,iVer) = lossc14p(iCol,iVer) * EXP(losscf2(iCol,iVer)) * & !S (1._r8 - EXP(losscf3(iCol,iVer))) / tENDiff(iCol,iVer) !S lossceN2Rot(iCol,iVer) = lossc15p(iCol,iVer) / sqrtTE(iCol,iVer) !S lossceO2Rot(iCol,iVer) = lossc16p(iCol,iVer) / sqrtTE(iCol,iVer) !S losscei(iCol,iVer) = lossc3p(iCol,iVer) / tE(iCol,iVer)**1.5_r8 !S lossg3(iCol,iVer) = lossceN2(iCol,iVer) + lossceO2(iCol,iVer) + lossceO1(iCol,iVer) + lossceN2v(iCol,iVer) & !S + lossceO2v(iCol,iVer) + lossceOf(iCol,iVer) + lossceO1D(iCol,iVer) & !S + lossceN2Rot(iCol,iVer) + lossceO2Rot(iCol,iVer) !S !S tE(iCol,iVer) = (wrk2(iCol,iVer)/ztodt * tE0(iCol,iVer) + (lossg3(iCol,iVer) * tN(iCol,iVer) & !S + losscei(iCol,iVer) * ti(iCol,iVer) + sourceg4(iCol,iVer)) & !S /(SIN(dipMag(iCol,iVer)))**2) / & !S (wrk2(iCol,iVer)/ztodt + (lossg3(iCol,iVer) + losscei(iCol,iVer))/(SIN(dipMag(iCol,iVer)))**2._r8) !S !S!------------------------------------------------------------------------------------------------------------------- !S! Limit maximum value of electron temperature to be a locally set maximum (temax) and minimum value of electron !S! temperature to be at least the neutral temperature for a given column and level !S!------------------------------------------------------------------------------------------------------------------- !S tE(iCol,iVer) = min(temax,tE(iCol,iVer)) !S tE(iCol,iVer) = max(tN(iCol,iVer),tE(iCol,iVer)) !S !S enddo !iVer loop !S !S endif ! Column not converged !S !S enddo ! End of column loop !S !S! fvitt -- make sqrtTE consistent with new tE !S sqrtTE(1:ncol,1:ionBot) = SQRT(tE(1:ncol,1:ionBot)) !S !S!----------------------------------------------------- !S! Calculate thermal conductivity of electron gas !S!----------------------------------------------------- !S do iCol = 1, ncol !S !S if (.not. colConv(iCol)) then !S !S do iVer = 1, ionBot !S !S f1Ted1 = 2.82E-17_r8 * sqrtTE(iCol,iVer) - 3.41E-21_r8 * tE(iCol,iVer)**1.5_r8 !S f1Ted2 = 2.2E-16_r8 + 7.92E-18_r8 * sqrtTE(iCol,iVer) !S f1Ted3 = 1.1E-16_r8 * (1._r8 + 5.7E-4_r8 * tE(iCol,iVer)) !S !S f1Te = ndensN2(iCol,iVer) / ndensE(iCol,iVer) * f1Ted1 + ndensO2(iCol,iVer) / & !S ndensE(iCol,iVer) * f1Ted2 + ndensO1(iCol,iVer) / ndensE(iCol,iVer) * f1Ted3 !S !S!----------------------------------------------------------------------------- !S! Calculate electron conductivity using parameters set in module and f1(Te) !S!----------------------------------------------------------------------------- !S Ke(iVer) = Kec1 * tE(iCol,iVer)**2.5_r8 / (1._r8 + Kec2 * tE(iCol,iVer)**2._r8 * f1Te) !S !S enddo !iVer loop !S !S!---------------------------------------------------------------------- !S! Get electron conductivity at interface levels to be used later !S!---------------------------------------------------------------------- !S do iVer = 2,ionBot !S Kei(iVer) = SQRT(Ke(iVer-1)*Ke(iVer)) !S enddo !S Kei(1) = 1.5_r8*Kei(2)-.5_r8*Kei(3) !S Kei(ionBotP) = 1.5_r8*Kei(ionBot)-.5_r8*Kei(ionBot-1) !S !S!------------------------------------------------------------------------------------------------------ !S! Derive subdiagonal, superdiagonal, and diagonal as input to solver for electron temperature tendency !S!------------------------------------------------------------------------------------------------------ !S do iVer = 2, ionBot-1 !S subDiag(iVer) = -Kei(iVer) / delZi(iCol,iVer) / delZ(iCol,iVer) !S superDiag(iVer) = -Kei(iVer+1) / delZi(iCol,iVer+1) / delZ(iCol,iVer) !S diag(iVer) = wrk2(iCol,iVer)/ztodt-subDiag(iVer)-superDiag(iVer) !S rHS(iVer) = tE(iCol,iVer) * wrk2(iCol,iVer)/ztodt !S enddo !iVer loop !S !S!------------------------------------------------------------------------------------- !S! Calculate diagonal, superdiagonal, and right hand side upper boundary values !S!------------------------------------------------------------------------------------- !S superDiag(1) = -Kei(2) / delZi(iCol,2) / delZ(iCol,1) !S diag(1) = wrk2(iCol,1)/ztodt - superDiag(1) !S rHS(1) = tE(iCol,1) * wrk2(iCol,1) / ztodt - Feub(iCol) / delZ(iCol,1) !S !S!--------------------------------------------------------------------------------------------- !S! Calculate subdiagonal, diagonal, superdiagonal, and right hand side lower boundary values !S!--------------------------------------------------------------------------------------------- !S subDiag(ionBot) = -Kei(ionBot) / delZi(iCol,ionBot) / delZ(iCol,ionBot) !S superDiag(ionBot) = -Kei(ionBotP) / delZi(iCol,ionBotP) / delZ(iCol,ionBot) !S diag(ionBot) = wrk2(iCol,ionBot)/ztodt-subDiag(ionBot)-superDiag(ionBot) !S rHS(ionBot) = tE(iCol,ionBot) * wrk2(iCol,ionBot)/ztodt -superDiag(ionBot)*tN(iCol,ionBotP) !S !S!------------------------------------------------- !S! Call solver to get electron temperature update !S!------------------------------------------------- !S call tridag(subDiag,diag,superDiag,rHS,tETemp,ionBot) !S !S tE(iCol,1:ionBot) = tETemp(1:ionBot) !S do iVer = 1,ionBot !S tE(iCol,iVer) = min(temax,tE(iCol,iVer)) !S tE(iCol,iVer) = max(tN(iCol,iVer),tE(iCol,iVer)) !S enddo !S!--------------------------------------------------------------------------------------------------------- !S! Calculate ion temperature from electron temperature, ion-neutral and electron-ion loss terms, neutral !S! temperature, mass density and joule heating. Set minimum value to neutral temperature and maximum !S! value to electron temperature for each column and vertical level !S!--------------------------------------------------------------------------------------------------------- !S do iVer = 1,ionBot !S ti(iCol,iVer) = (losscei(iCol,iVer) * tE(iCol,iVer) + losscin(iCol,iVer) * tN(iCol,iVer) + & !S rho(iCol,iVer) * qjoule(iCol,iVer))/(losscei(iCol,iVer) + losscin(iCol,iVer)) !S ti(iCol,iVer) = max(tN(iCol,iVer),ti(iCol,iVer)) !S ti(iCol,iVer) = min(tE(iCol,iVer),ti(iCol,iVer)) !S enddo !S !S!-------------------------------------------------------------------------------------------------------- !S! Check for convergence which is a change of electron temperature ratio to previous loop for all levels !S! and columns of less than 0.05K. Had to modify this to do convergence check on each column since !S! checking all columns in a chunk gives different answers depending on number of tasks and tasks per node. !S!-------------------------------------------------------------------------------------------------------- !S if (ALL(ABS(tE(iCol,1:ionBot) / tEPrevI(iCol,1:ionBot) - 1._r8) < 0.05_r8)) then !S !S colConv(iCol) = .true. !S !S endif !S !S endif ! Column not converged !S !S enddo ! iCol loop !S!-------------------------------------------------------------- !S! Check to see if all columns have converged and set flag !S!-------------------------------------------------------------- !S if (ALL(colConv(1:ncol))) converged = .true. !S !S enddo ! End of iteration loop !S !S write(iulog,*)'lchnk, Number of tempei iterations, converged flag ', lchnk, iter, converged !S !S!-------------------------------------------------------------------------------------------------------- !S! Calculate electron-neutral heating and electron-ion Coulomb heating. Then update dry static energy. !S!-------------------------------------------------------------------------------------------------------- !S do iVer = 1, ionBot !S do iCol = 1, ncol !S sqrtTE(iCol,iVer) = SQRT(tE(iCol,iVer)) !S lossceN2(iCol,iVer) = lossc4p(iCol,iVer) * (1._r8 - lossc5 * tE(iCol,iVer)) * tE(iCol,iVer) !S lossceO2(iCol,iVer) = lossc6p(iCol,iVer) * (1._r8 + lossc7 * sqrtTE(iCol,iVer)) * sqrtTE(iCol,iVer) !S lossceO1(iCol,iVer) = lossc8p(iCol,iVer) * (1._r8 + lossc9 * tE(iCol,iVer)) * sqrtTE(iCol,iVer) !S enddo !S enddo !S qen(1:ncol,1:ionBot) = (lossceN2(1:ncol,1:ionBot)+lossceO2(1:ncol,1:ionBot)+lossceO1(1:ncol,1:ionBot)) * & !S (tE(1:ncol,1:ionBot)-tN(1:ncol,1:ionBot)) / rho(1:ncol,1:ionBot) !S 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) !S 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 !S !S call outfld ('TE' , tE , pcols, lchnk) !S call outfld ('TI' , ti , pcols, lchnk) !S call outfld ('LOSS_g3' , lossg3 , pcols, lchnk) !S call outfld ('LOSS_EI' , losscei , pcols, lchnk) !S call outfld ('LOSS_IN' , losscin , pcols, lchnk) !S call outfld ('QIN' , qei , pcols, lchnk) !S !S return !S !S end subroutine ionos_tempei !=============================================================================== subroutine ionos_optransvert(ncol, istate, oPBot) !----------------------------------------------------------------------- ! Routine to compute O+ ion transport in the vertical direction through ! ambipolar diffusion. !----------------------------------------------------------------------- !S use cam_history, only : addfld, add_default, phys_decomp !S use phys_grid, only : get_lat_p, get_lon_p, get_rlat_p,get_rlon_p !S use physconst, only : gravit, avogad ! Gravity (m/s2) and Avogadro's number (molecules/kmole) !S use cam_control_mod, only : nsrest ! Flag to determine if this is an initial run or a restart/branch !S use time_manager, only : get_nstep ! Routine to get current time step !S use physconst, only : rairv, mbarv ! Constituent dependent rair and mbar !S use efield, only : ed1, ed2 ! x and y components of the electric field !S use mo_apex, only : bnorth, beast, bmag, alatm ! Magnetic field components and magnitude (nT) and geomagnetic latitude (radians) !S use exbdrift, only : map_mag2geo, cal_mlt ! Methods to convert e field from mag to geo coords !S use constituents, only : cnst_get_ind, cnst_mw ! Routines to get molecular weights for constituents !S use short_lived_species, only : slvd_index,slvd_pbf_ndx => pbf_idx ! Routines to access short lived species in storage and pbuf !S use charge_neutrality, only : charge_fix 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-------------------------------- !S type(physics_buffer_desc), pointer :: pbuf(:) ! physics buffer !S type(physics_state), intent(inout) :: state ! physics state structure !S type(physics_ptend), intent(inout) :: ptend ! parameterization tendency structure type(ionos_state), intent(inout) :: istate ! ionosphere state structure !S real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) !S integer, intent(in) :: lchnk ! chunk identifier 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,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 (?) !S real(r8), dimension(pcols) :: ed1_geo ! zonal electric field on geographic grid [V/m] !S real(r8), dimension(pcols) :: ed2_geo ! meridional electric field on geographic grid [V/m] 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] !S real(r8), dimension(:,:), pointer :: tE ! Pointer to electron temperature in pbuf (K) !S real(r8), dimension(:,:), pointer :: tI ! Pointer to ion temperature in pbuf (K) 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 !S real(r8), dimension(pcols,pver) :: mmrOp ! O plus ion mass mixing ratio kg/kg real(r8), dimension(pcols,pver) :: tP ! Plasma temperature (K) !S real(r8), dimension(pcols,pver) :: pMid ! Midpoint pressure (Pa) !S real(r8), dimension(pcols,pver) :: tN ! Neutral temperature (K) !S real(r8), dimension(pcols,pver) :: uN ! Neutral zonal wind velocity (m/s) !S real(r8), dimension(pcols,pver) :: vN ! Neutral meridional wind velocity (m/s) !S real(r8), dimension(pcols,pverp) :: pInt ! Interface pressure (Pa) 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) :: rairvi ! Constituent dependent gas constant on interface levels !S real(r8), dimension(pcols,pver) :: omega ! Vertical velocity (Pa/s?) 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) !S real(r8), dimension(pcols,pver) :: ndensOp ! 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+ !----------------------------------------------------------------------------------------------------------------- write(iulog,*)'lchnk, Inside optransvert ' !--------------------------------------------------------------------------------------------------------- ! Initialize arrays to zero !--------------------------------------------------------------------------------------------------------- !S ed1_geo(:) = 0._r8 !S ed2_geo(:) = 0._r8 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 !S pMid(:,:) = 0._r8 !S tN(:,:) = 0._r8 !S uN(:,:) = 0._r8 !S vN(:,:) = 0._r8 !S pInt(:,:) = 0._r8 tNInt(:,:) = 0._r8 tIInt(:,:) = 0._r8 tEInt(:,:) = 0._r8 rairvi(:,:) = 0._r8 !S omega(:,:) = 0._r8 wN(:,:) = 0._r8 ndensN2(:,:) = 0._r8 ndensO2(:,:) = 0._r8 ndensO1(:,:) = 0._r8 !S ndensOp(:,:) = 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 = 1._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) !S ndensOp(1:ncol,1:pver) = istate%ndensOp(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) zenAngD(1:ncol) = istate%zenAngD(1:ncol) !---------------------------------------------------------------- ! Get O+ mass !---------------------------------------------------------------- !S call cnst_get_ind( 'O', indxO1 ) !S rmassOP = cnst_mw(indxO1) rmassOP = 16._r8 !---------------------------------------------------------------- ! Get variables needed from the physics state structure !---------------------------------------------------------------- !S pInt(1:ncol,1:pverp) = state%pint(1:ncol,1:pverp) !S pMid(1:ncol,1:pver) = state%pmid(1:ncol,1:pver) !S omega(1:ncol,1:pver) = state%omega(1:ncol,1:pver) !S tN(1:ncol,1:pver) = state%t(1:ncol,1:pver) !S uN(1:ncol,1:pver) = state%u(1:ncol,1:pver) !S vN(1:ncol,1:pver) = state%v(1:ncol,1:pver) !------------------------------------------------------------------------------------------------------------------- ! 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. !------------------------------------------------------------------------------------------------------------------- !S indxTe = pbuf_get_index( 'ElecTemp' ) !S call pbuf_get_field(pbuf, indxTe, tE ) !S indxTi = pbuf_get_index( 'IonTemp' ) !S call pbuf_get_field(pbuf, indxTe, tI ) !S nstep = get_nstep() nstep = 1 !S if(nsrest == 0 .and. nstep == 0) then !S !S ti(1:ncol,1:pver) = tN(1:ncol,1:pver) !S tE(1:ncol,1:pver) = tN(1:ncol,1:pver) !S !S endif !S call cnst_get_ind( 'e', indxE, abort=.false. ) !S write(iulog,*)'optransvert lchnk,indxE ',lchnk,indxE !S !S if (indxE < 0) then !S sindxE = slvd_index( 'e' ) !S! write(iulog,*)'optransvert lchnk,sindxE,slvd_pbf_ndx ',lchnk,sindxE,slvd_pbf_ndx !S if (sindxE > 0) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrPE, start=(/1,1,sindxE/), kount=(/pcols,pver,1/) ) !S! write(iulog,*)'optransvert after pbuf_get_field E lchnk ',lchnk !S mmrE(1:ncol,1:pver) = mmrPE(1:ncol,1:pver) !S endif !S else !S mmrE(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxE) !S endif !S !S call cnst_get_ind( 'Op', indxOp, abort=.false. ) !S! write(iulog,*)'optransvert lchnk,indxOp ',lchnk,indxOp !S if (indxOp < 0) then !S sindxOp = slvd_index( 'Op' ) !S! write(iulog,*)'optransvert lchnk,sindxOp,slvd_pbf_ndx ',lchnk,sindxOp,slvd_pbf_ndx !S if (sindxOp > 0) then !S call pbuf_get_field(pbuf, slvd_pbf_ndx, mmrPOp, start=(/1,1,sindxOp/), kount=(/pcols,pver,1/) ) !S! write(iulog,*)'optransvert after pbuf_get_field Op lchnk ',lchnk !S mmrOp(1:ncol,1:pver) = mmrPOp(1:ncol,1:pver) !S endif !S else !S mmrOp(1:ncol,1:pver) = state%q(1:ncol,1:pver,indxOp) !S endif ! write(iulog,*)'lchnk, min/max after pbuf get mmrE(1:oPBot) ', lchnk, MINVAL(mmrE(1,1:oPBot)), MAXVAL(mmrE(1,1:oPBot)) ! write(iulog,*)'lchnk, min/max after pbuf get mmrOp(1:oPBot) ', lchnk, MINVAL(mmrOp(1,1:oPBot)), MAXVAL(mmrOp(1,1:oPBot)) !------------------------------------------------- ! Get geomagnetic latitudes for this chunk !------------------------------------------------- !S geoMagLat(1:ncol) = alatm(1:ncol,lchnk) !-------------------------------------------------- ! Calculate cosine/sine of declination angle !-------------------------------------------------- !S decMag(1:ncol) = -atan2( beast(1:ncol,lchnk),bnorth(1:ncol,lchnk) ) 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 !------------------------------------------------------------------------------------------------------------ !S call cal_mlt( mlt, lchnk, ncol ) !------------------------------------------------------------------------------------------------------------------ ! Convert e_field module zonal and meridional electric field components to geographic coords using exbdrift module ! map_mag2geo method !------------------------------------------------------------------------------------------------------------------ !S call map_mag2geo( mlt, lchnk, ncol, ed1_geo, ed2_geo, epot_geo ) !------------------------------------------------------------------------------------- ! 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 interface level !-------------------------- do iCol = 1, ncol tP(iCol,oPBotP) = 1.5_r8 * tP(iCol,oPBot) - 0.5_r8 * tP(iCol,oPBot-1) enddo !------------------------------------------------------------------------------------------ ! Calculate thickness terms for vertical derivative !------------------------------------------------------------------------------------------ do iVer = 1, oPBot do iCol = 1, ncol !S delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rairv(iCol,iVer,lchnk) * tN(iCol,iVer) / pMid(iCol,iVer) / gravit delZ(iCol,iVer) = (pInt(iCol,iVer+1) - pInt(iCol,iVer)) * rair * tN(iCol,iVer) / pMid(iCol,iVer) / gravity enddo enddo do iVer = 2, oPBotP ! Assuming oPBotP < pverp do iCol = 1, ncol !S delZi(iCol,iVer) = (pMid(iCol,iVer) - pMid(iCol,iVer-1)) * rairvi(iCol,iVer) * & !S tNInt(iCol,iVer) / pInt(iCol,iVer) / gravit 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 !----------------------------------------------------------------------------- ! 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) ! write(iulog,*)'lchnk, Before vertical loop optransvert ' !------------------------------------------------------------------ ! Calculate pressure scale height and vertical velocity !------------------------------------------------------------------ do iVer = 2,oPBot do iCol = 1,ncol !S pScaleHeight = .5_r8*(rairv(iCol,iVer,lchnk)*tN(iCol,iVer)+rairv(iCol,iVer-1,lchnk)*tN(iCol,iVer-1))/gravit pScaleHeight = .5_r8*(rair*tN(iCol,iVer)+rair*tN(iCol,iVer-1))/gravity ! if (iVer < 3) write(iulog,*)'lchnk,iVer,iCol,pScaleHeight,rairv(iCol,iVer,lchnk),tN(iCol,iVer),gravit ', lchnk,iVer,iCol,pScaleHeight,rairv(iCol,iVer,lchnk),tN(iCol,iVer),gravit wN(iCol,iVer) = -omega(iCol,iVer) * pScaleHeight ! if (iVer < 3) write(iulog,*)'lchnk,iVer,iCol,wN(iCol,iVer),state%omega(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight ', lchnk,iVer,iCol,wN(iCol,iVer),state%omega(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight write(iulog,*)'iVer,iCol,wN(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight,rair,tN(iCol,iVer),gravity ', iVer,iCol,wN(iCol,iVer),omega(iCol,iVer),pMid(iCol,iVer),pScaleHeight,rair,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 !----------------------------------------------------------------------------- !S bMagT = bmag(iCol,lchnk)*1.e-9_r8 bMagT = bmag*1.e-9_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) !--------------------------------------------------------------------------------- !S c1(iCol,iVer) = -eHPB(iCol) / bMagT * cos(dipMag(iCol,iVer)) * 100._r8 + & !S ( ( vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * sin(decMag(iCol)) ) * & !S cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ) * & !S sin(dipMag(iCol,iVer)) * 100._r8 + & !S dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + & !S 0.5_r8 * rmassOp / avogad * gravit / kboltz / tP(iCol,iVer) * .01_r8 ) * & !S 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 !-------------------------------------------------------------- ! Compute C2 coefficient for matrix formation ! dA(cm2/s) !-------------------------------------------------------------- c2(iCol,iVer) = dA(iCol,iVer) * sin(dipMag(iCol,iVer))**2 ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, c1(iCol,iVer),eHPB(iCol) / bMagT * cos(dipMag(iCol,iVer)) ', lchnk, iVer, iCol,c1(iCol,iVer), eHPB(iCol)/ bMagT * cos(dipMag(iCol,iVer)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, vN(iCol,iVer) * cos(decMag(iCol)) ', lchnk, iVer, iCol, (vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * cos(decMag(iCol))) * cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) * sin(dipMag(iCol,iVer)) * 100._r8 ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, vN(iCol,iVer) * cos(decMag(iCol)) ', lchnk, iVer, iCol, vN(iCol,iVer) * cos(decMag(iCol)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, uN(iCol,iVer) * cos(decMag(iCol)) ', lchnk, iVer, iCol, uN(iCol,iVer) * cos(decMag(iCol)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, cos(dipMag(iCol,iVer)) ', lchnk, iVer, iCol, cos(dipMag(iCol,iVer)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ', lchnk, iVer, iCol, wN(iCol,iVer) * sin(dipMag(iCol,iVer)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, wN(iCol,iVer) ', lchnk, iVer, iCol, wN(iCol,iVer) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, sin(dipMag(iCol,iVer)) ', lchnk, iVer, iCol, sin(dipMag(iCol,iVer)) ! if (iVer < 3) write(iulog,*)'lchnk, iVer, iCol, sin(dipMag(iCol,iVer)) * 100._r8 ', lchnk, iVer, iCol, sin(dipMag(iCol,iVer)) * 100._r8 ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, eHPB(iCol) / bMagT * cos(dipMag(iCol,iVer)) ', lchnk, iVer, iCol, eHPB(iCol)/ bMagT * cos(dipMag(iCol,iVer)) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, vN(iCol,iVer) * cos(decMag(iCol)) ', lchnk, iVer, iCol, (vN(iCol,iVer) * cos(decMag(iCol)) + uN(iCol,iVer) * cos(decMag(iCol))) * cos(dipMag(iCol,iVer)) - wN(iCol,iVer) * sin(dipMag(iCol,iVer)) * sin(dipMag(iCol,iVer)) * 100._r8 ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, decMag(iCol), uN(iCol,iVer) * cos(decMag(iCol)) ', lchnk, iVer, iCol, decMag(iCol), uN(iCol,iVer) * cos(decMag(iCol)) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, wN(iCol,iVer) ', lchnk, iVer, iCol, wN(iCol,iVer) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, dA(iCol,iVer) term ', lchnk, iVer, iCol, dA(iCol,iVer) * ( delTpDelZ / tP(iCol,iVer) + 0.5_r8 * rmassOp / avogad * gravit / kboltz / tP(iCol,iVer) * .01_r8 ) * sin(dipMag(iCol,iVer))**2 ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, avogad, dA(iCol,iVer) term pt1 ', lchnk, iVer, iCol, avogad, dA(iCol,iVer) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, dA(iCol,iVer) term pt2 ', lchnk, iVer, iCol, delTpDelZ / tP(iCol,iVer) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, dA(iCol,iVer) term pt3 ', lchnk, iVer, iCol, 0.5_r8 * rmassOp / avogad * gravit / kboltz / tP(iCol,iVer) * .01_r8 ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, dA(iCol,iVer) term pt4 ', lchnk, iVer, iCol, sin(dipMag(iCol,iVer))**2 ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, delTpDelZ, rmassOP, tP(iCol,iVer) ', lchnk, iVer, iCol, delTpDelZ, rmassOP, tP(iCol,iVer) ! if (c1(iCol,iVer) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, c1(iCol,iVer) ', lchnk, iVer, iCol, c1(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, eHPB(iCol), bMagT, dipMag(iCol,iVer) ', lchnk, iVer, iCol, eHPB(iCol), bMagT, dipMag(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, vN(iCol,iVer) ', lchnk, iVer, iCol, vN(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, decMag(iCol), uN(iCol,iVer) ', lchnk, iVer, iCol, decMag(iCol), uN(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, wN(iCol,iVer) ', lchnk, iVer, iCol, wN(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, dA(iCol,iVer) ', lchnk, iVer, iCol, dA(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, delTpDelZ, rmassOP, tP(iCol,iVer) ', lchnk, iVer, iCol, delTpDelZ, rmassOP, tP(iCol,iVer) ! if (iVer > 1 .and. c1(iCol,iVer-1) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, c1(iCol,iVer) ', lchnk, iVer, iCol, c1(iCol,iVer) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, eHPB(iCol), bMagT, min/max dipMag(iCol,1:oPBot) ', lchnk, iVer, iCol, eHPB(iCol), bMagT, MINVAL(dipMag(iCol,1:oPBot)), MAXVAL(dipMag(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, min/max vN(iCol,1:oPBot) ', lchnk, iVer, iCol, MINVAL(vN(iCol,1:oPBot)), MAXVAL(vN(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, decMag(iCol), min/max uN(iCol,1:oPBot) ', lchnk, iVer, iCol, decMag(iCol), MINVAL(uN(iCol,1:oPBot)), MAXVAL(uN(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, min/max wN(iCol,1:oPBot) ', lchnk, iVer, iCol, MINVAL(wN(iCol,1:oPBot)), MAXVAL(wN(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, min/max dA(iCol,1:oPBot) ', lchnk, iVer, iCol, MINVAL(dA(iCol,1:oPBot)), MAXVAL(dA(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, delTpDelZ, rmassOP, min/max tP(iCol,1:oPBot) ', lchnk, iVer, iCol, delTpDelZ, rmassOP, MINVAL(tP(iCol,1:oPBot)), MAXVAL(tP(iCol,1:oPBot)) ! if (MAXVAL(c1(iCol,1:oPBot)) > 1.E+30_r8) write(iulog,*)'lchnk, iVer, iCol, min/max c1(iCol,1:oPBot) ', lchnk, iVer, iCol, MINVAL(c1(iCol,1:oPBot)), MAXVAL(c1(iCol,1:oPBot)) enddo enddo write(iulog,*)'lchnk, After vertical loop optransvert ' write(iulog,*)'eHPB(1), bMagT, cos(decMag(1) ', eHPB(iCol), bMagT , cos(decMag(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,*)'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 !----------------------------------------------------------------- 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 !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 !S c41 = -dA(iCol,1) / 2._r8 / tP(iCol,1) * sin(dipMag(iCol,1))**2 * ( 2._r8 * delTpDelZ1 + rmassOp / avogad * gravit / kboltz * .01_r8 ) c41 = -dA(iCol,1) / 2._r8 / tP(iCol,1) * sin(dipMag(iCol,1))**2 * ( 2._r8 * delTpDelZ1 + rmassOp / avogad * gravity / kboltz * .01_r8 ) !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 = 1._r8 - 2._r8 * delZ(iCol,1) / c31 * c41 !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) 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) !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) rHS(1) = rHS(1) - subDiag(1) * c5 !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) !S if (newDensOp(1) < 0._r8) write(iulog,*)'lchnk, iCol, lat, lon, newDensOp(1) ', lchnk, iCol, get_lat_p(lchnk,iCol),get_lon_p(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)) !S newMMROp(1:oPBot) = newDensOp(1:oPBot) / mbarv(iCol,1:oPBot,lchnk) * rmassOp / pMid(iCol,1:oPBot) * (kboltz * tN(iCol,1:oPBot)) / 1.E-06_r8 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 before update ptend%q(iCol,1:oPBot,indxOp) ', lchnk, MINVAL(ptend%q(iCol,1:oPBot,indxOp)), MAXVAL(ptend%q(iCol,1:oPBot,indxOp)) !S if (sindxOp > 0) then !S mmrPOp(iCol,1:oPBot) = newMMROp(1:oPBot) !S else !S ptend%q(iCol,1:oPBot,indxOp) = ptend%q(iCol,1:oPBot,indxOp) + (newMMROp(1:oPBot) - state%q(iCol,1:oPBot,indxOp)) / ztodt !S endif ! write(iulog,*)'lchnk, min/max after update ptend%q(iCol,1:oPBot,indxOp) ', lchnk, MINVAL(ptend%q(iCol,1:oPBot,indxOp)), MAXVAL(ptend%q(iCol,1:oPBot,indxOp)) ! 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