sources_planetheating.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_planetheating.f90                                         #
!#                                                                           #
!# Copyright (C) 2013-2024                                                   #
!# Jannes Klee <jklee@astrophysik.uni-kiel.de>                               #
!# Tobias Illenseer <tillense@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_planetheating_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 = "heating of planetary atmosphere"
  REAL, PARAMETER   :: au      = 1.49597870691e+11    ! astr. unit [m]     !
  !--------------------------------------------------------------------------!
  TYPE, EXTENDS(sources_cooling) :: sources_planetheating
    TYPE(marray_base), ALLOCATABLE &
                      :: cos1,sin1
    REAL              :: radflux
    REAL              :: albedo
    REAL              :: sm_axis
    REAL              :: excentricity
    REAL              :: true_anomaly
    REAL              :: omegasun
    REAL              :: year
    REAL              :: theta0, phi0
  CONTAINS
    PROCEDURE :: initsources
!     PROCEDURE :: SetOutput
    PROCEDURE :: infosources
    PROCEDURE :: updatecooling => updateplanetheating ! overwrite inherited method from sources_cooling
    final :: finalize
  END TYPE
  PUBLIC :: &
       ! types
       sources_planetheating
  !--------------------------------------------------------------------------!
 
CONTAINS
 
 
  SUBROUTINE initsources(this,Mesh,Physics,Fluxes,config,IO)
    USE physics_euler_mod, ONLY : physics_euler
    USE geometry_spherical_mod, ONLY : geometry_spherical
    USE geometry_spherical_planet_mod, ONLY : geometry_spherical_planet
    USE constants_si_mod, ONLY : constants_si
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetheating), 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
    !------------------------------------------------------------------------!
    INTEGER                         :: err
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitSources_base(stype,source_name)
 
    CALL this%Info("############################ ATTENTION ############################")
    CALL this%Warning("sources_planetheating::InitSources","Untested module, use with care!")
 
    ! some sanity checks
    SELECT TYPE (phys => physics)
    CLASS IS(physics_euler)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("InitSources_planetheating","physics not supported")
    END SELECT
 
    SELECT TYPE (const => physics%constants)
    CLASS IS(constants_si)
      ! do nothing
    CLASS DEFAULT
      CALL this%Error("InitSources_planetheating","only SI units supported")
    END SELECT
 
    ALLOCATE(this%Q,this%cos1,this%sin1)
    this%Q     = marray_base()
    this%cos1  = marray_base(2) !TODO check if 2 is really necessary
    this%sin1  = marray_base(2) !TODO check if 2 is really necessary
 
    ! Courant number, i.e. safety factor for numerical stability
    CALL getattr(config, "cvis", this%cvis, 0.1)
    ! intensity of the planets star at 1 AU
    CALL getattr(config, "radflux", this%radflux)
    ! albedo of the planet
    CALL getattr(config, "albedo", this%albedo)
    ! nummerical excentricity
    CALL getattr(config, "excentricity", this%excentricity, 0.0)
    ! semi-major-axis
    CALL getattr(config, "semimajoraxis", this%sm_axis)
    ! day-night omega
    CALL getattr(config, "omegasun", this%omegasun, 0.0)
    ! year
    CALL getattr(config, "year", this%year)
    ! beginning points
    CALL getattr(config, "theta0", this%theta0, 0.0)
    CALL getattr(config, "phi0", this%phi0, 0.0)
 
    SELECT TYPE (mgeo => mesh%Geometry)
    TYPE IS(geometry_spherical)
      ! TODO: ATTENTION CHECK THAT THIS REGION IS CORRECT AFTER TRANSITION TO NEW FOSITE
      this%cos1%data4d(:,:,:,1) = cos(mesh%bcenter(:,:,:,2))
!      this%cos1(:,:,2) = COS(Mesh%bcenter(:,:,2))
      this%sin1%data4d(:,:,:,1) = sin(mesh%bcenter(:,:,:,2))
!      this%sin1(:,:,2) = SIN(Mesh%bcenter(:,:,2))
    TYPE IS(geometry_spherical_planet)
      ! TODO: ATTENTION CHECK THAT THIS REGION IS CORRECT AFTER TRANSITION TO NEW FOSITE
      this%cos1%data4d(:,:,:,1) = cos(mesh%bcenter(:,:,:,1))
!      this%cos1(:,:,2) = COS(Mesh%bcenter(:,:,2))
      this%sin1%data4d(:,:,:,1) = sin(mesh%bcenter(:,:,:,1))
