subroutine chem
use modglobal, only : lchem,k1,JNO2,dt,rk3step,ib,ie,ihc,ih,jb,je,jhc,jh,kb,ke,khc,kh
use modfields, only : svp,svm,sv0,IIc
implicit none
real, dimension(ib-ihc:ie+ihc,jb-jhc:je+jhc,kb:ke+khc) :: dummyNO, dummyNO2, dummyO3
if (lchem .eqv. .false.) return
! forward Euler
! [NO]
! svp(:,:,:,1) = svp(:,:,:,1) + 30.006 * ( JNO2 * (svm(:,:,:,2)/46.005) - k1 * (svm(:,:,:,1) / 30.006) * (svm(:,:,:,3)/47.997) )
! [NO2]
! svp(:,:,:,2) = svp(:,:,:,2) + 46.005 * (-JNO2 * (svm(:,:,:,2)/46.005) + k1 * (svm(:,:,:,1) / 30.006) * (svm(:,:,:,3)/47.997) )
! [O3]
! svp(:,:,:,3) = svp(:,:,:,3) + 47.997 * ( JNO2 * (svm(:,:,:,2)/46.005) - k1 * (svm(:,:,:,1) / 30.006) * (svm(:,:,:,3)/47.997) )
if (.not. rk3step==3) return
! convert into mol/m^3
dummyNO = IIc*sv0(:,:,kb:ke+khc,1)/30.006
dummyNO2= IIc*sv0(:,:,kb:ke+khc,2)/46.005
dummyO3 = IIc*sv0(:,:,kb:ke+khc,3)/47.997
!backward Euler, semi-implicit
! sv0(:,:,kb:ke+khc,1) = 30.006 * ( ( (sv0(:,:,kb:ke+khc,1)/30.006) + JNO2 * dummyNO2 * dt) / (1. + k1 * dummyO3 * dt) )
! sv0(:,:,kb:ke+khc,2) = 46.005 * ( ( (sv0(:,:,kb:ke+khc,2)/46.005) + k1 * dummyNO * dummyO3 * dt ) / (1. + JNO2 * dt) )
! sv0(:,:,kb:ke+khc,3) = 47.997 * ( ( (sv0(:,:,kb:ke+khc,3)/47.997) + JNO2 * dummyNO2 * dt) / (1. + k1 * dummyNO * dt) )
!backward Euler, fully implicit. Derivation at (/projects/Chemistry/FullyImplicit.mw)
sv0(:,:,kb:ke+khc,1) = 30.006 * ( (sv0(:,:,kb:ke+khc,1)/30.006) + ( dt * (- k1 * dummyNO * dummyO3 + JNO2 * dummyNO2) ) / &
( 1. + ( ( dummyNO + dummyO3 ) * k1 + JNO2 ) * dt ) )
sv0(:,:,kb:ke+khc,2) = 46.005 * ( (sv0(:,:,kb:ke+khc,2)/46.005) - ( dt * (- k1 * dummyNO * dummyO3 + JNO2 * dummyNO2) ) / &
( 1. + ( ( dummyNO + dummyO3 ) * k1 + JNO2 ) * dt ) )
sv0(:,:,kb:ke+khc,3) = 47.997 * ( (sv0(:,:,kb:ke+khc,3)/47.997) + ( dt * (- k1 * dummyNO * dummyO3 + JNO2 * dummyNO2) ) / &
( 1. + ( ( dummyNO + dummyO3 ) * k1 + JNO2 ) * dt ) )
! alternative method in Zhong 2017 eqn. 7, analytical solution!
end subroutine chem