!
      module msis_module
!
! calculate background temperature from NRL-MSIS00
!
      use params_module,only: nlon,nlonp2,nlonp4,nlat
      
      implicit none
      real,allocatable,dimension(:,:) :: msis_tn,   ! background temperature from MSIS
     |    msis_z                      ! background geopotential height for pressure level
      logical, parameter :: msis = .true.
      real, parameter :: ap_msis = 4. ! from M. Hagan
      integer, parameter :: mass = 0  ! only temperature from MSIS
!
      logical, parameter :: debug = .false.
!
      contains
!-----------------------------------------------------------------------
      subroutine msis_bndry(z,lon0,lon1,lev0,lev1,lat0,lat1)
!
! driver routione for NRL-MSIS00
! calculate the z at all the longitudes to get the longitudinal averaged
! mean (with annual tides), and if no annual tides are used, then calculate 
! the z at the full domain to get the global mean
! neutral temperature is not averaged
!
      use input_module,only: f107,f107a,tideann
      use init_module, only: iyear,iday,uthr
      use cons_module,only: ylatg,ylong,rtd, ! geographic latitude/longitudes (radians)
     |     p0,   !  5.0e-4 microbars standard pressure
     |     r0 ,  !  re+h0 with h0 = 80 km [cm]
     |     dphi,  ! delta lat (pi/nlat)  
     |     dlamda ! delta lon (2pi/nlon) 
      use params_module,only: zibot ! bottom interface level
!
! Routine to add fields to secondary histories:
      use addfld_module,only: addfld
      
#ifdef MPI
      use mpi_module,only: mp_periodic_f2d
#endif

! Args:
      integer,intent(in) :: lon0,lon1,lev0,lev1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(in) :: z   
!
!      
      integer :: i,j,iyd,icentury,iutsec,lonbeg,lonend,
     |   latbeg,latend
      real :: alt,glon,glat,slt,press,sum_z,sum_glm,sum,fac
      real :: z_loc(1:nlonp4,1:nlat), tn_loc(1:nlonp4,1:nlat)
      real :: d(8),tt(2) ! msis output tt(2) temperature
      real :: sw(25)      
      data (sw(i),i=1,25)/25*1./
!      
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = lon1-2
      
!     if(tideann == 0) then   ! calculate global mean
!       latbeg = 1
!       latend = nlat
!     else
!       latbeg = lat0
!       latend = lat1     
!     endif
      latbeg = 1
      latend = nlat
!      
      icentury = int(float(iyear)/100.)
      iyd = (iyear- icentury*100)*1000+iday  ! YYDDD 
      iutsec = int(uthr*3600)      
! 
! pressure level in mB for MSIS
      press = p0*exp(-1.*zibot)*1.e-3   ! 1.e-3 to convert from micro- to millibar
!      
! set switches for zonal mean (same as SR setatmos in GSWM)
      sw(7)=0.		! diurnal
      sw(8)=0.		! semidiurnal
      sw(14)=0. 	! terdiurnal
      sw(10)=1. 	! all ut/long. effects
      call tselec(sw)	! turn on/off particular variations
!     
      if(debug) write(6,*) 'msis: ',iyd,uthr,f107a,f107,ap_msis
!      
      sum_glm = 0.               ! global mean
      sum     = 0.
      fac     = r0*r0*dphi*dlamda
!      
      do j = latbeg,latend
	glat = ylatg(j)*rtd      ! geo.latitude rad to deg
	sum_z = 0.               ! zonal mean
        do i = 3,nlonp2          ! over all longitude to get zonal mean
	  glon = ylong(i-2)*rtd  ! geo.longitude rad to deg
	  slt  = uthr + glon/15. ! apparent solar local time
!	  
! MSIS for certain height alt	
! output background TN
!          alt = z(1,i,j)*1.e-5 ! z-full level  [cm->km]  
!          call gtd7(iyd,iutsec,alt,glat,glon,slt,f107a,f107,ap_msis,
!     |        mass,d,tt)
     	  
! MSIS for pressure level press at lower boundary
! output background TN + Z
          call ghp7(iyd,iutsec,alt,glat,glon,slt,f107a,f107,ap_msis,
     |        d,tt,press)
!     
          sum_z      = sum_z + alt*1.e5
!         sum        = sum + fac*cos(ylatg(j))               ! for global mean
!         sum_glm    = sum_glm + alt*1.e5*fac*cos(ylatg(j))  ! for global mean
          tn_loc(i,j)= tt(2)
	  
          if(debug) write(6,*)'msis_loc ',glon,glat,alt,slt,msis_tn(i,j)
        enddo
	
	sum_z = sum_z/nlon  ! longitudinal averaged value
 	z_loc(:,j) = sum_z   
	
      enddo
!
! copied local field into subdomain msis fields      
!     if(tideann == 0) then
!       fac = sum_glm/sum
!       msis_z(lonbeg:lonend,lat0:lat1) = fac
!       z_loc(lonbeg:lonend,lat0:lat1)  = fac 
!     else
!       msis_z(lonbeg:lonend,lat0:lat1) = z_loc(lonbeg:lonend,lat0:lat1)
!     endif
!
! 4/20/06 btf: use zonal mean msis Z:
      msis_z(lonbeg:lonend,lat0:lat1) = z_loc(lonbeg:lonend,lat0:lat1)
!
      msis_tn(lonbeg:lonend,lat0:lat1)= tn_loc(lonbeg:lonend,lat0:lat1)
	
! periodic points (no mpi since we have local global field)
      if(lon0 == 1) then
        do j=lat0,lat1
          do i=1,2
            msis_tn(i,j) = tn_loc(nlon+i,j)
            msis_z(i,j)  = z_loc(nlon+i,j)
          enddo
        enddo
      endif
      
      if(lon1 == nlonp4) then
        do j=lat0,lat1
          do i=1,2
            msis_tn(nlonp2+i,j)= tn_loc(i+2,j)
            msis_z(nlonp2+i,j) = z_loc(i+2,j)
          enddo
        enddo
      endif 
!      
      call addfld('MSIS_TN','TN back(msis)','K',
     |    msis_tn(lon0:lon1,lat0:lat1),'lon',lon0,lon1,'lat',
     |    lat0,lat1,0)
      call addfld('MSIS_Z','Z(msis)','cm',msis_z(lon0:lon1,lat0:lat1),
     |    'lon',lon0,lon1,'lat',lat0,lat1,0)
!      
 99   return
      end subroutine msis_bndry
!-----------------------------------------------------------------------
      subroutine alloc_msis(lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
      allocate(msis_tn(lon0-2:lon1+2,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_msis: error allocating',
     |  ' msis_tn: stat=',i3)") istat
     
      allocate(msis_z(lon0-2:lon1+2,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_msis: error allocating',
     |  ' msis_z: stat=',i3)") istat
     
      end subroutine alloc_msis
!-----------------------------------------------------------------------

      end module msis_module