sources_planetcooling.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_planetcooling.f90                                         #
!#                                                                           #
!# Copyright (C) 2011-2024                                                   #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!# Jannes Klee      <jklee@astrophysik.uni-kiel.de>                          #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with this program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE sources_planetcooling_mod
  USE sources_base_mod
  USE sources_cooling_mod
  USE physics_base_mod
  USE fluxes_base_mod
  USE mesh_base_mod
  USE marray_compound_mod
  USE marray_base_mod
  USE logging_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
  CHARACTER(LEN=32), PARAMETER :: source_name = "cooling of planetary atmosphere"
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(sources_cooling) :: sources_planetcooling
    TYPE(marray_base), ALLOCATABLE &
                      :: t_s, &         !< temperature at the surface
                         rho_s, &       !< density at the surface
                         p_s
    REAL              :: t_0
    REAL              :: radflux
    REAL              :: albedo
    REAL              :: distance
    REAL              :: gacc
    REAL              :: tau_inf
    REAL              :: const1, const2
  CONTAINS
    PROCEDURE :: initsources
    PROCEDURE :: infosources
    PROCEDURE :: updatecooling
    final :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
  PUBLIC :: &
       ! types
       sources_planetcooling
  !--------------------------------------------------------------------------!
 
CONTAINS
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY : physics_euler
    USE constants_si_mod, ONLY : constants_si
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetcooling), INTENT(INOUT) :: this
    CLASS(mesh_base),    INTENT(IN) :: Mesh
    CLASS(physics_base), INTENT(IN) :: Physics
    CLASS(fluxes_base),  INTENT(IN) :: Fluxes
    TYPE(dict_typ),      POINTER    :: config,IO
    INTEGER           :: stype
    !------------------------------------------------------------------------!
    REAL              :: effradflux, q
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitSources_base(stype,source_name)
 
    CALL this%Info("############################ ATTENTION ############################")
    CALL this%Warning("sources_planetcooling::InitSources","Untested module, use with care!")
 
    ! some sanity checks
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("sources_planetheating::InitSources","physics not supported")
    END SELECT
 
    SELECT TYPE (const => physics%constants)
    CLASS IS(constants_si)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("sources_planetheating::InitSources","only SI units supported")
    END SELECT
 
    ALLOCATE(this%Q,this%T_s)
    this%Q = marray_base()
    this%T_s = marray_base()
 
    ! Courant number, i.e. safety factor for numerical stability
    CALL getattr(config, "cvis", this%cvis, 0.9)
    ! radiant energy flux density of stellar radiation at 1 AU distance
    ! from the central star [W/m^2]
    CALL getattr(config, "radflux", this%radflux)
    ! Bond albedo of the planet
    CALL getattr(config, "albedo", this%albedo,0.0)
    ! temporal mean distance between planet and star
    CALL getattr(config, "mean_distance", this%distance,physics%constants%AU)
    ! equilibrium temperature
    CALL getattr(config,"T_0", this%T_0)
    ! gravitational acceleration
    CALL getattr(config, "gacc", this%gacc)
 
    ! set initial time < 0
    this%time = -1.0
 
    this%Q%data1d(:) = 0.0
    this%T_s%data1d(:)   = this%T_0
 
    ! initial calculation of optical thickness
    SELECT TYPE(phys => physics)
    CLASS IS(physics_euler)
      ! effective mean radiant flux density at mean distance
      effradflux = (1.-this%albedo)*this%radflux*(phys%constants%AU/this%distance)**2
      ! compute the effective optical depth that leads to the equilibrium
      ! mean surface temperature T_0 of the planet
      q = phys%gamma/(phys%gamma-1.0)
      this%tau_inf= ((4*phys%Constants%SB*gamma(1.+4./q)/effradflux)**(0.25) * this%T_0)**q
      ! conversion factor cs^2 -> T
      this%const1 = phys%mu/(phys%gamma*phys%Constants%RG)
      this%const2 = phys%Constants%SB*this%tau_inf**(-4./q) * gamma(1.0+4./q)
    END SELECT
 
    ! register cooling function for output
    CALL this%SetOutput(mesh,config,io)
 
    ! call InitSources in base
    CALL this%InfoSources(mesh)
  END SUBROUTINE initsources
 
 
  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetcooling), INTENT(IN) :: this
    CLASS(mesh_base),             INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: param_str
    !------------------------------------------------------------------------!
    WRITE (param_str,'(ES10.2)') this%radflux
    CALL this%Info("            mean stellar flux: " // trim(param_str) // " W/m^2")
    WRITE (param_str,'(ES10.2)') this%T_0
    CALL this%Info("            mean equil. temp.: " // trim(param_str) // " K")
    WRITE (param_str,'(ES10.2)') this%tau_inf
    CALL this%Info("            optical depth:     " // trim(param_str))
  END SUBROUTINE infosources
 
 
  SUBROUTINE updatecooling(this,Mesh,Physics,time,pvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    USE functions, ONLY : lngamma
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetcooling), INTENT(INOUT) :: this
    CLASS(mesh_base),             INTENT(IN)    :: Mesh
    CLASS(physics_euler),         INTENT(IN)    :: Physics
    REAL,                         INTENT(IN)    :: time
    CLASS(statevector_euler),     INTENT(IN)    :: pvar
    !------------------------------------------------------------------------!
    INTEGER           :: i,j,k
    REAL              :: Qcool
    REAL              :: const1,const2       ! needed due to performance issues
    !------------------------------------------------------------------------!
    ! calculate the cooling function if time has changed
    IF (time.NE.this%time) THEN
      ! calculate planet-surface temperature using integrated
      ! pressure and density
      this%T_s%data1d(:) = this%const1 * pvar%pressure%data1d(:)/pvar%density%data1d(:)
 
      ! calculate surface density and pressure if output is requested
      ! (irrelevant for cooling function)
      IF (ALLOCATED(this%P_s)) this%P_s%data1d(:) = pvar%density%data1d(:)*this%gacc
      IF (ALLOCATED(this%rho_s)) this%rho_s%data1d(:) = physics%mu/(physics%Constants%RG) &
        *this%P_s%data1d(:)/this%T_s%data1d(:)
 
      ! cooling source
      WHERE (mesh%without_ghost_zones%mask1d(:))
        this%Q%data1d(:) = this%const2*this%T_s%data1d(:)**4
      ELSEWHERE
        this%Q%data1d(:) = 0.0
      END WHERE
 
      ! update time
      this%time=time
    END IF
  END SUBROUTINE updatecooling
 
 
  SUBROUTINE setoutput(this,Mesh,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetcooling), INTENT(INOUT) :: this
    CLASS(mesh_base), INTENT(IN) :: Mesh
    TYPE(dict_typ),POINTER  :: config,IO
    !------------------------------------------------------------------------!
    INTEGER              :: valwrite
    !------------------------------------------------------------------------!
 
    !cooling source term
    CALL getattr(config, "output/Qcool", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "Qcool", &
              this%Q%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
 
    ! temperature
    CALL getattr(config, "output/T_s", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "T_s", &
              this%T_s%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
 
    ! surface pressure
    CALL getattr(config, "output/P_s", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
      ALLOCATE(this%P_s)
      this%P_s = marray_base()
      this%P_s%data1d(:)   = 0.0
      CALL setattr(io, "P_s", &
              this%P_s%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    END IF
 
    ! surface density
    CALL getattr(config, "output/rho_s", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
      ALLOCATE(this%rho_s)
      this%rho_s = marray_base()
      this%rho_s%data1d(:) = 0.0
      CALL setattr(io, "rho_s", &
              this%rho_s%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    END IF
  END SUBROUTINE setoutput
 
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_planetcooling), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    DEALLOCATE(this%T_s)
    IF (ALLOCATED(this%P_s)) DEALLOCATE(this%P_s)
    IF (ALLOCATED(this%rho_s)) DEALLOCATE(this%rho_s)
    ! this%Q is deallocated in parent class destructor
  END SUBROUTINE finalize
 
END MODULE sources_planetcooling_mod