!      this%sin1(:,:,2) = SIN(Mesh%bcenter(:,:,2))
    CLASS DEFAULT
      CALL this%Error("InitSources_planetheating","only spherical geometry allowed")
    END SELECT
 
    ! set initial time < 0
    ! TODO: Check were this time is set? Is it correct? Why not timedisc time?
    this%time = -1.0
 
    ! initialize arrays
    this%Q%data1d(:) = 0.0
 
    ! register ouput arrays
    CALL this%SetOutput(mesh,config,io)
 
    ! call InitSources in base
    CALL this%InfoSources(mesh)
  END SUBROUTINE initsources
 
  SUBROUTINE setoutput(this,Mesh,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetheating), INTENT(INOUT) :: this
    CLASS(mesh_base), INTENT(IN) :: Mesh
    TYPE(dict_typ), POINTER     :: config,IO
    !------------------------------------------------------------------------!
    INTEGER              :: valwrite
    !------------------------------------------------------------------------!
    ! register heating source term for output
    CALL getattr(config, "output/Qstar", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "Qstar", &
              this%Q%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
  END SUBROUTINE setoutput
 
 
  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetheating), INTENT(IN) :: this
    CLASS(mesh_base),             INTENT(IN) :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: param_str
    REAL              :: tmp_out
    REAL, PARAMETER   :: AU      = 1.49597870691e+11    ! astr. unit [m]     !
    !------------------------------------------------------------------------!
    WRITE (param_str,'(ES10.2)') this%year/(3.6e3*24*365.2425)
    CALL this%Info("            sid. year:         " // trim(param_str) // " yr")
 
    IF (abs(this%omegasun).GT.tiny(this%omegasun)) THEN
       WRITE (param_str,'(ES10.2)') abs(1./(this%omegasun*3.6e3*24.0))
       CALL this%Info("            day:               " // trim(param_str) // " d")
    ELSE
       CALL this%Info("            day:               " // "infinite (tidally locked)")
    END IF
 
    WRITE (param_str,'(ES10.2)') this%sm_axis/au
    CALL this%Info("            semi-major axis:   " // trim(param_str) // " au")
 
    WRITE (param_str,'(ES10.2)') this%excentricity
    CALL this%Info("            excentricity:      " // trim(param_str))
 
    WRITE (param_str,'(ES10.2)') this%sm_axis/au*(1.+0.5*this%excentricity**2)
    CALL this%Info("            mean distance:     " // trim(param_str) // " au")
  END SUBROUTINE infosources
 
 
  SUBROUTINE updateplanetheating(this,Mesh,Physics,time,pvar)
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    USE gravity_binary_mod, ONLY : kepler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_planetheating), 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              :: theta2,phi1,phi2,r,mean_anomaly,true_anomaly
    !------------------------------------------------------------------------!
    ! heating by the star
    IF (time.NE.this%time) THEN
      ! calculate the new distance (and true anomaly) by solving Kepler's equation
      CALL kepler(time,this%year,this%excentricity,this%sm_axis,r,true_anomaly)
 
      DO k=mesh%KMIN,mesh%KMAX
        DO j=mesh%JMIN,mesh%JMAX
          DO i=mesh%IMIN,mesh%IMAX
            !--------------------------------------------------------------!
            phi1=mesh%bcenter(i,j,k,2) + 2.*pi*this%omegasun*this%time
            ! transformation
            theta2 = acos(this%cos1%data4d(i,j,k,1)*cos(this%theta0*&
                     cos(2*pi*this%time/this%year)) + &
                     this%sin1%data4d(i,j,k,1)*sin(this%theta0*&
                     cos(2*pi*this%time/this%year))*cos(phi1-this%phi0))
 
            ! modulo function to prevent phi2>2*PI,phi2<0
            phi2   = modulo(atan2(this%sin1%data4d(i,j,k,1)*sin(phi1-this%phi0),&
                     (-this%cos1%data4d(i,j,k,1)*sin(this%theta0*cos(2.*pi*this%time/&
                     this%year)) + this%sin1%data4d(i,j,k,1)*cos(phi1-this%phi0)*&
                     cos(this%theta0*cos(2*pi*this%time/this%year)))),2*pi)
 
            ! calculate heating source                                     !
            !--------------------------------------------------------------!
            IF (phi2.LE.pi/2..OR.phi2.GE.3./2.*pi) THEN
              ! for sin,cos have a look at spherical trigonometry
              ! the albedo can contain things like scattering etc. (not implemented),
              ! radflux is given at 1 AU
              this%Q%data3d(i,j,k) = (1.-this%albedo)*this%radflux * (au/r)**2 &
                                            *sin(theta2)*cos(phi2)
            ELSE
              this%Q%data3d(i,j,k) = 0.0
            END IF
          END DO
        END DO
      END DO
      this%time = time
    END IF
  END SUBROUTINE updateplanetheating
 
 
  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    TYPE(sources_planetheating), INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    DEALLOCATE(this%cos1,this%sin1)
    ! this%Q is deallocated in parent class destructor
  END SUBROUTINE finalize
 
END MODULE sources_planetheating_mod