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