geometry.f90 Source File

!> author: Sam Hatfield, Fred Kucharski, Franco Molteni
!  date: 01/05/2019
!  For storing all variables related to the model's grid space.

module geometry
    use types, only : p
    use params

    implicit none

    private
    public ModGeometry_t

    !> Geometry module object
    type ModGeometry_t
        logical :: mod_geometry_initialized = .false.

        ! Vertical level parameters
        real(p), allocatable, dimension(:) :: hsg !! Half sigma levels. Dimension(kx + 1).
        real(p), allocatable, dimension(:) :: dhs   !! Sigma level thicknesses. Dimension(kx).
        real(p), allocatable, dimension(:) :: fsg   !! Full sigma levels. Dimension(kx).
        real(p), allocatable, dimension(:) :: dhsr  !! 1/(2*sigma level thicknesses). Dimension(kx).
        real(p), allocatable, dimension(:) :: fsgr  !! akap/(2*full sigma levels). Dimension(kx).

        ! Functions of latitude and longitude
        real(p), allocatable, dimension(:) :: radang   !! Latitudes in radians. Dimension(il).
        real(p), allocatable, dimension(:) :: coriol   !! Coriolis parameter as a function of latitude. Dimension(il).
        real(p), allocatable, dimension(:) :: sia      !! sine(latitude). Dimension(il).
        real(p), allocatable, dimension(:) :: coa      !! cosine(latitude). Dimension(il).
        real(p), allocatable, dimension(:) :: sia_half !! sine(latitude) over one hemisphere only. Dimension(iy).
        real(p), allocatable, dimension(:) :: coa_half !! cosine(latitude) over one hemisphere only. Dimension(il).
        real(p), allocatable, dimension(:) :: cosg     !! Same as coa (TODO: remove). Dimension(il).
        real(p), allocatable, dimension(:) :: cosgr    !! 1/coa. Dimension(il).
        real(p), allocatable, dimension(:) :: cosgr2   !! 1/coa^2. Dimension(il).

        ! Functions of sigma and latitude (initial. in INPHYS)
        real(p), allocatable, dimension(:) :: sigl   !! Logarithm of full-level sigma. Dimension(kx).
        real(p), allocatable, dimension(:) :: sigh   !! Half-level sigma. Dimension(0:kx).

        !> g/(d_sigma p0) : to convert fluxes of u,v,q into d(u,v,q)/dt .
        !! Dimension(kx).
        real(p), allocatable, dimension(:) :: grdsig

        !> g/(d_sigma p0 c_p): to convert energy fluxes into dT/dt
        !! Dimension(kx).
        real(p), allocatable, dimension(:) :: grdscp

        !> Weights for vertical interpolation. Dimension(kx, 2).
        real(p), allocatable, dimension(:, :) :: wvi

    contains
        procedure :: initialize => ModGeometry_initialize
        procedure :: delete => ModGeometry_delete
        !        procedure :: fourier => ModFourier_fourier
        !        procedure :: fourier_inv => ModFourier_fourier_inv
    end type


