LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_apex.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 11 | 140 | 7.9 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 4 | 25.0 %

              Line data    Source code
           1             : module mo_apex
           2             :
           3             : !-------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !
           6             : !   Calculate apex coordinates and magnetic field magnitudes
           7             : !   at global geographic grid for year of current model run.
           8             : !
           9             : ! Method:
          10             : !
          11             : !   The magnetic field parameters output by this module are time and height
          12             : !     independent. They are chunked for waccm physics, i.e., allocated as
          13             : !     (pcols,begchunk:endchunk)
          14             : !   Interface sub apexmag is called once per run from sub inti.
          15             : !     Sub apexmag may be called for years 1900 through 2005.
          16             : !   This module is dependent on routines in apex_subs.F (modified IGRF model).
          17             : !   Apex_subs has several authors, but has been modified and maintained
          18             : !     in recent years by Roy Barnes (bozo@ucar.edu).
          19             : !   Subs apxmka and apxmall are called with the current lat x lon grid
          20             : !     resolution.
          21             : !
          22             : ! Author: Ben Foster, foster@ucar.edu (Nov, 2003)
          23             : !-------------------------------------------------------------------------------
          24             :
          25             :    use shr_kind_mod,    only: r8 => shr_kind_r8
          26             :    use ppgrid,          only: pcols, begchunk, endchunk          ! physics grid
          27             :    use cam_abortutils,  only: endrun
          28             :    use cam_logfile,     only: iulog
          29             :    use spmd_utils,      only: masterproc
          30             :    use apex,            only: apex_mka, apex_mall, apex_dypol, apex_set_igrf
          31             :    use apex,            only: apex_beg_yr, apex_end_yr
          32             :    implicit none
          33             :
          34             :    private
          35             :    public :: mo_apex_readnl
          36             :    public :: mo_apex_init
          37             :    public :: mo_apex_init1
          38             :    public :: alatm, alonm, bnorth, beast, bdown, bmag
          39             :    public :: d1vec, d2vec, colatp, elonp
          40             :    public :: maglon0 ! geographic longitude at the equator where geomagnetic longitude is zero (radians)
          41             :
          42             :    ! year to initialize apex
          43             :    real(r8), public, protected :: geomag_year = -1._r8
          44             :    logical, public, protected :: geomag_year_updated = .true.
          45             :
          46             :    integer :: fixed_geomag_year = -1
          47             :
          48             : !-------------------------------------------------------------------------------
          49             : ! Magnetic field output arrays, chunked for physics:
          50             : ! (these are allocated (pcols,begchunk:endchunk) by sub allocate_arrays)
          51             : !-------------------------------------------------------------------------------
          52             :    real(r8), protected, allocatable, dimension(:,:) :: & ! (pcols,begchunk:endchunk)
          53             :      alatm,  & ! apex mag latitude at each geographic grid point (radians)
          54             :      alonm,  & ! apex mag longitude at each geographic grid point (radians)
          55             :      bnorth, & ! northward component of magnetic field
          56             :      beast,  & ! eastward component of magnetic field
          57             :      bdown,  & ! downward component of magnetic field
          58             :      bmag      ! magnitude of magnetic field
          59             :    real(r8), protected, allocatable, dimension(:,:,:) :: & ! (3,pcols,begchunk:endchunk)
          60             :      d1vec,    & ! base vectors more-or-less magnetic eastward direction
          61             :      d2vec       ! base vectors more-or-less magnetic downward/equatorward direction
          62             :    real(r8), protected :: &
          63             :      colatp,   & ! geocentric colatitude of geomagnetic dipole north pole (deg)
          64             :      elonp       ! East longitude of geomagnetic dipole north pole (deg)
          65             :
          66             :    real(r8), protected :: maglon0
          67             :
          68             :    character(len=256) :: igrf_geomag_coefs_file = 'igrf_geomag_coefs_file'
          69             :
          70             : contains
          71             :
          72             : !======================================================================
          73             : !======================================================================
          74        1536 : subroutine mo_apex_readnl(nlfile)
          75             :
          76             :   use namelist_utils, only : find_group_name
          77             :   use units,          only : getunit, freeunit
          78             :   use spmd_utils,     only : mpicom, masterprocid, mpi_integer, mpi_character
          79             :
          80             :   character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          81             :
          82             :   ! Local variables
          83             :   integer :: unitn, ierr
          84             :   character(len=*), parameter :: subname = 'mo_apex_readnl'
          85             :
          86             :   namelist /geomag_nl/ fixed_geomag_year, igrf_geomag_coefs_file
          87             :
          88             :   ! Read namelist
          89        1536 :   if (masterproc) then
          90           2 :      unitn = getunit()
          91           2 :      open( unitn, file=trim(nlfile), status='old' )
          92           2 :      call find_group_name(unitn, 'geomag_nl', status=ierr)
          93           2 :      if (ierr == 0) then
          94           0 :         read(unitn, geomag_nl, iostat=ierr)
          95           0 :         if (ierr /= 0) then
          96           0 :            call endrun(subname // ':: ERROR reading namelist')
          97             :         end if
          98             :      end if
          99           2 :      close(unitn)
         100           2 :      call freeunit(unitn)
         101             :   end if
         102             :
         103             :   ! Broadcast namelist variables
         104        1536 :   call mpi_bcast(fixed_geomag_year,  1, mpi_integer, masterprocid, mpicom, ierr)
         105        1536 :   call mpi_bcast(igrf_geomag_coefs_file,  len(igrf_geomag_coefs_file), mpi_character, masterprocid, mpicom, ierr)
         106             :
         107        1536 : end subroutine mo_apex_readnl
         108             :
         109             : !======================================================================
         110             : !======================================================================
         111           0 : subroutine mo_apex_init1()
         112             :    use time_manager,  only: get_curr_date
         113             :    use phys_grid,     only: get_grid_dims
         114             :
         115             :    integer  :: i, j, ist          ! indices
         116             :
         117             :    integer :: nglats
         118             :    integer :: nglons
         119             :    integer, parameter :: ngalts = 2             ! number of altitudes
         120             :
         121           0 :    real(r8), allocatable :: gridlats(:)
         122           0 :    real(r8), allocatable :: gridlons(:)
         123             :    real(r8) :: gridalts(ngalts)                   ! altitudes passed to apxmka
         124             :
         125             :    integer :: ngcols, hdim1_d, hdim2_d
         126             :    integer :: yr, mon, day, sec
         127             :
         128             :    ! read the IGRF coefs from file
         129           0 :    call apex_set_igrf( igrf_geomag_coefs_file )
         130             :
         131           0 :    if (fixed_geomag_year>0) then
         132           0 :       yr = fixed_geomag_year
         133             :    else
         134           0 :       call get_curr_date(yr, mon, day, sec)
         135             :    end if
         136             :
         137           0 :    if ( yr < apex_beg_yr )   yr = apex_beg_yr
         138           0 :    if ( yr > apex_end_yr-1 ) yr = apex_end_yr-1
         139             :
         140           0 :    if (.not.(yr > geomag_year)) then
         141           0 :       geomag_year_updated = .false.
         142             :       return
         143             :    else
         144           0 :       geomag_year_updated = .true.
         145             :    endif
         146             :
         147           0 :    geomag_year = dble(yr)+0.5_r8
         148             :
         149             : !-------------------------------------------------------------------------------
         150             : ! Center min, max altitudes about 130 km
         151             : !-------------------------------------------------------------------------------
         152           0 :    gridalts(:ngalts) =  (/ 90._r8, 170._r8 /)
         153             :
         154             : !-------------------------------------------------------------------------------
         155             : ! Initialize APEX with a regular lat/lon grid ...
         156             : ! (Note apex_mka expects longitudes in -180 -> +180)
         157             : !-------------------------------------------------------------------------------
         158           0 :    call get_grid_dims(hdim1_d,hdim2_d)
         159           0 :    ngcols = hdim1_d*hdim2_d
         160           0 :    if (     ngcols < 1000 ) then ! 10-degrees
         161           0 :       nglats = 19
         162           0 :       nglons = 37
         163           0 :    elseif ( ngcols < 10000 ) then ! 5-degrees
         164           0 :       nglats = 37
         165           0 :       nglons = 73
         166           0 :    elseif ( ngcols < 20000 ) then ! 2-degree
         167           0 :       nglats = 91
         168           0 :       nglons = 181
         169           0 :    elseif ( ngcols < 100000 ) then ! 1-degree
         170           0 :       nglats = 181
         171           0 :       nglons = 361
         172             :    else                            ! half-degee
         173           0 :       nglats = 361
         174           0 :       nglons = 721
         175             :    endif
         176             :
         177           0 :    allocate ( gridlats(nglats), gridlons(nglons) )
         178           0 :    do i = 1,nglons
         179           0 :       gridlons(i) = -180._r8 + dble(i-1)*360._r8/(nglons-1)
         180             :    enddo
         181           0 :    do j = 1,nglats
         182           0 :       gridlats(j) = -90._r8 + dble(j-1)*180._r8/(nglats-1)
         183             :    enddo
         184             :
         185             :    call apex_mka( geomag_year, gridlats, gridlons, gridalts, &
         186           0 :                   nglats, nglons, ngalts, ist )
         187             :
         188           0 :    if( ist /= 0 ) then
         189           0 :       write(iulog,"(/,'>>> mo_apex_init: Error from apxmka: ist=',i5)") ist
         190           0 :       call endrun("mo_apex_init: Error from apxmka")
         191             :    end if
         192             :
         193           0 :    deallocate( gridlats, gridlons )
         194             :
         195           0 :    if (masterproc) then
         196           0 :       if (fixed_geomag_year<1) then
         197           0 :          write(iulog, "('mo_apex_init: model yr,mon,day,sec ',4i6)") yr, mon, day, sec
         198             :       endif
         199           0 :       write(iulog, "('mo_apex_init: nglons,nglats ', 2i6)") nglons, nglats
         200             :    endif
         201             :
         202           0 : end subroutine mo_apex_init1
         203             :
         204             : !======================================================================
         205             : !======================================================================
         206           0 : subroutine mo_apex_init(phys_state)
         207             : !-------------------------------------------------------------------------------
         208             : ! Driver for apex code to calculate apex magnetic coordinates at
         209             : !   current geographic spatial resolution for given year. This calls
         210             : !   routines in apex_subs.F.
         211             : !
         212             : ! This is called once per run from sub inti.
         213             : !-------------------------------------------------------------------------------
         214             :
         215           0 :    use physconst,only : pi
         216             :    use physics_types, only: physics_state
         217             :    use epp_ionization,only: epp_ionization_setmag
         218             :
         219             :    ! Input/output arguments
         220             :    type(physics_state), intent(in), dimension(begchunk:endchunk) :: phys_state
         221             :
         222             : !-------------------------------------------------------------------------------
         223             : ! Local variables
         224             : !-------------------------------------------------------------------------------
         225             :    real(r8), parameter :: re    = 6.378165e8_r8 ! earth radius (cm)
         226             :    real(r8), parameter :: h0    = 9.0e6_r8      ! base height (90 km)
         227             :    real(r8), parameter :: hs    = 1.3e7_r8
         228             :    real(r8), parameter :: eps   = 1.e-6_r8      ! epsilon
         229             :    real(r8), parameter :: cm2km = 1.e-5_r8
         230             :
         231             :    integer  :: c, i, ist           ! indices
         232             :    integer  :: ncol
         233             :
         234             :    real(r8) :: alt, hr, alon, alat, & ! apxmall args
         235             :                vmp, w, d, be3, sim, xlatqd, f, si, collat, collon
         236             :
         237             : !-------------------------------------------------------------------------------
         238             : ! Non-scalar arguments returned by APXMALL:
         239             : !-------------------------------------------------------------------------------
         240             :    real(r8) :: bhat(3)
         241             :    real(r8) :: d3(3)
         242             :    real(r8) :: e1(3), e2(3), e3(3)
         243             :    real(r8) :: f1(2), f2(2)
         244             :
         245             :    real(r8) :: bg(3), d1g(3), d2g(3), bmg
         246             :
         247             :    real(r8) :: rdum
         248             :
         249           0 :    real(r8) :: maglat(pcols,begchunk:endchunk)
         250             :
         251             :    real(r8), parameter :: rtd = 180._r8/pi     ! radians to degrees
         252             :    real(r8), parameter :: dtr = pi/180._r8     ! degrees to radians
         253             :
         254           0 :    call mo_apex_init1()
         255           0 :    if ((.not.geomag_year_updated) .and. (allocated(alatm))) return
         256             :
         257             : !-------------------------------------------------------------------------------
         258             : ! Allocate output arrays
         259             : !-------------------------------------------------------------------------------
         260           0 :    call allocate_arrays()
         261             :
         262           0 :    alt   = hs*cm2km    ! altitude for apxmall (km)
         263           0 :    hr    = alt         ! reference altitude (km)
         264             :
         265             : !------------------------------------------------------------------------------
         266             : ! Apex coords alon, alat are returned for each geographic grid point:
         267             : ! first form global arrays
         268             : !------------------------------------------------------------------------------
         269           0 :    do c = begchunk, endchunk
         270           0 :       ncol = phys_state(c)%ncol
         271           0 :       do i = 1,ncol
         272           0 :          collat = phys_state(c)%lat(i)*rtd  ! latitude of current column (deg)
         273           0 :          collon = phys_state(c)%lon(i)*rtd  ! latitude of current column (deg)
         274           0 :          if ( collon < -180._r8 ) collon = collon+360._r8
         275           0 :          if ( collon >  180._r8 ) collon = collon-360._r8
         276             :          call apex_mall(                                  &
         277             :            collat, collon, alt, hr,           & ! Inputs
         278           0 :            bg, bhat, bmag(i,c), si,     & ! Mag Fld
         279             :            alon, alat,                                  & ! Apex lon,lat output
         280             :            vmp, w, d, be3, sim, d1vec(:,i,c), d2vec(:,i,c),  d3, e1, e2, e3, & ! Mod Apex
         281           0 :            xlatqd, f, f1, f2, ist )                       ! Qsi-Dpl
         282           0 :          if( ist /= 0 ) then
         283           0 :            write(iulog,"(/,'>>> mo_apex_init: Error from apxmall: ist=',i4)") ist
         284           0 :            call endrun('mo_apex_init: Error from apxmall')
         285             :          end if
         286           0 :          beast (i,c) = bg(1)
         287           0 :          bnorth(i,c) = bg(2)
         288           0 :          bdown (i,c) = -bg(3)
         289           0 :          alonm (i,c) = alon*dtr       ! mag lons (radians)
         290           0 :          alatm (i,c) = alat*dtr       ! mag lats (radians)
         291           0 :          maglat(i,c) = alat  ! mag lats (degrees)
         292             :       enddo
         293             :    enddo
         294             :
         295             :    ! find geograghic latitude ( maglon0 ) where the geomagnetic latitude is zero at the equator
         296             :    ! by first extracting the geographic coordinates at zero degrees longitude ...
         297           0 :    collat = 0._r8
         298           0 :    collon = 0._r8
         299             :    call apex_mall(                                     &
         300             :         collat, collon, alt, hr,                       & ! Inputs
         301             :         bg, bhat, bmg, si,                             & ! Mag Fld
         302             :         alon, alat,                                    & ! Apex lon,lat output
         303             :         vmp, w, d, be3, sim, d1g, d2g, d3, e1, e2, e3, & ! Mod Apex
         304           0 :         xlatqd, f, f1, f2, ist )                         ! Qsi-Dpl
         305             :
         306           0 :    if( ist /= 0 ) then
         307           0 :       write(iulog,"(/,'>>> mo_apex_init: Error from apxmall: ist=',i4)") ist
         308           0 :       call endrun('mo_apex_init: Error from apxmall')
         309             :    end if
         310             :
         311           0 :    maglon0 = -alon*dtr ! (radians) geograghic latitude where the geomagnetic latitude is zero
         312             :                        ! where longitude ranges from -180E to 180E
         313             :
         314           0 :    call apex_dypol( colatp, elonp, rdum )       ! get geomagnetic dipole north pole
         315             :
         316           0 :    if (masterproc) then
         317           0 :       write(iulog, "('mo_apex_init: colatp,elonp ', 2f12.6)") colatp, elonp
         318           0 :       write(iulog, "('mo_apex_init: Calculated apex magnetic coordinates for year AD ',f8.2)") geomag_year
         319             :    endif
         320             :
         321           0 :    call epp_ionization_setmag(maglat)
         322             :
         323           0 : end subroutine mo_apex_init
         324             :
         325           0 : subroutine allocate_arrays
         326             : !------------------------------------------------------------------------------
         327             : ! Allocate module output arrays for chunked physics grid.
         328             : !------------------------------------------------------------------------------
         329             :
         330             : !------------------------------------------------------------------------------
         331             : ! local variables
         332             : !------------------------------------------------------------------------------
         333             :   integer :: istat  ! status of allocate statements
         334             :
         335           0 :   if (.not.allocated(alatm)) then
         336           0 :     allocate(alatm(pcols,begchunk:endchunk),stat=istat)
         337           0 :     if (istat /= 0) then
         338           0 :       write(iulog,"('>>> allocate_arrays: allocate of alatm failed: istat=',i5)") istat
         339           0 :       call endrun
         340             :     end if
         341             :   end if
         342             :
         343           0 :   if (.not.allocated(alonm)) then
         344           0 :     allocate(alonm(pcols,begchunk:endchunk),stat=istat)
         345           0 :     if (istat /= 0) then
         346           0 :       write(iulog,"('>>> allocate_arrays: allocate of alonm failed: istat=',i5)") istat
         347           0 :       call endrun
         348             :     end if
         349             :   end if
         350             :
         351           0 :   if (.not.allocated(bnorth)) then
         352           0 :     allocate(bnorth(pcols,begchunk:endchunk),stat=istat)
         353           0 :     if (istat /= 0) then
         354           0 :       write(iulog,"('>>> allocate_arrays: allocate of bnorth failed: istat=',i5)") istat
         355           0 :       call endrun
         356             :     end if
         357             :   end if
         358             :
         359           0 :   if (.not.allocated(beast)) then
         360           0 :     allocate(beast(pcols,begchunk:endchunk),stat=istat)
         361           0 :     if (istat /= 0) then
         362           0 :       write(iulog,"('>>> allocate_arrays: allocate of beast failed: istat=',i5)") istat
         363           0 :       call endrun
         364             :     end if
         365             :   end if
         366             :
         367           0 :   if (.not.allocated(bdown)) then
         368           0 :     allocate(bdown(pcols,begchunk:endchunk),stat=istat)
         369           0 :     if (istat /= 0) then
         370           0 :       write(iulog,"('>>> allocate_arrays: allocate of bdown failed: istat=',i5)") istat
         371           0 :       call endrun
         372             :     end if
         373             :   end if
         374             :
         375           0 :   if (.not.allocated(bmag)) then
         376           0 :     allocate(bmag(pcols,begchunk:endchunk),stat=istat)
         377           0 :     if (istat /= 0) then
         378           0 :       write(iulog,"('>>> allocate_arrays: allocate of bmag failed: istat=',i5)") istat
         379           0 :       call endrun
         380             :     end if
         381             :   end if
         382           0 :   if (.not.allocated(d1vec)) then
         383           0 :     allocate(d1vec(3,pcols,begchunk:endchunk),stat=istat)
         384           0 :     if (istat /= 0) then
         385           0 :       write(iulog,"('>>> allocate_arrays: allocate of d1vec failed: istat=',i5)") istat
         386           0 :       call endrun
         387             :     endif
         388             :   endif
         389             :
         390           0 :   if (.not.allocated(d2vec)) then
         391           0 :     allocate(d2vec(3,pcols,begchunk:endchunk),stat=istat)
         392           0 :     if (istat /= 0) then
         393           0 :       write(iulog,"('>>> allocate_arrays: allocate of d2vec failed: istat=',i5)") istat
         394           0 :       call endrun
         395             :     endif
         396             :   endif
         397             :
         398           0 : end subroutine allocate_arrays
         399             :
         400             : end module mo_apex

Generated by: LCOV version 1.14