Line data    Source code

       1             :       module mo_msis_ubc
       2             : !---------------------------------------------------------------
       3             : !       ... msis upper bndy values
       4             : !---------------------------------------------------------------
       5             : 
       6             :       use shr_kind_mod,     only: r8 => shr_kind_r8
       7             :       use constituents,     only: pcnst
       8             : 
       9             :       use cam_abortutils,   only: endrun
      10             :       use cam_logfile,      only: iulog
      11             :       use cam_history,      only: addfld, horiz_only, outfld
      12             : 
      13             :       implicit none
      14             : 
      15             :       private
      16             :       public  :: msis_ubc_inti, get_msis_ubc, msis_timestep_init
      17             : 
      18             :       save
      19             : 
      20             :       integer                :: ndx_n=-1, ndx_h=-1, ndx_o=-1, ndx_o2=-1 ! n, h, o, o2 spc indicies
      21             :       real(r8), allocatable  :: msis_ubc(:,:,:)                       ! module array for msis ub values (kg/kg)
      22             : 
      23             :       logical                :: zonal_average         = .false.       ! use zonal averaged tgcm values
      24             : 
      25             :       contains
      26             : 
      27           0 :       subroutine msis_ubc_inti( zonal_avg_in, n_ndx_in,h_ndx_in,o_ndx_in,o2_ndx_in )
      28             : !------------------------------------------------------------------
      29             : !       ... initialize upper boundary values
      30             : !------------------------------------------------------------------
      31             : 
      32             :       use ppgrid, only : pcols, begchunk, endchunk
      33             : 
      34             : !------------------------------------------------------------------
      35             : !       ... dummy args
      36             : !------------------------------------------------------------------
      37             :       logical, intent(in) :: zonal_avg_in  ! zonal averaging switch
      38             :       integer, intent(in) :: n_ndx_in,h_ndx_in,o_ndx_in,o2_ndx_in
      39             : 
      40             : !------------------------------------------------------------------
      41             : !       ... local variables
      42             : !------------------------------------------------------------------
      43             :       integer  :: astat
      44             :       real(r8) :: msis_switches(25) = 1._r8
      45             : 
      46           0 :       zonal_average = zonal_avg_in
      47             : 
      48           0 :       if (h_ndx_in>0) then
      49           0 :          ndx_h = h_ndx_in
      50             :       endif
      51           0 :       if (n_ndx_in>0) then
      52           0 :          ndx_n = n_ndx_in
      53             :       endif
      54           0 :       if (o_ndx_in>0) then
      55           0 :          ndx_o = o_ndx_in
      56             :       endif
      57           0 :       if (o2_ndx_in>0) then
      58           0 :          ndx_o2 = o2_ndx_in
      59             :       endif
      60             : 
      61             : !------------------------------------------------------------------
      62             : !       ... allocate msis ubc array
      63             : !------------------------------------------------------------------
      64           0 :       allocate( msis_ubc(pcols,6,begchunk:endchunk),stat=astat )
      65           0 :       if( astat /= 0 ) then
      66           0 :          write(iulog,*) 'msis_ubc_inti: failed to allocate msis_ubc; error = ',astat
      67           0 :          call endrun('msis_ubc_inti: failed to allocate msis_ubc')
      68             :       end if
      69             : 
      70           0 :       if( zonal_average ) then
      71           0 :          msis_switches(7:8)   = 0._r8
      72           0 :          msis_switches(10:14) = 0._r8
      73             :       end if
      74             : 
      75             : !------------------------------------------------------------------
      76             : !       ... initialize msis switches
      77             : !------------------------------------------------------------------
      78           0 :       call tselec( msis_switches )
      79             : 
      80           0 :       call addfld( 'MSIS_T', horiz_only, 'A', 'K',     'T upper boundary condition from MSIS')
      81           0 :       call addfld( 'MSIS_H', horiz_only, 'A', 'kg/kg', 'H upper boundary condition from MSIS')
      82           0 :       call addfld( 'MSIS_N', horiz_only, 'A', 'kg/kg', 'N upper boundary condition from MSIS')
      83           0 :       call addfld( 'MSIS_O', horiz_only, 'A', 'kg/kg', 'O upper boundary condition from MSIS')
      84           0 :       call addfld( 'MSIS_O2',horiz_only, 'A', 'kg/kg', 'O2 upper boundary condition from MSIS')
      85             : 
      86           0 :       end subroutine msis_ubc_inti
      87             : 
      88           0 :       subroutine msis_timestep_init( ap, f107p_in, f107a_in )
      89             : !--------------------------------------------------------------------
      90             : !       ... get the upper boundary values for h, n, o, o2 and temp
      91             : !--------------------------------------------------------------------
      92             : 
      93             :       use ppgrid,       only : pcols, begchunk, endchunk
      94             :       use constituents, only : cnst_mw
      95             :       use time_manager, only : get_curr_date, get_calday, get_curr_calday
      96             :       use phys_grid,    only : get_ncols_p, get_rlon_all_p, get_rlat_all_p
      97             :       use ref_pres,     only : ptop_ref
      98             :       use spmd_utils,   only : masterproc
      99             :       use physconst,    only : pi
     100             :       use cam_control_mod,only : lambm0, eccen, mvelpp, obliqr
     101             :       use shr_orb_mod,    only : shr_orb_decl
     102             : 
     103             : !--------------------------------------------------------------------
     104             : !       ... dummy args
     105             : !--------------------------------------------------------------------
     106             :       real(r8), intent(in)    ::  ap
     107             :       real(r8), intent(in)    ::  f107p_in ! previous day
     108             :       real(r8), intent(in)    ::  f107a_in
     109             : 
     110             : !--------------------------------------------------------------------
     111             : !       ... local variables
     112             : !--------------------------------------------------------------------
     113             :       real(r8), parameter :: mass_switch = 48._r8
     114             :       real(r8), parameter :: pa2mb       = 1.e-2_r8       ! pascal to mb
     115             :       real(r8), parameter :: amu_fac     = 1.65979e-24_r8 ! g/amu
     116             :       real(r8), parameter :: r2d         = 180._r8/pi
     117             :       integer  ::  i, c, ncol
     118             :       integer  ::  yr, mon, day, tod
     119             :       integer  ::  yrday
     120             :       integer  ::  date
     121             :       real(r8) ::  alt, solar_time, ut, rtod, doy
     122             :       real(r8) ::  msis_press
     123             :       real(r8) ::  msis_ap(7)
     124             :       real(r8) ::  msis_temp(2)
     125             :       real(r8) ::  msis_conc(9)
     126             :       real(r8) ::  rlons(pcols)
     127             :       real(r8) ::  rlats(pcols)
     128             :       real(r8) ::  dnom(pcols)
     129             :       real(r8) ::  pint(pcols)       ! top interface pressure (Pa)
     130             :       real(r8) ::  calday, delta, esfact
     131             :       real(r8) ::  f107p, f107a
     132             : 
     133             :       !--------------------------------------------------------------------
     134             :       ! ... get values from msis
     135             :       !--------------------------------------------------------------------
     136             : 
     137           0 :       call get_curr_date( yr, mon, day, tod )
     138           0 :       tod = 0
     139           0 :       rtod       = tod
     140           0 :       ut         = rtod/3600._r8
     141           0 :       date       = 10000*yr + 100*mon + day
     142           0 :       doy        = get_calday( date, tod )
     143           0 :       msis_ap(:) = 0._r8
     144           0 :       msis_ap(1) = ap
     145           0 :       pint(:)    = ptop_ref
     146             : 
     147           0 :       calday = get_curr_calday()
     148             : 
     149           0 :       esfact = 1._r8
     150           0 :       call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr, delta, esfact )
     151             : 
     152           0 :       f107p = esfact*f107p_in
     153           0 :       f107a = esfact*f107a_in
     154             : 
     155             : #ifdef MSIS_DIAGS
     156             :       if( masterproc ) then
     157             :          write(iulog,*) '===================================='
     158             :          write(iulog,*) 'msis_timestep_init: diagnostics'
     159             :          write(iulog,*) 'yr,mon,day,tod,date,ut,doy,esfact = ', yr, mon, day, tod, date, ut, doy, esfact
     160             :          write(iulog,*) '===================================='
     161             :       end if
     162             : #endif
     163           0 :       chunk_loop : do c = begchunk,endchunk
     164           0 :          ncol = get_ncols_p( c )
     165           0 :          call get_rlat_all_p( c, ncol, rlats )
     166           0 :          call get_rlon_all_p( c, ncol, rlons )
     167           0 :          rlons(:ncol) = r2d * rlons(:ncol)
     168           0 :          rlats(:ncol) = r2d * rlats(:ncol)
     169           0 :          yrday = mod( yr,100 ) * 1000 + int( doy )
     170           0 :          column_loop : do i = 1,ncol
     171           0 :             solar_time = ut + rlons(i)/15._r8
     172           0 :             msis_press = pint(i)*pa2mb
     173             :             call ghp7( yrday, rtod, alt, rlats(i), rlons(i), &
     174             :                  solar_time, f107a, f107p, msis_ap, msis_conc, &
     175             :                  msis_temp, msis_press )
     176           0 :             msis_ubc(i,1,c) = msis_temp(2)              ! temp (K)
     177             : #ifdef MSIS_DIAGS
     178             :             write(iulog,*) '===================================='
     179             :             write(iulog,*) 'msis_timestep_init: diagnostics for col,chnk = ',i,c
     180             :             write(iulog,*) 'yrday, rtod, alt,press = ',yrday,rtod,alt,msis_press
     181             :             write(iulog,*) 'msis_temp = ',msis_temp(2)
     182             : #endif
     183           0 :             msis_ubc(i,2,c) = msis_conc(7)           ! h (molec/cm^3)
     184           0 :             msis_ubc(i,3,c) = msis_conc(8)           ! n (molec/cm^3)
     185           0 :             msis_ubc(i,4,c) = msis_conc(2)           ! o (molec/cm^3)
     186           0 :             msis_ubc(i,5,c) = msis_conc(4)           ! o2 (molec/cm^3)
     187           0 :             msis_ubc(i,6,c) = msis_conc(6)           ! total atm dens (g/cm^3)
     188             : 
     189             : #ifdef MSIS_DIAGS
     190             :             write(iulog,*) 'msis h,n,o,o2,m = ',msis_ubc(i,2:6,c)
     191             :             write(iulog,*) '===================================='
     192             : #endif
     193             :          end do column_loop
     194             : 
     195             :          !--------------------------------------------------------------------
     196             :          !      ... transform from molecular density to mass mixing ratio
     197             :          !--------------------------------------------------------------------
     198           0 :          dnom(:ncol) = amu_fac/msis_ubc(:ncol,6,c)
     199           0 :          if( ndx_h > 0 ) then
     200           0 :             msis_ubc(:ncol,2,c) = cnst_mw(ndx_h)*msis_ubc(:ncol,2,c)*dnom(:ncol)
     201             :          end if
     202           0 :          if( ndx_n > 0 ) then
     203           0 :             msis_ubc(:ncol,3,c) = cnst_mw(ndx_n)*msis_ubc(:ncol,3,c)*dnom(:ncol)
     204             :          end if
     205           0 :          if( ndx_o > 0 ) then
     206           0 :             msis_ubc(:ncol,4,c) = cnst_mw(ndx_o)*msis_ubc(:ncol,4,c)*dnom(:ncol)
     207             :          end if
     208           0 :          if( ndx_o2 > 0 ) then
     209           0 :             msis_ubc(:ncol,5,c) = cnst_mw(ndx_o2)*msis_ubc(:ncol,5,c)*dnom(:ncol)
     210             :          end if
     211             :       end do chunk_loop
     212             : 
     213           0 :       end subroutine msis_timestep_init
     214             : 
     215           0 :       subroutine get_msis_ubc( lchunk, ncol, temp, mmr )
     216             : !--------------------------------------------------------------------
     217             : !       ... get the upper boundary values for h, n, o, o2 and temp
     218             : !--------------------------------------------------------------------
     219             : 
     220           0 :       use ppgrid,       only : pcols
     221             : 
     222             : !--------------------------------------------------------------------
     223             : !       ... dummy args
     224             : !--------------------------------------------------------------------
     225             :       integer, intent(in)     :: lchunk            ! chunk id
     226             :       integer, intent(in)     :: ncol              ! columns in chunk
     227             :       real(r8), intent(out)   :: temp(pcols)       ! msis temperature at top interface (K)
     228             :       real(r8), intent(inout) :: mmr(pcols,pcnst)  ! msis concentrations at top interface (kg/kg)
     229             : 
     230             : !--------------------------------------------------------------------
     231             : !       ... set model ubc values from msis
     232             : !--------------------------------------------------------------------
     233           0 :       temp(:ncol) = msis_ubc(:ncol,1,lchunk)
     234             : 
     235           0 :       call outfld( 'MSIS_T', msis_ubc(:ncol,1,lchunk), ncol, lchunk)
     236           0 :       call outfld( 'MSIS_H', msis_ubc(:ncol,2,lchunk), ncol, lchunk)
     237           0 :       call outfld( 'MSIS_N', msis_ubc(:ncol,3,lchunk), ncol, lchunk)
     238           0 :       call outfld( 'MSIS_O', msis_ubc(:ncol,4,lchunk), ncol, lchunk)
     239           0 :       call outfld( 'MSIS_O2',msis_ubc(:ncol,5,lchunk), ncol, lchunk)
     240             : 
     241           0 :       if( ndx_h > 0 ) then
     242           0 :          mmr(:ncol,ndx_h) = msis_ubc(:ncol,2,lchunk)
     243             :       end if
     244           0 :       if( ndx_n > 0 ) then
     245           0 :          mmr(:ncol,ndx_n) = msis_ubc(:ncol,3,lchunk)
     246             :       end if
     247           0 :       if( ndx_o > 0 ) then
     248           0 :          mmr(:ncol,ndx_o) = msis_ubc(:ncol,4,lchunk)
     249             :       end if
     250           0 :       if( ndx_o2 > 0 ) then
     251           0 :          mmr(:ncol,ndx_o2) = msis_ubc(:ncol,5,lchunk)
     252             :       end if
     253             : 
     254           0 :       end subroutine get_msis_ubc
     255             : 
     256             :       end module mo_msis_ubc