contains
    !> Initializes all of the model geometry variables.
    subroutine ModGeometry_initialize(this)
        use physical_constants, only : akap, omega, grav, cp, p0
        class(ModGeometry_t), intent(inout) :: this

        integer j, jj, k

        if (this%mod_geometry_initialized) then
            !Do nothing, the module is already initialized.
            return
        end if

        ! Allocate variables -----------------------------------------------
        allocate(this%hsg(kx + 1))
        allocate(this%dhs(kx), this%fsg(kx), this%dhsr(kx), this%fsgr(kx))
        allocate(this%radang(il), this%coriol(il), this%sia(il), this%coa(il))
        allocate(this%sia_half(iy), this%coa_half(il))
        allocate(this%cosg(il), this%cosgr(il), this%cosgr2(il))
        allocate(this%sigl(kx), this%sigh(0:kx))
        allocate(this%grdsig(kx), this%grdscp(kx), this%wvi(kx, 2))
        ! End allocation ---------------------------------------------------

        ! Definition of model levels
        ! Half (vertical velocity) levels
        if (kx == 5) then
            this%hsg(:6) = (/0.000, 0.150, 0.350, 0.650, 0.900, 1.000/)
        else if (kx == 7) then
            this%hsg(:8) = (/0.020, 0.140, 0.260, 0.420, 0.600, 0.770, 0.900, 1.000/)
        else if (kx == 8) then
            this%hsg(:9) = (/0.000, 0.050, 0.140, 0.260, 0.420, 0.600, 0.770, 0.900, 1.000/)
        end if

        ! Layer thicknesses and full (u,v,T) levels
        do k = 1, kx
            this%dhs(k) = this%hsg(k + 1) - this%hsg(k)
            this%fsg(k) = 0.5 * (this%hsg(k + 1) + this%hsg(k))
        end do

        ! Additional functions of sigma
        do k = 1, kx
            this%dhsr(k) = 0.5 / this%dhs(k)
            this%fsgr(k) = akap / (2. * this%fsg(k))
        end do

        ! Horizontal functions

        ! Latitudes and functions of latitude
        ! NB: J=1 is Southernmost point!
        do j = 1, iy
            jj = il + 1 - j
            this%sia_half(j) = cos(3.141592654 * (j - 0.25) / (il + 0.5))
            this%coa_half(j) = sqrt(1.0 - this%sia_half(j)**2.0)
            this%sia(j) = -this%sia_half(j)
            this%sia(jj) = this%sia_half(j)
            this%coa(j) = this%coa_half(j)
            this%coa(jj) = this%coa_half(j)
            this%radang(j) = -asin(this%sia_half(j))
            this%radang(jj) = asin(this%sia_half(j))
        end do

        ! Expand cosine and its reciprocal to cover both hemispheres
        do j = 1, iy
            jj = il + 1 - j
            this%cosg(j) = this%coa_half(j)
            this%cosg(jj) = this%coa_half(j)
            this%cosgr(j) = 1. / this%coa_half(j)
            this%cosgr(jj) = 1. / this%coa_half(j)
            this%cosgr2(j) = 1. / (this%coa_half(j) * this%coa_half(j))
            this%cosgr2(jj) = 1. / (this%coa_half(j) * this%coa_half(j))
        end do

        this%coriol = 2.0 * omega * this%sia

        ! 1.2 Functions of sigma and latitude
        this%sigh(0) = this%hsg(1)

        do k = 1, kx
            this%sigl(k) = log(this%fsg(k))
            this%sigh(k) = this%hsg(k + 1)
            this%grdsig(k) = grav / (this%dhs(k) * p0)
            this%grdscp(k) = this%grdsig(k) / cp
        end do

        ! Weights for vertical interpolation at half-levels(1,kx) and surface
        ! Note that for phys.par. half-lev(k) is between full-lev k and k+1
        ! Fhalf(k) = Ffull(k)+this%wvi(K,2)*(Ffull(k+1)-Ffull(k))
        ! Fsurf = Ffull(kx)+this%wvi(kx,2)*(Ffull(kx)-Ffull(kx-1))
        do k = 1, kx - 1
            this%wvi(k, 1) = 1. / (this%sigl(k + 1) - this%sigl(k))
            this%wvi(k, 2) = (log(this%sigh(k)) - this%sigl(k)) * this%wvi(k, 1)
        end do

        this%wvi(kx, 1) = 0.
        this%wvi(kx, 2) = (log(0.99) - this%sigl(kx)) * this%wvi(kx - 1, 1)

        this%mod_geometry_initialized = .true.
    end subroutine

    !> Delete the Geometry instance
    subroutine ModGeometry_delete(this)
        class(ModGeometry_t), intent(inout) :: this

        if (this%mod_geometry_initialized) then
            ! Allocate variables -----------------------------------------------
            deallocate(this%hsg, this%dhs, this%fsg, this%dhsr, this%fsgr)
            deallocate(this%radang, this%coriol, this%sia, this%coa)
            deallocate(this%sia_half, this%coa_half, this%cosg, this%cosgr, this%cosgr2)
            deallocate(this%sigl, this%sigh, this%grdsig, this%grdscp, this%wvi)
            ! End allocation ---------------------------------------------------
            this%mod_geometry_initialized = .false.
        end if
    end subroutine
end module