LCOV - code coverage report
Current view: | top level - physics/cam - lunar_tides.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 20 | 44 | 45.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 3 | 100.0 %

              Line data    Source code
           1             : ! Module to add lunar tide forcing for the M2
           2             : ! gravitational lunar tide based on the empirical
           3             : ! formula of Champman and Lindzen (1970) and as
           4             : ! implemented in Pedatella, Liu, and Richmond (2012, JGR)
           5             :
           6             : module lunar_tides
           7             :   use shr_kind_mod,   only: r8=>shr_kind_r8
           8             :   use physics_types,  only: physics_state, physics_ptend, physics_ptend_init
           9             :   use phys_control,   only: use_simple_phys
          10             :   use spmd_utils,     only: masterproc
          11             :   use cam_abortutils, only: endrun
          12             :   use cam_logfile,    only: iulog
          13             :
          14             :   implicit none
          15             :   private
          16             :
          17             :   public :: lunar_tides_readnl
          18             :   public :: lunar_tides_init
          19             :   public :: lunar_tides_tend
          20             :
          21             :   ! lunar tides forcing option
          22             :   logical  :: apply_lunar_tides = .false.
          23             :
          24             : contains
          25             :
          26             :   !==========================================================================
          27             :   !==========================================================================
          28     1490712 :   subroutine lunar_tides_readnl(nlfile)
          29             :     use namelist_utils, only: find_group_name
          30             :     use units,          only: getunit, freeunit
          31             :     use spmd_utils,     only: mpicom, mstrid=>masterprocid, mpi_logical
          32             :
          33             :     ! File containing namelist input.
          34             :     character(len=*), intent(in) :: nlfile
          35             :
          36             :     ! Local variables
          37             :     integer :: unitn, ierr
          38             :     character(len=*), parameter :: sub = 'lunar_tides_readnl'
          39             :
          40             :     namelist /lunar_tides_opts/ apply_lunar_tides
          41             :
          42        1536 :     if (use_simple_phys) return
          43             :
          44        1536 :     if (masterproc) then
          45           2 :        unitn = getunit()
          46           2 :        open( unitn, file=trim(nlfile), status='old' )
          47           2 :        call find_group_name(unitn, 'lunar_tides_opts', status=ierr)
          48           2 :        if (ierr == 0) then
          49           0 :           read(unitn, lunar_tides_opts, iostat=ierr)
          50           0 :           if (ierr /= 0) then
          51           0 :              call endrun(sub // ':: ERROR reading namelist')
          52             :           end if
          53             :        end if
          54           2 :        close(unitn)
          55           2 :        call freeunit(unitn)
          56             :     end if
          57             :
          58        1536 :     call mpi_bcast(apply_lunar_tides, 1, mpi_logical, mstrid, mpicom, ierr)
          59        1536 :     if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: apply_lunar_tides")
          60             :
          61        1536 :     if (masterproc) then
          62           2 :        write(iulog,*) 'lunar_tides_readnl: apply_lunar_tides: ',apply_lunar_tides
          63             :     end if
          64             :
          65             :   end subroutine lunar_tides_readnl
          66             :
          67             :   !==========================================================================
          68             :   !==========================================================================
          69        1536 :   subroutine lunar_tides_init
          70             :     use cam_history, only: addfld
          71             :     use time_manager,only: timemgr_get_calendar_cf
          72             :
          73        1536 :     if (apply_lunar_tides) then
          74           0 :        if (timemgr_get_calendar_cf().ne.'gregorian') then
          75           0 :           call endrun('lunar_tides_init: calendar must be gregorian')
          76             :        endif
          77           0 :        call addfld('UT_LUNAR', (/ 'lev' /), 'A','m/s2','Zonal wind tendency due to lunar tides')
          78           0 :        call addfld('VT_LUNAR', (/ 'lev' /), 'A','m/s2','Meridional wind tendency due to lunar tides')
          79             :     end if
          80             :
          81        1536 :   end subroutine lunar_tides_init
          82             :
          83             :   !==========================================================================
          84             :   !==========================================================================
          85    62545392 :   subroutine lunar_tides_tend( state, ptend )
          86        1536 :     use time_manager, only: get_curr_date, get_julday
          87             :     use physconst,    only: pi, rearth
          88             :     use ppgrid,       only: pver
          89             :     use cam_history,  only: outfld
          90             :
          91             :     type(physics_state), intent(in) :: state
          92             :     type(physics_ptend), intent(out):: ptend
          93             :
          94             :     integer  :: tod,yr,mm,dd
          95             :     real(r8) :: jd,nu,lt,lun_lt
          96             :
          97             :     integer :: i, k
          98             :
          99             :     real(r8), parameter :: deg2hrs = 1._r8/15._r8
         100             :     real(r8), parameter :: rad2deg = 180._r8/pi
         101             :     real(r8), parameter :: rad2hrs = rad2deg*deg2hrs
         102             :     real(r8), parameter :: tod2hrs = 24._r8/86400._r8
         103             :     real(r8), parameter :: hrs2rad = 1._r8/rad2hrs
         104             :
         105     1489176 :     if (apply_lunar_tides) then
         106             :
         107           0 :        call physics_ptend_init(ptend, state%psetcols, "Lunar Tides", lu=.true., lv=.true. )
         108             :
         109             :        ! calculate the current date:
         110           0 :        call get_curr_date(yr,mm,dd,tod)
         111             :        ! convert date to Julian centuries
         112           0 :        jd = get_julday(yr,mm,dd,tod)
         113             :        ! calculation relies on time from noon on December 31, 1899, so
         114             :        ! subtract 2415020, which corresponds to the Julian date for Dec. 31 1899.
         115           0 :        jd = jd - 2415020._r8
         116           0 :        jd = jd / 36525._r8 ! convert to julian centuries
         117             :
         118             :        ! Calculate the lunar local time (nu) based on the the time
         119             :        ! in Julian centuries using the formula given in Chapman and Lindzen (1970)
         120           0 :        nu = -9.26009_r8 + 445267.12165_r8*jd+0.00168_r8*jd*jd !nu in degrees
         121             :
         122           0 :        do i=1,state%ncol
         123             :           ! solar local time (hours)
         124           0 :           lt = real(tod,kind=r8)*tod2hrs + state%lon(i)*rad2hrs
         125             :
         126             :           ! lunar local time
         127           0 :           lun_lt = lt - nu*deg2hrs ! hours
         128           0 :           lun_lt = lun_lt*hrs2rad ! radians
         129             :
         130           0 :           do k=1,pver
         131             :              ! Calculate the M2 lunar tide forcing in the zonal and meridional directions.
         132             :              ! The forcing is calculated based on the gradient of the M2 tidal
         133             :              ! potential, which is given in Chapman and Lindzen (1970).
         134             :              ! Additional details on the derivation of the forcing are in
         135             :              ! Pedatella, Liu, and Richmond (2012)
         136           0 :              ptend%u(i,k) = (-1._r8/((state%zm(i,k)+rearth)*cos(state%lat(i))))*2.456_r8*3._r8 *  &
         137           0 :                   ((state%zm(i,k)+rearth)/rearth)**2*cos(state%lat(i))*cos(state%lat(i))*2._r8*sin(2._r8*lun_lt)
         138           0 :              ptend%v(i,k) = (1._r8/(state%zm(i,k)+rearth))*2.456_r8*3._r8 * &
         139           0 :                   ((state%zm(i,k)+rearth)/rearth)**2*cos(2._r8*lun_lt)*2._r8*cos(state%lat(i))*sin(state%lat(i))
         140             :           end do
         141             :        end do
         142             :
         143           0 :        call outfld('UT_LUNAR', ptend%u(:state%ncol,:), state%ncol, state%lchnk)
         144           0 :        call outfld('VT_LUNAR', ptend%v(:state%ncol,:), state%ncol, state%lchnk)
         145             :
         146             :     end if
         147             :
         148     1489176 :   end subroutine lunar_tides_tend
         149             :
         150             : end module lunar_tides

Generated by: LCOV version 1.14