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 22:39:59 |  | 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