LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - calc_pressure.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 8 34 23.5 %
Date: 2024-12-17 17:57:11 Functions: 1 2 50.0 %

          Line data    Source code
       1             : !===============================================================================
       2             : module calc_pressure
       3             : 
       4             :   implicit none
       5             : 
       6             :   public :: init_pressure,   & ! Procedure(s)
       7             :             calculate_thvm
       8             : 
       9             :   private ! default scope
      10             : 
      11             :   contains
      12             : 
      13             :   !=============================================================================
      14           0 :   subroutine init_pressure( gr, thvm, p_sfc, &
      15           0 :                             p_in_Pa, exner, p_in_Pa_zm, exner_zm )
      16             : 
      17             :     ! Description:
      18             :     ! Calculates the initial pressure according to the hydrostatic
      19             :     ! approximation.  Combining the moist equation of state and the hydrostatic
      20             :     ! approximation, the change of pressure with respect to height can be
      21             :     ! calculated based on theta_v, such that:
      22             :     !
      23             :     ! dp/dz = - p * grav / ( Rd * theta_v * exner );
      24             :     !
      25             :     ! where exner = ( p / p0 )^(Rd/Cp);
      26             :     !
      27             :     ! and where p0 is a reference pressure of 100000 Pa.
      28             :     !
      29             :     ! The integral equation is set up to integrate over p on the left-hand side
      30             :     ! and integrate over z on the right-hand side.  The equation is:
      31             :     !
      32             :     ! INT(p1:p2) p^(Rd/Cp-1) dp
      33             :     ! = - p0^(Rd/Cp) * ( grav / Rd ) * INT(z1:z2) (1/thvm) dz.
      34             :     !
      35             :     ! The value of mean theta_v (thvm) is calculated at each thermodynamic grid
      36             :     ! level, and linear interpolation is used in the integral equation for all
      37             :     ! altitude in-between successive thermodynamic levels, such that:
      38             :     !
      39             :     ! thvm(z) = ( ( thvm2 - thvm1 ) / ( z2 - z1 ) ) * ( z - z1 ) + thvm1.
      40             :     !
      41             :     ! The integrals are solved, and the results for pressure can be rewritten
      42             :     ! in terms of exner, such that:
      43             :     !
      44             :     ! exner2 - exner1
      45             :     !   | - ( grav / Cp )
      46             :     !   |   * ( ( z2 - z1 ) / ( thvm2 - thvm1 ) ) * ln( thvm2 / thvm1 );
      47             :     ! = | where thvm2 /= thvm1;
      48             :     !   |
      49             :     !   | - ( grav / Cp ) * ( z2 - z1 ) / thvm; where thvm2 = thvm1 (= thvm).
      50             :     ! 
      51             :     ! The value of pressure (exner) can be calculated using the above equation
      52             :     ! at all vertical levels once the value of pressure is known at one level.
      53             :     ! Since the surface pressure is known at the initial time, that allows
      54             :     ! pressure to be calculated for the rest of the vertical profile.
      55             : 
      56             :     ! References:
      57             :     !-----------------------------------------------------------------------
      58             : 
      59             :     use grid_class, only: &
      60             :         grid, & ! Type
      61             :         zt2zm    ! Procedure(s)
      62             : 
      63             :     use constants_clubb, only: &
      64             :         one,   & ! 1
      65             :         Cp,    & ! Specific heat of dry air                    [J/(kg K)]
      66             :         kappa, & ! Rd/Cp                                       [-]
      67             :         p0,    & ! Reference pressure of 100000 Pa             [Pa]
      68             :         grav,  & ! Acceleration of gravity (9.81 m/s^2)        [m/s^2]
      69             :         zero_threshold
      70             : 
      71             :     use clubb_precision, only: &
      72             :         core_rknd    ! Variable(s)
      73             : 
      74             :     implicit none
      75             : 
      76             :     type (grid), target, intent(in) :: gr
      77             : 
      78             :     ! Input Variables
      79             :     real( kind = core_rknd ), dimension(gr%nz), intent(in) :: &
      80             :       thvm    ! Mean theta_v (thermodynamic levels)                [K]
      81             : 
      82             :     real( kind = core_rknd ), intent(in) :: &
      83             :       p_sfc    ! Surface pressure                                   [Pa]
      84             : 
      85             :     ! Output Variables
      86             :     real( kind = core_rknd ), dimension(gr%nz), intent(out) :: &
      87             :       p_in_Pa,    & ! Pressure (thermodynamic levels)        [Pa]
      88             :       exner,      & ! Exner function (thermodynamic levels)  [-]
      89             :       p_in_Pa_zm, & ! Pressure on momentum levels            [Pa]
      90             :       exner_zm      ! Exner function on momentum levels      [-]
      91             : 
      92             :     ! Local Variables
      93             :     real( kind = core_rknd ), dimension(gr%nz) :: &
      94           0 :       thvm_zm    ! Mean theta_v interpolated to momentum grid levels  [K]
      95             : 
      96             :     real( kind = core_rknd ), parameter :: &
      97             :       g_ov_Cp = grav / Cp  ! g / Cp  [K/m]
      98             : 
      99             :     integer :: k  ! Vertical level index
     100             : 
     101             : 
     102             :     ! The pressure (and exner) at the lowest momentum level is the surface
     103             :     ! pressure (and exner based on the surface pressure).
     104           0 :     p_in_Pa_zm(1) = p_sfc
     105           0 :     exner_zm(1) = ( p_sfc / p0 )**kappa
     106             : 
     107             :     ! Set the pressure (and exner) at the lowest thermodynamic level, which is
     108             :     ! below the model lower boundary, to their values at the model lower
     109             :     ! boundary or surface.
     110           0 :     p_in_Pa(1) = p_in_Pa_zm(1)
     111           0 :     exner(1) = exner_zm(1)
     112             : 
     113             :     ! Interpolate theta_v to momentum levels.
     114           0 :     thvm_zm = zt2zm( gr, thvm, zero_threshold )
     115             : 
     116             :     ! Calculate exner at all other thermodynamic and momentum grid levels.
     117             :     ! exner2
     118             :     ! = exner1
     119             :     !     | ( grav / Cp )
     120             :     !     | * ( ( z2 - z1 ) / ( thvm2 - thvm1 ) ) * ln( thvm2 / thvm1 );
     121             :     !   - | where thvm2 /= thvm1;
     122             :     !     |
     123             :     !     | ( grav / Cp ) * ( z2 - z1 ) / thvm; where thvm2 = thvm1 (= thvm).
     124             : 
     125             :     ! Calculate exner at thermodynamic level 2 (first thermodynamic grid level
     126             :     ! that is above the lower boundary).
     127           0 :     if ( abs( thvm(2) - thvm_zm(1) ) > epsilon( thvm ) * thvm(2) ) then
     128             : 
     129             :        exner(2) &
     130             :        = exner_zm(1) &
     131           0 :          - g_ov_Cp * ( gr%zt(1,2) - gr%zm(1,1) ) / ( thvm(2) - thvm_zm(1) ) &
     132           0 :            * log( thvm(2) / thvm_zm(1) )
     133             : 
     134             :     else ! thvm(2) = thvm_zm(1)
     135             : 
     136           0 :        exner(2) = exner_zm(1) - g_ov_Cp * ( gr%zt(1,2) - gr%zm(1,1) ) / thvm(2)
     137             : 
     138             :     endif
     139             : 
     140             :     ! Calculate pressure on thermodynamic levels.
     141           0 :     p_in_Pa(2) = p0 * exner(2)**(one/kappa)
     142             : 
     143             :     ! Loop over all other thermodynamic vertical grid levels.
     144           0 :     do k = 3, gr%nz, 1
     145             : 
     146             :        ! Calculate exner on thermodynamic levels.
     147           0 :        if ( abs( thvm(k) - thvm(k-1) ) > epsilon( thvm ) * thvm(k) ) then
     148             : 
     149             :           exner(k) &
     150             :           = exner(k-1) &
     151           0 :             - g_ov_Cp * ( gr%zt(1,k) - gr%zt(1,k-1) ) / ( thvm(k) - thvm(k-1) ) &
     152           0 :               * log( thvm(k) / thvm(k-1) )
     153             : 
     154             :        else ! thvm(k+1) = thvm(k)
     155             : 
     156           0 :           exner(k) = exner(k-1) - g_ov_Cp * ( gr%zt(1,k) - gr%zt(1,k-1) ) / thvm(k)
     157             : 
     158             :        endif
     159             : 
     160             :        ! Calculate pressure on thermodynamic levels.
     161           0 :        p_in_Pa(k) = p0 * exner(k)**(one/kappa)
     162             : 
     163             :     enddo ! k = 2, gr%nz, 1
     164             : 
     165             :     ! Loop over all momentum grid levels.
     166           0 :     do k = 2, gr%nz, 1
     167             : 
     168             :        ! Calculate exner on momentum levels.
     169           0 :        if ( abs( thvm(k) - thvm_zm(k) ) > epsilon( thvm ) * thvm(k) ) then
     170             : 
     171             :           exner_zm(k) &
     172             :           = exner(k) &
     173           0 :             - g_ov_Cp * ( gr%zm(1,k) - gr%zt(1,k) ) / ( thvm_zm(k) - thvm(k) ) &
     174           0 :               * log( thvm_zm(k) / thvm(k) )
     175             : 
     176             :        else ! thvm(k) = thvm_zm(k)
     177             : 
     178             :           exner_zm(k) &
     179           0 :           = exner(k) - g_ov_Cp * ( gr%zm(1,k) - gr%zt(1,k) ) / thvm_zm(k)
     180             : 
     181             :        endif
     182             : 
     183             :        ! Calculate pressure on momentum levels.
     184           0 :        p_in_Pa_zm(k) = p0 * exner_zm(k)**(one/kappa)
     185             : 
     186             :     enddo ! k = 2, gr%nz, 1
     187             : 
     188             : 
     189           0 :     return
     190             : 
     191             :   end subroutine init_pressure
     192             : 
     193             :   !=============================================================================
     194     8935056 :   subroutine calculate_thvm( nz, ngrdcol, &
     195     8935056 :                              thlm, rtm, rcm, exner, thv_ds_zt, &
     196     8935056 :                              thvm) 
     197             :   ! Description:
     198             :   ! Calculates mean theta_v based on a linearized approximation to the theta_v
     199             :   ! equation.  This equation also includes liquid water loading.
     200             : 
     201             :   ! References:
     202             :   !-----------------------------------------------------------------------
     203             : 
     204             :     use constants_clubb, only: &
     205             :         Lv,  & ! Latent Heat of Vaporizaion    [J/kg]
     206             :         Cp,  & ! Specific Heat of Dry Air      [J/(kg K)]
     207             :         ep1, & ! Rv/Rd - 1                     [-]
     208             :         ep2    ! Rv/Rd                         [-]
     209             : 
     210             :     use clubb_precision, only: &
     211             :         core_rknd
     212             : 
     213             :     implicit none
     214             : 
     215             :     ! ---------------------------- Input Variables ----------------------------
     216             :     integer, intent(in) :: &
     217             :       nz, &
     218             :       ngrdcol
     219             : 
     220             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in) :: &
     221             :       thlm,      & ! Mean theta_l (thermodynamic levels)          [K]
     222             :       rtm,       & ! Mean total water (thermodynamic levels)      [kg/kg]
     223             :       rcm,       & ! Mean cloud water (thermodynamic levels)      [kg/kg]
     224             :       exner,     & ! Exner function (thermodynamic levels)        [-]
     225             :       thv_ds_zt    ! Reference theta_v on thermodynamic levels    [K]
     226             : 
     227             :     ! ---------------------------- Return Variable ----------------------------
     228             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(out) :: &
     229             :       thvm    ! Mean theta_v (thermodynamic levels)    [K]
     230             : 
     231             :     ! ---------------------------- Local Variables ----------------------------
     232             :     integer :: i, k
     233             : 
     234             :     ! ---------------------------- Begin Code ----------------------------
     235             : 
     236             :     !$acc data copyin( thlm, rtm, rcm, exner, thv_ds_zt ) & 
     237             :     !$acc     copyout( thvm )
     238             : 
     239             :     ! Calculate mean theta_v
     240             :     !$acc parallel loop gang vector collapse(2) default(present)
     241   768414816 :     do k = 1, nz
     242 12690480816 :       do i = 1, ngrdcol
     243 23844132000 :         thvm(i,k) = thlm(i,k) + ep1 * thv_ds_zt(i,k) * rtm(i,k) &
     244 36525677760 :                + ( Lv / ( Cp * exner(i,k) ) - ep2 * thv_ds_zt(i,k) ) * rcm(i,k)
     245             :       end do
     246             :     end do
     247             :     !$acc end parallel loop
     248             :     
     249             :     !$acc end data
     250             : 
     251     8935056 :     return
     252             : 
     253             :   end subroutine calculate_thvm
     254             :   !=============================================================================
     255             : 
     256             : end module calc_pressure 

Generated by: LCOV version 1.14