LCOV - code coverage report
Current view: | top level - atmos_phys/schemes/tropopause_find - tropopause_find.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 216 | 359 | 60.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 11 | 16 | 68.8 %

              Line data    Source code
           1             : ! This module is used to diagnose the location of the tropopause. Multiple
           2             : ! algorithms are provided, some of which may not be able to identify a
           3             : ! tropopause in all situations. To handle these cases, an analytic
           4             : ! definition and a climatology are provided that can be used to fill in
           5             : ! when the original algorithm fails. The tropopause temperature and
           6             : ! pressure are determined and can be output to the history file.
           7             : !
           8             : ! Author: Charles Bardeen
           9             : ! Created: April, 2009
          10             : !
          11             : ! CCPP-ized: Haipeng Lin, August 2024
          12             : module tropopause_find
          13             :   !---------------------------------------------------------------
          14             :   ! ... variables for the tropopause module
          15             :   !---------------------------------------------------------------
          16             :
          17             :   use ccpp_kinds,           only : kind_phys
          18             :
          19             :   implicit none
          20             :
          21             :   private
          22             :
          23             :   ! CCPP-compliant subroutines
          24             :   public :: tropopause_find_init
          25             :   public :: tropopause_find_run
          26             :
          27             :   ! "Wrapped" routine for use in old CAM for backward compatibility.
          28             :   ! Also called by tropopause_find_run driver routine.
          29             :   public :: tropopause_findWithBackup
          30             :
          31             :   ! Switches for tropopause method
          32             :   public :: TROP_ALG_NONE, TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
          33             :   public :: TROP_ALG_STOBIE, TROP_ALG_HYBSTOB, TROP_ALG_TWMO, TROP_ALG_WMO
          34             :   public :: TROP_ALG_CPP
          35             :   public :: NOTFOUND
          36             :
          37             :   save
          38             :
          39             :   ! These parameters define an enumeration to be used to define the primary
          40             :   ! and backup algorithms to be used with the tropopause_find() method. The
          41             :   ! backup algorithm is meant to provide a solution when the primary algorithm
          42             :   ! fails. The algorithms that can't fail (i.e., always find a tropopause) are:
          43             :   ! TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE, and TROP_ALG_STOBIE.
          44             :   integer, parameter    :: TROP_ALG_NONE      = 1    ! Don't evaluate
          45             :   integer, parameter    :: TROP_ALG_ANALYTIC  = 2    ! Analytic Expression
          46             :   integer, parameter    :: TROP_ALG_CLIMATE   = 3    ! Climatology
          47             :   integer, parameter    :: TROP_ALG_STOBIE    = 4    ! Stobie Algorithm
          48             :   integer, parameter    :: TROP_ALG_TWMO      = 5    ! WMO Definition, Reichler et al. [2003]
          49             :   integer, parameter    :: TROP_ALG_WMO       = 6    ! WMO Definition
          50             :   integer, parameter    :: TROP_ALG_HYBSTOB   = 7    ! Hybrid Stobie Algorithm
          51             :   integer, parameter    :: TROP_ALG_CPP       = 8    ! Cold Point Parabolic
          52             :   integer, parameter    :: TROP_ALG_CHEMTROP  = 9    ! Chemical tropopause
          53             :
          54             :   integer, parameter    :: default_primary    = TROP_ALG_TWMO        ! default primary algorithm
          55             :   integer, parameter    :: default_backup     = TROP_ALG_CLIMATE     ! default backup algorithm
          56             :
          57             :   integer, parameter    :: NOTFOUND = -1
          58             :
          59             :   real(kind_phys), parameter :: ALPHA  = 0.03_kind_phys
          60             :
          61             :   ! physical constants
          62             :   ! These constants are set in module variables rather than as parameters
          63             :   ! to support the aquaplanet mode in which the constants have values determined
          64             :   ! by the experiment protocol
          65             :   real(kind_phys) :: cnst_kap     ! = cappa
          66             :   real(kind_phys) :: cnst_faktor  ! = -gravit/rair
          67             :   real(kind_phys) :: cnst_rga     ! = 1/gravit
          68             :   real(kind_phys) :: cnst_ka1     ! = cnst_kap - 1._kind_phys
          69             :   real(kind_phys) :: cnst_rad2deg ! = 180/pi
          70             :
          71             : !================================================================================================
          72             : contains
          73             : !================================================================================================
          74             :
          75             : !> \section arg_table_tropopause_find_init Argument Table
          76             : !! \htmlinclude tropopause_find_init.html
          77        1536 :   subroutine tropopause_find_init(cappa, rair, gravit, pi, errmsg, errflg)
          78             :
          79             :     real(kind_phys), intent(in)         :: cappa   ! R/Cp
          80             :     real(kind_phys), intent(in)         :: rair    ! Dry air gas constant (J K-1 kg-1)
          81             :     real(kind_phys), intent(in)         :: gravit  ! Gravitational acceleration (m s-2)
          82             :     real(kind_phys), intent(in)         :: pi      ! Pi
          83             :
          84             :     character(len=512), intent(out) :: errmsg
          85             :     integer,            intent(out) :: errflg
          86             :
          87        1536 :     errmsg = ' '
          88        1536 :     errflg = 0
          89             :
          90             :     ! define physical constants
          91        1536 :     cnst_kap     = cappa
          92        1536 :     cnst_faktor  = -gravit/rair
          93        1536 :     cnst_rga     = 1._kind_phys/gravit              ! Reciprocal of gravit (s2 m-1)
          94        1536 :     cnst_ka1     = cnst_kap - 1._kind_phys
          95        1536 :     cnst_rad2deg = 180._kind_phys/pi               ! radians to degrees conversion factor
          96             :
          97        1536 :   end subroutine tropopause_find_init
          98             :
          99             :   ! "Driver" routine for tropopause_find. Identifies the tropopause using several methods
         100             :   ! and populates them into the model state.
         101             :   ! Most methods use climatological tropopause as a backup, and as such is guaranteed to
         102             :   ! find a tropopause in all columns. Others are explicitly single-method and used by
         103             :   ! other parameterizations as-is with NOTFOUND values being intentional.
         104             : !> \section arg_table_tropopause_find_run Argument Table
         105             : !! \htmlinclude tropopause_find_run.html
         106           0 :   subroutine tropopause_find_run(ncol, pver, fillvalue, lat, pint, pmid, t, zi, zm, phis, &
         107           0 :                                  calday, tropp_p_loc, tropp_days, &
         108           0 :                                  tropLev, tropP, tropT, tropZ, & ! Default primary+backup (twmo+climate)
         109           0 :                                  tropLev_twmo, tropP_twmo, tropT_twmo, tropZ_twmo, & ! Primary only (twmo)
         110           0 :                                  tropLev_clim, tropP_clim, tropT_clim, tropZ_clim, & ! Climate-only
         111           0 :                                  tropLev_hybstob, tropP_hybstob, tropT_hybstob, tropZ_hybstob, & ! Hybridstobie + climate backup
         112           0 :                                  tropLev_cpp, tropP_cpp, tropT_cpp, tropZ_cpp, & ! Cold point only
         113           0 :                                  tropLev_chem, tropP_chem, tropT_chem, tropZ_chem, & ! Chemical tropopause only
         114           0 :                                  hstobie_trop, hstobie_linoz, hstobie_tropop, & ! Hybridstobie only for chemistry diagnostics
         115           0 :                                  scheme_name, errmsg, errflg)
         116             :
         117             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         118             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         119             :     real(kind_phys), intent(in)         :: fillvalue     ! Fill value for diagnostic outputs
         120             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         121             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         122             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         123             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         124             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         125             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         126             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         127             :
         128             :     real(kind_phys), intent(in)         :: calday        ! Day of year including fraction of day
         129             :
         130             :     ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
         131             :     real(kind_phys), intent(in)         :: tropp_p_loc(:,:)
         132             :     real(kind_phys), intent(in)         :: tropp_days(:) ! Day-of-year for climo data, 12
         133             :
         134             :     integer,         intent(out)     :: tropLev(:)         ! tropopause level index
         135             :     real(kind_phys), intent(out)     :: tropP(:)           ! tropopause pressure (Pa)
         136             :     real(kind_phys), intent(out)     :: tropT(:)           ! tropopause temperature (K)
         137             :     real(kind_phys), intent(out)     :: tropZ(:)           ! tropopause height (m)
         138             :
         139             :     integer,         intent(out)     :: tropLev_twmo(:)    ! lapse-rate tropopause level index
         140             :     real(kind_phys), intent(out)     :: tropP_twmo(:)      ! lapse-rate tropopause pressure (Pa)
         141             :     real(kind_phys), intent(out)     :: tropT_twmo(:)      ! lapse-rate tropopause temperature (K)
         142             :     real(kind_phys), intent(out)     :: tropZ_twmo(:)      ! lapse-rate tropopause height (m)
         143             :
         144             :     integer,         intent(out)     :: tropLev_clim(:)    ! climatology-backed tropopause level index
         145             :     real(kind_phys), intent(out)     :: tropP_clim(:)      ! climatology-backed tropopause pressure (Pa)
         146             :     real(kind_phys), intent(out)     :: tropT_clim(:)      ! climatology-backed tropopause temperature (K)
         147             :     real(kind_phys), intent(out)     :: tropZ_clim(:)      ! climatology-backed tropopause height (m)
         148             :
         149             :     integer,         intent(out)     :: tropLev_hybstob(:) ! hybridstobie climatology-backed tropopause level index
         150             :     real(kind_phys), intent(out)     :: tropP_hybstob(:)   ! hybridstobie climatology-backed tropopause pressure (Pa)
         151             :     real(kind_phys), intent(out)     :: tropT_hybstob(:)   ! hybridstobie climatology-backed tropopause temperature (K)
         152             :     real(kind_phys), intent(out)     :: tropZ_hybstob(:)   ! hybridstobie climatology-backed tropopause height (m)
         153             :
         154             :     integer,         intent(out)     :: tropLev_cpp(:)     ! cold point tropopause level index
         155             :     real(kind_phys), intent(out)     :: tropP_cpp(:)       ! cold point tropopause pressure (Pa)
         156             :     real(kind_phys), intent(out)     :: tropT_cpp(:)       ! cold point tropopause temperature (K)
         157             :     real(kind_phys), intent(out)     :: tropZ_cpp(:)       ! cold point tropopause height (m)
         158             :
         159             :     integer,         intent(out)     :: tropLev_chem(:)    ! chemical tropopause level index
         160             :     real(kind_phys), intent(out)     :: tropP_chem(:)      ! chemical tropopause pressure (Pa)
         161             :     real(kind_phys), intent(out)     :: tropT_chem(:)      ! chemical tropopause temperature (K)
         162             :     real(kind_phys), intent(out)     :: tropZ_chem(:)      ! chemical tropopause height (m)
         163             :
         164             :     ! Optional output arguments for hybridstobie with chemistry
         165             :     real(kind_phys), intent(out)   :: hstobie_trop(:,:)   ! Lowest level with strat. chem
         166             :     real(kind_phys), intent(out)   :: hstobie_linoz(:,:)  ! Lowest possible Linoz level
         167             :     real(kind_phys), intent(out)   :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
         168             :
         169             :     character(len=64),  intent(out) :: scheme_name
         170             :     character(len=512), intent(out) :: errmsg
         171             :     integer,            intent(out) :: errflg
         172             :
         173           0 :     scheme_name = 'tropopause_find'
         174           0 :     errmsg = ' '
         175           0 :     errflg = 0
         176             :
         177             :     ! Obtain the primary output, which is TWMO + climate
         178             :     call tropopause_findWithBackup( &
         179             :          ncol           = ncol, &
         180             :          pver           = pver, &
         181             :          fillvalue      = fillvalue, &
         182             :          lat            = lat, &
         183             :          pint           = pint, &
         184             :          pmid           = pmid, &
         185             :          t              = t, &
         186             :          zi             = zi, &
         187             :          zm             = zm, &
         188             :          phis           = phis, &
         189             :          calday         = calday, &
         190             :          tropp_p_loc    = tropp_p_loc, &
         191             :          tropp_days     = tropp_days, &
         192             :          tropLev        = tropLev, &
         193             :          tropP          = tropP, &
         194             :          tropT          = tropT, &
         195             :          tropZ          = tropZ, &
         196             :          primary        = default_primary, &
         197             :          backup         = default_backup, &
         198             :          errmsg         = errmsg, &
         199             :          errflg         = errflg &
         200           0 :     )
         201             :
         202             :     ! Any other intended outputs
         203             :     ! Primary (TWMO) only
         204             :     call tropopause_findWithBackup( &
         205             :          ncol           = ncol, &
         206             :          pver           = pver, &
         207             :          fillvalue      = fillvalue, &
         208             :          lat            = lat, &
         209             :          pint           = pint, &
         210             :          pmid           = pmid, &
         211             :          t              = t, &
         212             :          zi             = zi, &
         213             :          zm             = zm, &
         214             :          phis           = phis, &
         215             :          calday         = calday, &
         216             :          tropp_p_loc    = tropp_p_loc, &
         217             :          tropp_days     = tropp_days, &
         218             :          tropLev        = tropLev_twmo, &
         219             :          tropP          = tropP_twmo, &
         220             :          tropT          = tropT_twmo, &
         221             :          tropZ          = tropZ_twmo, &
         222             :          primary        = default_primary, &
         223             :          backup         = TROP_ALG_NONE, &
         224             :          errmsg         = errmsg, &
         225             :          errflg         = errflg &
         226           0 :     )
         227             :
         228             :     ! Climatology only
         229             :     call tropopause_findWithBackup( &
         230             :          ncol           = ncol, &
         231             :          pver           = pver, &
         232             :          fillvalue      = fillvalue, &
         233             :          lat            = lat, &
         234             :          pint           = pint, &
         235             :          pmid           = pmid, &
         236             :          t              = t, &
         237             :          zi             = zi, &
         238             :          zm             = zm, &
         239             :          phis           = phis, &
         240             :          calday         = calday, &
         241             :          tropp_p_loc    = tropp_p_loc, &
         242             :          tropp_days     = tropp_days, &
         243             :          tropLev        = tropLev_clim, &
         244             :          tropP          = tropP_clim, &
         245             :          tropT          = tropT_clim, &
         246             :          tropZ          = tropZ_clim, &
         247             :          primary        = TROP_ALG_CLIMATE, &
         248             :          backup         = TROP_ALG_NONE, &
         249             :          errmsg         = errmsg, &
         250             :          errflg         = errflg &
         251           0 :     )
         252             :
         253             :     ! Cold point (CPP) only
         254             :     call tropopause_findWithBackup( &
         255             :          ncol           = ncol, &
         256             :          pver           = pver, &
         257             :          fillvalue      = fillvalue, &
         258             :          lat            = lat, &
         259             :          pint           = pint, &
         260             :          pmid           = pmid, &
         261             :          t              = t, &
         262             :          zi             = zi, &
         263             :          zm             = zm, &
         264             :          phis           = phis, &
         265             :          calday         = calday, &
         266             :          tropp_p_loc    = tropp_p_loc, &
         267             :          tropp_days     = tropp_days, &
         268             :          tropLev        = tropLev_cpp, &
         269             :          tropP          = tropP_cpp, &
         270             :          tropT          = tropT_cpp, &
         271             :          tropZ          = tropZ_cpp, &
         272             :          primary        = TROP_ALG_CPP, &
         273             :          backup         = TROP_ALG_NONE, &
         274             :          errmsg         = errmsg, &
         275             :          errflg         = errflg &
         276           0 :     )
         277             :
         278             :     ! Hybridstobie with climatology-backed
         279             :     call tropopause_findWithBackup( &
         280             :          ncol           = ncol, &
         281             :          pver           = pver, &
         282             :          fillvalue      = fillvalue, &
         283             :          lat            = lat, &
         284             :          pint           = pint, &
         285             :          pmid           = pmid, &
         286             :          t              = t, &
         287             :          zi             = zi, &
         288             :          zm             = zm, &
         289             :          phis           = phis, &
         290             :          calday         = calday, &
         291             :          tropp_p_loc    = tropp_p_loc, &
         292             :          tropp_days     = tropp_days, &
         293             :          tropLev        = tropLev_hybstob, &
         294             :          tropP          = tropP_hybstob, &
         295             :          tropT          = tropT_hybstob, &
         296             :          tropZ          = tropZ_hybstob, &
         297             :          primary        = TROP_ALG_HYBSTOB, &
         298             :          backup         = TROP_ALG_CLIMATE, &
         299             :          hstobie_trop   = hstobie_trop, &    ! Only used if TROP_ALG_HYBSTOB
         300             :          hstobie_linoz  = hstobie_linoz, &   ! Only used if TROP_ALG_HYBSTOB
         301             :          hstobie_tropop = hstobie_tropop, &  ! Only used if TROP_ALG_HYBSTOB
         302             :          errmsg         = errmsg, &
         303             :          errflg         = errflg &
         304           0 :     )
         305             :
         306             :     ! Chemical tropopause (used for chemistry)
         307             :     call tropopause_findWithBackup( &
         308             :          ncol           = ncol, &
         309             :          pver           = pver, &
         310             :          fillvalue      = fillvalue, &
         311             :          lat            = lat, &
         312             :          pint           = pint, &
         313             :          pmid           = pmid, &
         314             :          t              = t, &
         315             :          zi             = zi, &
         316             :          zm             = zm, &
         317             :          phis           = phis, &
         318             :          calday         = calday, &
         319             :          tropp_p_loc    = tropp_p_loc, &
         320             :          tropp_days     = tropp_days, &
         321             :          tropLev        = tropLev_chem, &
         322             :          tropP          = tropP_chem, &
         323             :          tropT          = tropT_chem, &
         324             :          tropZ          = tropZ_chem, &
         325             :          primary        = TROP_ALG_CHEMTROP, &
         326             :          backup         = TROP_ALG_CLIMATE, &
         327             :          errmsg         = errmsg, &
         328             :          errflg         = errflg &
         329           0 :     )
         330             :
         331           0 :   end subroutine tropopause_find_run
         332             :
         333             :   ! Searches all the columns and attempts to identify the tropopause.
         334             :   ! Two routines can be specifed, a primary routine which is tried first and a
         335             :   ! backup routine which will be tried only if the first routine fails. If the
         336             :   ! tropopause can not be identified by either routine, then a NOTFOUND is returned
         337             :   ! for the tropopause level, temperature and pressure.
         338    15644088 :   subroutine tropopause_findWithBackup(ncol, pver, fillvalue, lat, pint, pmid, t, zi, zm, phis, &
         339    15644088 :                                        calday, tropp_p_loc, tropp_days, &
         340    15644088 :                                        tropLev, tropP, tropT, tropZ, &
         341    15644088 :                                        hstobie_trop, hstobie_linoz, hstobie_tropop, &
         342             :                                        primary, backup, &
         343    15644088 :                                        errmsg, errflg)
         344             :
         345             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         346             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         347             :     real(kind_phys), intent(in)         :: fillvalue     ! Fill value for diagnostic outputs
         348             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         349             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         350             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         351             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         352             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         353             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         354             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         355             :
         356             :     real(kind_phys), intent(in)         :: calday        ! Day of year including fraction of day
         357             :
         358             :     ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
         359             :     real(kind_phys), intent(in)         :: tropp_p_loc(:,:)
         360             :     real(kind_phys), intent(in)         :: tropp_days(:) ! Day-of-year for climo data, 12
         361             :
         362             :     integer,                   intent(out)     :: tropLev(:)          ! tropopause level index
         363             :     real(kind_phys), optional, intent(out)     :: tropP(:)            ! tropopause pressure (Pa)
         364             :     real(kind_phys), optional, intent(out)     :: tropT(:)            ! tropopause temperature (K)
         365             :     real(kind_phys), optional, intent(out)     :: tropZ(:)            ! tropopause height (m)
         366             :
         367             :     ! Optional output arguments for hybridstobie with chemistry
         368             :     real(kind_phys), optional, intent(inout)   :: hstobie_trop(:,:)   ! Lowest level with strat. chem
         369             :     real(kind_phys), optional, intent(inout)   :: hstobie_linoz(:,:)  ! Lowest possible Linoz level
         370             :     real(kind_phys), optional, intent(inout)   :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
         371             :
         372             :     ! primary and backup are no longer optional arguments for CCPP-compliance.
         373             :     ! specify defaults when calling (TWMO, CLIMO)
         374             :     integer, intent(in)       :: primary                   ! primary detection algorithm
         375             :     integer, intent(in)       :: backup                    ! backup detection algorithm
         376             :
         377             :     character(len=512), intent(out) :: errmsg
         378             :     integer,            intent(out) :: errflg
         379             :
         380    15644088 :     errmsg = ' '
         381    15644088 :     errflg = 0
         382             :
         383             :     ! Initialize the results to a missing value, so that the algorithms will
         384             :     ! attempt to find the tropopause for all of them.
         385   261219888 :     tropLev(:) = NOTFOUND
         386    87124536 :     if(present(tropP)) tropP(:) = fillvalue
         387    87124536 :     if(present(tropT)) tropT(:) = fillvalue
         388    87124536 :     if(present(tropZ)) tropZ(:) = fillvalue
         389             :
         390             :     ! Try to find the tropopause using the primary algorithm.
         391    15644088 :     if (primary /= TROP_ALG_NONE) then
         392             :       call tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         393             :                                 calday, tropp_p_loc, tropp_days, &
         394             :                                 primary, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, &
         395             :                                 hstobie_trop=hstobie_trop, hstobie_linoz=hstobie_linoz, hstobie_tropop=hstobie_tropop, & ! only for HYBSTOB
         396    80465040 :                                 errmsg=errmsg, errflg=errflg)
         397             :     end if
         398             :
         399   259869870 :     if ((backup /= TROP_ALG_NONE) .and. any(tropLev(:) == NOTFOUND)) then
         400             :       call tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         401             :                                 calday, tropp_p_loc, tropp_days, &
         402             :                                 backup, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, &
         403       96260 :                                 errmsg=errmsg, errflg=errflg)
         404             :     end if
         405             :
         406    15644088 :   end subroutine tropopause_findWithBackup
         407             :
         408             :   ! Call the appropriate tropopause detection routine based upon the algorithm
         409             :   ! specifed.
         410             :   !
         411             :   ! NOTE: It is assumed that the output fields have been initialized by the
         412             :   ! caller, and only output values set to fillvalue will be detected.
         413    15692090 :   subroutine tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         414    15692090 :                                   calday, tropp_p_loc, tropp_days, &
         415    15692090 :                                   algorithm, tropLev, tropP, tropT, tropZ, &
         416    15692090 :                                   hstobie_trop, hstobie_linoz, hstobie_tropop, &
         417    15692090 :                                   errmsg, errflg)
         418             :
         419             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         420             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         421             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         422             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         423             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         424             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         425             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         426             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         427             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         428             :
         429             :     real(kind_phys), intent(in)         :: calday        ! Day of year including fraction of day
         430             :
         431             :     ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
         432             :     real(kind_phys), intent(in)         :: tropp_p_loc(:,:)
         433             :     real(kind_phys), intent(in)         :: tropp_days(:) ! Day-of-year for climo data, 12
         434             :
         435             :     integer,                   intent(in)      :: algorithm             ! detection algorithm
         436             :     integer,                   intent(inout)   :: tropLev(:)            ! tropopause level index
         437             :     real(kind_phys), optional, intent(inout)   :: tropP(:)              ! tropopause pressure (Pa)
         438             :     real(kind_phys), optional, intent(inout)   :: tropT(:)              ! tropopause temperature (K)
         439             :     real(kind_phys), optional, intent(inout)   :: tropZ(:)              ! tropopause height (m)
         440             :
         441             :     ! Optional output arguments for hybridstobie with chemistry
         442             :     real(kind_phys), optional, intent(inout)   :: hstobie_trop(:,:)   ! Lowest level with strat. chem
         443             :     real(kind_phys), optional, intent(inout)   :: hstobie_linoz(:,:)  ! Lowest possible Linoz level
         444             :     real(kind_phys), optional, intent(inout)   :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
         445             :
         446             :     character(len=512), intent(out) :: errmsg
         447             :     integer,            intent(out) :: errflg
         448             :
         449    15692090 :     errmsg = ' '
         450    15692090 :     errflg = 0
         451             :
         452             :     ! Dispatch the request to the appropriate routine.
         453    15692090 :     select case(algorithm)
         454             :       case(TROP_ALG_ANALYTIC)
         455           0 :         call tropopause_analytic(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
         456             :
         457             :       case(TROP_ALG_CLIMATE)
         458             :         call tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         459             :                                 calday, tropp_p_loc, tropp_days, &
         460       96260 :                                 tropLev, tropP, tropT, tropZ)
         461             :
         462             :       case(TROP_ALG_STOBIE)
         463           0 :         call tropopause_stobie(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
         464             :
         465             :       case(TROP_ALG_HYBSTOB)
         466           0 :         if(present(hstobie_trop) .and. present(hstobie_linoz) .and. present(hstobie_tropop)) then
         467             :           call tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
         468             :                                        tropLev, tropP, tropT, tropZ, &
         469           0 :                                        hstobie_trop, hstobie_linoz, hstobie_tropop)
         470             :         else
         471             :           call tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
         472           0 :                                        tropLev, tropP, tropT, tropZ)
         473             :         endif
         474             :
         475             :       case(TROP_ALG_TWMO)
         476     5962896 :         call tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
         477             :
         478             :       case(TROP_ALG_WMO)
         479           0 :         call tropopause_wmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
         480             :
         481             :       case(TROP_ALG_CPP)
         482     1489176 :         call tropopause_cpp(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
         483             :
         484             :       case(TROP_ALG_CHEMTROP)
         485             :         ! hplin: needs climatological arguments as calling tropopause_climate from within findChemTrop
         486             :         call tropopause_findChemTrop(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         487             :                                      calday, tropp_p_loc, tropp_days, &
         488             :                                      tropLev, tropP, tropT, tropZ, &
         489    41721696 :                                      errmsg, errflg)
         490             :
         491             :       case default
         492           0 :         errflg = 1
         493    15692090 :         write(errmsg,*) 'tropopause: Invalid detection algorithm (',  algorithm, ') specified.'
         494             :     end select
         495             :
         496    15692090 :   end subroutine tropopause_findUsing
         497             :
         498             :   ! This analytic expression closely matches the mean tropopause determined
         499             :   ! by the NCEP reanalysis and has been used by the radiation code.
         500           0 :   subroutine tropopause_analytic(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         501           0 :                                  tropLev, tropP, tropT, tropZ)
         502             :
         503             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         504             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         505             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         506             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         507             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         508             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         509             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         510             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         511             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         512             :     integer,                   intent(inout)   :: tropLev(:)             ! tropopause level index
         513             :     real(kind_phys), optional, intent(inout)   :: tropP(:)               ! tropopause pressure (Pa)
         514             :     real(kind_phys), optional, intent(inout)   :: tropT(:)               ! tropopause temperature (K)
         515             :     real(kind_phys), optional, intent(inout)   :: tropZ(:)               ! tropopause height (m)
         516             :
         517             :     ! Local Variables
         518             :     integer       :: i
         519             :     integer       :: k
         520             :     real(kind_phys)      :: tP                       ! tropopause pressure (Pa)
         521             :
         522             :     ! Iterate over all of the columns.
         523           0 :     do i = 1, ncol
         524             :
         525             :       ! Skip column in which the tropopause has already been found.
         526           0 :       if (tropLev(i) == NOTFOUND) then
         527             :
         528             :         ! Calculate the pressure of the tropopause.
         529           0 :         tP = (25000.0_kind_phys - 15000.0_kind_phys * (cos(lat(i)))**2)
         530             :
         531             :         ! Find the level that contains the tropopause.
         532           0 :         do k = pver, 2, -1
         533           0 :           if (tP >= pint(i, k)) then
         534           0 :             tropLev(i) = k
         535           0 :             exit
         536             :           end if
         537             :         end do
         538             :
         539             :         ! Return the optional outputs
         540           0 :         if (present(tropP)) tropP(i) = tP
         541             :
         542           0 :         if (present(tropT)) then
         543           0 :           tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
         544             :         end if
         545             :
         546           0 :         if (present(tropZ)) then
         547           0 :           tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
         548             :         end if
         549             :       end if
         550             :     end do
         551             :
         552           0 :   end subroutine tropopause_analytic
         553             :
         554             :   ! Determine the tropopause pressure from a climatology,
         555             :   ! interpolated to the current day of year and latitude.
         556      546502 :   subroutine tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         557      546502 :                                 calday, tropp_p_loc, tropp_days, tropLev, tropP, tropT, tropZ)
         558             :
         559             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         560             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         561             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         562             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa), pverp
         563             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         564             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         565             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         566             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         567             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         568             :
         569             :     real(kind_phys), intent(in)         :: calday        ! Day of year including fraction of day
         570             :
         571             :     ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
         572             :     real(kind_phys), intent(in)         :: tropp_p_loc(:,:)
         573             :     real(kind_phys), intent(in)         :: tropp_days(:) ! Day-of-year for climo data, 12
         574             :
         575             :     integer,                   intent(inout)  :: tropLev(:)            ! tropopause level index
         576             :     real(kind_phys), optional, intent(inout)  :: tropP(:)              ! tropopause pressure (Pa)
         577             :     real(kind_phys), optional, intent(inout)  :: tropT(:)              ! tropopause temperature (K)
         578             :     real(kind_phys), optional, intent(inout)  :: tropZ(:)              ! tropopause height (m)
         579             :
         580             :     ! Local Variables
         581             :     integer       :: i
         582             :     integer       :: k
         583             :     integer       :: m
         584             :     real(kind_phys)      :: tP                       ! tropopause pressure (Pa)
         585             :     real(kind_phys)      :: dels
         586             :     integer       :: last
         587             :     integer       :: next
         588             :
         589             :     ! If any columns remain to be indentified, then get the current
         590             :     ! day from the calendar.
         591             :
         592     1454751 :     if (any(tropLev == NOTFOUND)) then
         593             :
         594             :       !--------------------------------------------------------
         595             :       ! ... setup the time interpolation
         596             :       !--------------------------------------------------------
         597      273251 :       if( calday < tropp_days(1) ) then
         598      273251 :         next = 1
         599      273251 :         last = 12
         600      273251 :         dels = (365._kind_phys + calday - tropp_days(12)) / (365._kind_phys + tropp_days(1) - tropp_days(12))
         601           0 :       else if( calday >= tropp_days(12) ) then
         602           0 :         next = 1
         603           0 :         last = 12
         604           0 :         dels = (calday - tropp_days(12)) / (365._kind_phys + tropp_days(1) - tropp_days(12))
         605             :       else
         606           0 :         do m = 11,1,-1
         607           0 :            if( calday >= tropp_days(m) ) then
         608             :               exit
         609             :            end if
         610             :         end do
         611           0 :         last = m
         612           0 :         next = m + 1
         613           0 :         dels = (calday - tropp_days(m)) / (tropp_days(m+1) - tropp_days(m))
         614             :       end if
         615             :
         616      273251 :       dels = max( min( 1._kind_phys,dels ),0._kind_phys )
         617             :
         618             :
         619             :       ! Iterate over all of the columns.
         620     4566917 :       do i = 1, ncol
         621             :
         622             :         ! Skip column in which the tropopause has already been found.
         623     4566917 :         if (tropLev(i) == NOTFOUND) then
         624             :
         625             :         !--------------------------------------------------------
         626             :         ! ... get tropopause level from climatology
         627             :         !--------------------------------------------------------
         628             :           ! Interpolate the tropopause pressure.
         629     1313786 :           tP = tropp_p_loc(i,last) &
         630     2627572 :             + dels * (tropp_p_loc(i,next) - tropp_p_loc(i,last))
         631             :
         632             :           ! Find the associated level.
         633    61055808 :           do k = pver, 2, -1
         634    61055808 :             if (tP >= pint(i, k)) then
         635     1313786 :               tropLev(i) = k
         636     1313786 :               exit
         637             :             end if
         638             :           end do
         639             :
         640             :           ! Return the optional outputs
         641     1313786 :           if (present(tropP)) tropP(i) = tP
         642             :
         643     1313786 :           if (present(tropT)) then
         644      153298 :             tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
         645             :           end if
         646             :
         647     1313786 :           if (present(tropZ)) then
         648      153298 :             tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
         649             :           end if
         650             :         end if
         651             :       end do
         652             :     end if
         653             :
         654      273251 :   end subroutine tropopause_climate
         655             :
         656             :   !-----------------------------------------------------------------------
         657             :   !-----------------------------------------------------------------------
         658           0 :   subroutine tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
         659           0 :                                      tropLev, tropP, tropT, tropZ, &
         660           0 :                                      hstobie_trop, hstobie_linoz, hstobie_tropop)
         661             :
         662             :     !-----------------------------------------------------------------------
         663             :     ! Originally written by Philip Cameron-Smith, LLNL
         664             :     !
         665             :     !   Stobie-Linoz hybrid: the highest altitude of
         666             :     !          a) Stobie algorithm, or
         667             :     !          b) minimum Linoz pressure.
         668             :     !
         669             :     ! NOTE: the ltrop(i) gridbox itself is assumed to be a STRATOSPHERIC gridbox.
         670             :     !-----------------------------------------------------------------------
         671             :     !-----------------------------------------------------------------------
         672             :     !        ... Local variables
         673             :     !-----------------------------------------------------------------------
         674             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         675             :     integer,         intent(in)         :: pver          ! Number of vertical levelserp
         676             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         677             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         678             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         679             :
         680             :     integer,            intent(inout)   :: tropLev(:)             ! tropopause level index
         681             :     real(kind_phys), optional, intent(inout)   :: tropP(:)               ! tropopause pressure (Pa)
         682             :     real(kind_phys), optional, intent(inout)   :: tropT(:)               ! tropopause temperature (K)
         683             :     real(kind_phys), optional, intent(inout)   :: tropZ(:)               ! tropopause height (m)
         684             :
         685             :     ! Optional output arguments for hybridstobie with chemistry
         686             :     real(kind_phys), optional, intent(inout)   :: hstobie_trop(:,:)   ! Lowest level with strat. chem
         687             :     real(kind_phys), optional, intent(inout)   :: hstobie_linoz(:,:)  ! Lowest possible Linoz level
         688             :     real(kind_phys), optional, intent(inout)   :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
         689             :
         690             :     real(kind_phys),parameter  ::  min_Stobie_Pressure= 40.E2_kind_phys !For case 2 & 4.  [Pa]
         691             :     real(kind_phys),parameter  ::  max_Linoz_Pressure =208.E2_kind_phys !For case     4.  [Pa]
         692             :
         693             :     integer      :: i, k
         694             :     real(kind_phys)     :: stobie_min, shybrid_temp      !temporary variable for case 2 & 3.
         695           0 :     integer      :: ltrop_linoz(ncol)            !Lowest possible Linoz vertical level
         696           0 :     integer      :: ltrop_trop(ncol)             !Tropopause level for hybrid case.
         697             :     logical      :: ltrop_linoz_set               !Flag that lowest linoz level already found.
         698           0 :     real(kind_phys)     :: trop_output(ncol,pver)        !For output purposes only.
         699           0 :     real(kind_phys)     :: trop_linoz_output(ncol,pver)  !For output purposes only.
         700           0 :     real(kind_phys)     :: trop_trop_output(ncol,pver)   !For output purposes only.
         701             :
         702           0 :     ltrop_linoz(:) = 1  ! Initialize to default value.
         703           0 :     ltrop_trop(:) = 1   ! Initialize to default value.
         704             :
         705           0 :     LOOP_COL4: do i=1,ncol
         706             :
         707             :        ! Skip column in which the tropopause has already been found.
         708           0 :        not_found: if (tropLev(i) == NOTFOUND) then
         709             :
         710             :           stobie_min = 1.e10_kind_phys    ! An impossibly large number
         711             :           ltrop_linoz_set = .FALSE.
         712           0 :           LOOP_LEV: do k=pver,1,-1
         713           0 :              IF (pmid(i,k) < min_stobie_pressure) cycle
         714           0 :              shybrid_temp = ALPHA * t(i,k) - Log10(pmid(i,k))
         715             :              !PJC_NOTE: the units of pmid won't matter, because it is just an additive offset.
         716           0 :              IF (shybrid_temp<stobie_min) then
         717           0 :                 ltrop_trop(i)=k
         718           0 :                 stobie_min = shybrid_temp
         719             :              ENDIF
         720           0 :              IF (pmid(i,k) < max_Linoz_pressure .AND. .NOT. ltrop_linoz_set) THEN
         721           0 :                 ltrop_linoz(i) = k
         722           0 :                 ltrop_linoz_set = .TRUE.
         723             :              ENDIF
         724             :           enddo LOOP_LEV
         725             :
         726           0 :           tropLev(i) = MIN(ltrop_trop(i),ltrop_linoz(i))
         727             :
         728           0 :           if (present(tropP)) then
         729           0 :              tropP(i) = pmid(i,tropLev(i))
         730             :           endif
         731           0 :           if (present(tropT)) then
         732           0 :              tropT(i) = t(i,tropLev(i))
         733             :           endif
         734           0 :           if (present(tropZ)) then
         735           0 :              tropZ(i) = zm(i,tropLev(i))
         736             :           endif
         737             :
         738             :        endif not_found
         739             :
         740             :     enddo LOOP_COL4
         741             :
         742           0 :     trop_output(:,:)=0._kind_phys
         743           0 :     trop_linoz_output(:,:)=0._kind_phys
         744           0 :     trop_trop_output(:,:)=0._kind_phys
         745           0 :     do i=1,ncol
         746           0 :        if (tropLev(i)>0) then
         747           0 :           trop_output(i,tropLev(i))=1._kind_phys
         748           0 :           trop_linoz_output(i,ltrop_linoz(i))=1._kind_phys
         749           0 :           trop_trop_output(i,ltrop_trop(i))=1._kind_phys
         750             :        endif
         751             :     enddo
         752             :
         753           0 :     if(present(hstobie_trop)) then
         754           0 :       hstobie_trop(:,:) = trop_output(:,:)
         755             :     endif
         756             :
         757           0 :     if(present(hstobie_linoz)) then
         758           0 :       hstobie_linoz(:,:) = trop_linoz_output(:,:)
         759             :     endif
         760             :
         761           0 :     if(present(hstobie_tropop)) then
         762           0 :       hstobie_tropop(:,:) = trop_trop_output(:,:)
         763             :     endif
         764             :
         765           0 :   end subroutine tropopause_hybridstobie
         766             :
         767             :   ! This routine originates with Stobie at NASA Goddard, but does not have a
         768             :   ! known reference. It was supplied by Philip Cameron-Smith of LLNL.
         769             :   !
         770           0 :   subroutine tropopause_stobie(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         771           0 :                                tropLev, tropP, tropT, tropZ)
         772             :
         773             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         774             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         775             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         776             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         777             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         778             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         779             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         780             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         781             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         782             :
         783             :     integer,            intent(inout)   :: tropLev(:)             ! tropopause level index
         784             :     real(kind_phys), optional, intent(inout)   :: tropP(:)               ! tropopause pressure (Pa)
         785             :     real(kind_phys), optional, intent(inout)   :: tropT(:)               ! tropopause temperature (K)
         786             :     real(kind_phys), optional, intent(inout)   :: tropZ(:)               ! tropopause height (m)
         787             :
         788             :     ! Local Variables
         789             :     integer       :: i
         790             :     integer       :: k
         791             :     integer       :: tLev                     ! tropopause level
         792             :     real(kind_phys)      :: tP                       ! tropopause pressure (Pa)
         793           0 :     real(kind_phys)      :: stobie(pver)             ! stobie weighted temperature
         794             :     real(kind_phys)      :: sTrop                    ! stobie value at the tropopause
         795             :
         796             :     ! Iterate over all of the columns.
         797           0 :     do i = 1, ncol
         798             :
         799             :       ! Skip column in which the tropopause has already been found.
         800           0 :       if (tropLev(i) == NOTFOUND) then
         801             :
         802             :         ! Caclulate a pressure weighted temperature.
         803           0 :         stobie(:) = ALPHA * t(i,:) - log10(pmid(i, :))
         804             :
         805             :         ! Search from the bottom up, looking for the first minimum.
         806           0 :         tLev  = -1
         807             :
         808           0 :         do k = pver-1, 1, -1
         809             :
         810           0 :           if (pmid(i, k) <= 4000._kind_phys) then
         811             :             exit
         812             :           end if
         813             :
         814           0 :           if (pmid(i, k) >= 55000._kind_phys) then
         815             :             cycle
         816             :           end if
         817             :
         818           0 :           if ((tLev == -1) .or. (stobie(k) < sTrop)) then
         819           0 :             tLev  = k
         820           0 :             tP    = pmid(i, k)
         821           0 :             sTrop = stobie(k)
         822             :           end if
         823             :         end do
         824             :
         825           0 :         if (tLev /= -1) then
         826           0 :           tropLev(i) = tLev
         827             :
         828             :           ! Return the optional outputs
         829           0 :           if (present(tropP)) tropP(i) = tP
         830             :
         831           0 :           if (present(tropT)) then
         832           0 :             tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
         833             :           end if
         834             :
         835           0 :           if (present(tropZ)) then
         836           0 :             tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
         837             :           end if
         838             :         end if
         839             :       end if
         840             :     end do
         841             :
         842           0 :   end subroutine tropopause_stobie
         843             :
         844             :
         845             :   ! This routine is an implementation of Reichler et al. [2003] done by
         846             :   ! Reichler and downloaded from his web site. Minimal modifications were
         847             :   ! made to have the routine work within the CAM framework (i.e. using
         848             :   ! CAM constants and types).
         849             :   !
         850             :   ! NOTE: I am not a big fan of the goto's and multiple returns in this
         851             :   ! code, but for the moment I have left them to preserve as much of the
         852             :   ! original and presumably well tested code as possible.
         853             :   ! UPDATE: The most "obvious" substitutions have been made to replace
         854             :   ! goto/return statements with cycle/exit. The structure is still
         855             :   ! somewhat tangled.
         856             :   ! UPDATE 2: "gamma" renamed to "gam" in order to avoid confusion
         857             :   ! with the Fortran 2008 intrinsic. "level" argument removed because
         858             :   ! a physics column is not contiguous, so using explicit dimensions
         859             :   ! will cause the data to be needlessly copied.
         860             :   !
         861             :   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         862             :   !
         863             :   ! determination of tropopause height from gridded temperature data
         864             :   !
         865             :   ! Reichler, T., M. Dameris, and R. Sausen (2003),
         866             :   ! Determining the tropopause height from gridded data,
         867             :   ! Geophys. Res. Lett., 30, 2042, doi:10.1029/2003GL018240, 20.
         868             :   !
         869             :   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         870   222199200 :   subroutine twmo(t, p, plimu, pliml, gam, trp)
         871             :
         872             :     real(kind_phys), intent(in), dimension(:)      :: t, p
         873             :     real(kind_phys), intent(in)                    :: plimu, pliml, gam
         874             :     real(kind_phys), intent(out)                   :: trp
         875             :
         876             :     real(kind_phys), parameter                     :: deltaz = 2000.0_kind_phys
         877             :
         878             :     real(kind_phys)                                :: pmk, pm, a, b, tm, dtdp, dtdz
         879             :     real(kind_phys)                                :: ag, bg, ptph
         880             :     real(kind_phys)                                :: pm0, pmk0, dtdz0
         881             :     real(kind_phys)                                :: p2km, asum, aquer
         882             :     real(kind_phys)                                :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
         883             :     integer                                 :: level
         884             :     integer                                 :: icount, jj
         885             :     integer                                 :: j
         886             :
         887             :
         888   222199200 :     trp=-99.0_kind_phys                           ! negative means not valid
         889             :
         890             :     ! initialize start level
         891             :     ! dt/dz
         892   222199200 :     level = size(t)
         893   222199200 :     pmk= .5_kind_phys * (p(level-1)**cnst_kap+p(level)**cnst_kap)
         894   222199200 :     pm = pmk**(1/cnst_kap)
         895   222199200 :     a = (t(level-1)-t(level))/(p(level-1)**cnst_kap-p(level)**cnst_kap)
         896   222199200 :     b = t(level)-(a*p(level)**cnst_kap)
         897   222199200 :     tm = a * pmk + b
         898   222199200 :     dtdp = a * cnst_kap * (pm**cnst_ka1)
         899   222199200 :     dtdz = cnst_faktor*dtdp*pm/tm
         900             :
         901  9130974539 :     main_loop: do j=level-1,2,-1
         902             :       pm0 = pm
         903  9129507504 :       pmk0 = pmk
         904  9129507504 :       dtdz0  = dtdz
         905             :
         906             :       ! dt/dz
         907  9129507504 :       pmk= .5_kind_phys * (p(j-1)**cnst_kap+p(j)**cnst_kap)
         908  9129507504 :       pm = pmk**(1/cnst_kap)
         909  9129507504 :       a = (t(j-1)-t(j))/(p(j-1)**cnst_kap-p(j)**cnst_kap)
         910  9129507504 :       b = t(j)-(a*p(j)**cnst_kap)
         911  9129507504 :       tm = a * pmk + b
         912  9129507504 :       dtdp = a * cnst_kap * (pm**cnst_ka1)
         913  9129507504 :       dtdz = cnst_faktor*dtdp*pm/tm
         914             :       ! dt/dz valid?
         915  9129507504 :       if (dtdz.le.gam) cycle main_loop    ! no, dt/dz < -2 K/km
         916   998046624 :       if (pm.gt.plimu)   cycle main_loop    ! no, too low
         917             :
         918             :       ! dtdz is valid, calculate tropopause pressure
         919   345584039 :       if (dtdz0.lt.gam) then
         920   271777273 :         ag = (dtdz-dtdz0) / (pmk-pmk0)
         921   271777273 :         bg = dtdz0 - (ag * pmk0)
         922   271777273 :         ptph = exp(log((gam-bg)/ag)/cnst_kap)
         923             :       else
         924             :         ptph = pm
         925             :       endif
         926             :
         927   345584039 :       if (ptph.lt.pliml) cycle main_loop
         928   292328389 :       if (ptph.gt.plimu) cycle main_loop
         929             :
         930             :       ! 2nd test: dtdz above 2 km must not exceed gam
         931   292135934 :       p2km = ptph + deltaz*(pm/tm)*cnst_faktor     ! p at ptph + 2km
         932   292135934 :       asum = 0.0_kind_phys                                ! dtdz above
         933   292135934 :       icount = 0                                   ! number of levels above
         934             :
         935             :       ! test until apm < p2km
         936  1639332767 :       in_loop: do jj=j,2,-1
         937             :
         938  1639332767 :         pmk2 = .5_kind_phys * (p(jj-1)**cnst_kap+p(jj)**cnst_kap) ! p mean ^kappa
         939  1639332767 :         pm2 = pmk2**(1/cnst_kap)                           ! p mean
         940  1639332767 :         if(pm2.gt.ptph) cycle in_loop            ! doesn't happen
         941  1639332767 :         if(pm2.lt.p2km) exit in_loop             ! ptropo is valid
         942             :
         943  1418600602 :         a2 = (t(jj-1)-t(jj))                     ! a
         944  1418600602 :         a2 = a2/(p(jj-1)**cnst_kap-p(jj)**cnst_kap)
         945  1418600602 :         b2 = t(jj)-(a2*p(jj)**cnst_kap)          ! b
         946  1418600602 :         tm2 = a2 * pmk2 + b2                     ! T mean
         947  1418600602 :         dtdp2 = a2 * cnst_kap * (pm2**(cnst_kap-1))  ! dt/dp
         948  1418600602 :         dtdz2 = cnst_faktor*dtdp2*pm2/tm2
         949  1418600602 :         asum = asum+dtdz2
         950  1418600602 :         icount = icount+1
         951  1418600602 :         aquer = asum/float(icount)               ! dt/dz mean
         952             :
         953             :         ! discard ptropo ?
         954  1639332767 :         if (aquer.le.gam) cycle main_loop      ! dt/dz above < gam
         955             :
         956             :       enddo in_loop  ! test next level
         957             :
         958   220732165 :       trp = ptph
         959  9130974539 :       exit main_loop
         960             :     enddo main_loop
         961             :
         962   222199200 :   end subroutine twmo
         963             :
         964             :
         965             :   ! This routine uses an implementation of Reichler et al. [2003] done by
         966             :   ! Reichler and downloaded from his web site. This is similar to the WMO
         967             :   ! routines, but is designed for GCMs with a coarse vertical grid.
         968    28309824 :   subroutine tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
         969    28309824 :                              tropLev, tropP, tropT, tropZ)
         970             :
         971             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
         972             :     integer,         intent(in)         :: pver          ! Number of vertical levels
         973             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
         974             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
         975             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
         976             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
         977             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
         978             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
         979             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
         980             :     integer,            intent(inout)  :: tropLev(:)            ! tropopause level index
         981             :     real(kind_phys), optional, intent(inout)  :: tropP(:)              ! tropopause pressure (Pa)
         982             :     real(kind_phys), optional, intent(inout)  :: tropT(:)              ! tropopause temperature (K)
         983             :     real(kind_phys), optional, intent(inout)  :: tropZ(:)              ! tropopause height (m)
         984             :
         985             :     ! Local Variables
         986             :     real(kind_phys), parameter     :: gam    = -0.002_kind_phys         ! K/m
         987             :     real(kind_phys), parameter     :: plimu    = 45000._kind_phys         ! Pa
         988             :     real(kind_phys), parameter     :: pliml    = 7500._kind_phys          ! Pa
         989             :
         990             :     integer                 :: i
         991             :     integer                 :: k
         992             :     real(kind_phys)                :: tP                       ! tropopause pressure (Pa)
         993             :
         994             :     ! Iterate over all of the columns.
         995   236354112 :     do i = 1, ncol
         996             :
         997             :       ! Skip column in which the tropopause has already been found.
         998   236354112 :       if (tropLev(i) == NOTFOUND) then
         999             :
        1000             :         ! Use the routine from Reichler.
        1001   222199200 :         call twmo(t(i, :), pmid(i, :), plimu, pliml, gam, tP)
        1002             :
        1003             :         ! if successful, store of the results and find the level and temperature.
        1004   222199200 :         if (tP > 0) then
        1005             :
        1006             :           ! Find the associated level.
        1007  9214301784 :           do k = pver, 2, -1
        1008  9214301784 :             if (tP >= pint(i, k)) then
        1009   220732165 :               tropLev(i) = k
        1010   220732165 :               exit
        1011             :             end if
        1012             :           end do
        1013             :
        1014             :           ! Return the optional outputs
        1015   220732165 :           if (present(tropP)) tropP(i) = tP
        1016             :
        1017   220732165 :           if (present(tropT)) then
        1018    46446604 :             tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
        1019             :           end if
        1020             :
        1021   220732165 :           if (present(tropZ)) then
        1022    46446604 :             tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
        1023             :           end if
        1024             :         end if
        1025             :       end if
        1026             :     end do
        1027             :
        1028    14154912 :   end subroutine tropopause_twmo
        1029             :
        1030             :   ! This routine implements the WMO definition of the tropopause (WMO, 1957; Seidel and Randel, 2006).
        1031             :   ! Seidel, D. J., and W. J. Randel (2006),
        1032             :   ! Variability and trends in the global tropopause estimated from radiosonde data,
        1033             :   ! J. Geophys. Res., 111, D21101, doi:10.1029/2006JD007363.
        1034             :   !
        1035             :   ! This requires that the lapse rate be less than 2 K/km for an altitude range
        1036             :   ! of 2 km. The search starts at the surface and stops the first time this
        1037             :   ! criteria is met.
        1038           0 :   subroutine tropopause_wmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
        1039           0 :                             tropLev, tropP, tropT, tropZ)
        1040             :
        1041             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
        1042             :     integer,         intent(in)         :: pver          ! Number of vertical levels
        1043             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
        1044             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
        1045             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1046             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
        1047             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
        1048             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
        1049             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
        1050             :
        1051             :     integer,            intent(inout)  :: tropLev(:)            ! tropopause level index
        1052             :     real(kind_phys), optional, intent(inout)  :: tropP(:)              ! tropopause pressure (Pa)
        1053             :     real(kind_phys), optional, intent(inout)  :: tropT(:)              ! tropopause temperature (K)
        1054             :     real(kind_phys), optional, intent(inout)  :: tropZ(:)              ! tropopause height (m)
        1055             :
        1056             :     ! Local Variables
        1057             :     real(kind_phys), parameter    :: ztrop_low   = 5000._kind_phys        ! lowest tropopause level allowed (m)
        1058             :     real(kind_phys), parameter    :: ztrop_high  = 20000._kind_phys       ! highest tropopause level allowed (m)
        1059             :     real(kind_phys), parameter    :: max_dtdz    = 0.002_kind_phys        ! max dt/dz for tropopause level (K/m)
        1060             :     real(kind_phys), parameter    :: min_trop_dz = 2000._kind_phys        ! min tropopause thickness (m)
        1061             :
        1062             :     integer                 :: i
        1063             :     integer                 :: k
        1064             :     integer                 :: k2
        1065             :     real(kind_phys)                :: tP                           ! tropopause pressure (Pa)
        1066             :     real(kind_phys)                :: dt
        1067             :
        1068             :     ! Iterate over all of the columns.
        1069           0 :     do i = 1, ncol
        1070             :
        1071             :       ! Skip column in which the tropopause has already been found.
        1072           0 :       if (tropLev(i) == NOTFOUND) then
        1073             :
        1074           0 :         kloop: do k = pver-1, 2, -1
        1075             :
        1076             :           ! Skip levels below the minimum and stop if nothing is found
        1077             :           ! before the maximum.
        1078           0 :           if (zm(i, k) < ztrop_low) then
        1079             :             cycle kloop
        1080           0 :           else if (zm(i, k) > ztrop_high) then
        1081             :             exit kloop
        1082             :           end if
        1083             :
        1084             :           ! Compare the actual lapse rate to the threshold
        1085           0 :           dt = t(i, k) - t(i, k-1)
        1086             :
        1087           0 :           if (dt <= (max_dtdz * (zm(i, k-1) - zm(i, k)))) then
        1088             :
        1089             :             ! Make sure that the lapse rate stays below the threshold for the
        1090             :             ! specified range.
        1091           0 :             k2loop: do k2 = k-1, 2, -1
        1092           0 :               if ((zm(i, k2) - zm(i, k)) >= min_trop_dz) then
        1093           0 :                 tP = pmid(i, k)
        1094           0 :                 tropLev(i) = k
        1095           0 :                 exit k2loop
        1096             :               end if
        1097             :
        1098           0 :               dt = t(i, k) - t(i, k2)
        1099           0 :               if (dt > (max_dtdz * (zm(i, k2) - zm(i, k)))) then
        1100             :                 exit k2loop
        1101             :               end if
        1102             :            end do k2loop
        1103             :
        1104           0 :            if (tropLev(i) == NOTFOUND) then
        1105             :               cycle kloop
        1106             :            else
        1107             :
        1108             :               ! Return the optional outputs
        1109           0 :               if (present(tropP)) tropP(i) = tP
        1110             :
        1111           0 :               if (present(tropT)) then
        1112           0 :                 tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
        1113             :               end if
        1114             :
        1115           0 :               if (present(tropZ)) then
        1116           0 :                 tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
        1117             :               end if
        1118             :
        1119             :               exit kloop
        1120             :             end if
        1121             :           end if
        1122             :         end do kloop
        1123             :       end if
        1124             :     end do
        1125             :
        1126           0 :   end subroutine tropopause_wmo
        1127             :
        1128             :
        1129             :   ! This routine searches for the cold point tropopause, and uses a parabolic
        1130             :   ! fit of the coldest point and two adjacent points to interpolate the cold point
        1131             :   ! between model levels.
        1132     2978352 :   subroutine tropopause_cpp(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
        1133     2978352 :                             tropLev, tropP, tropT, tropZ)
        1134             :
        1135             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
        1136             :     integer,         intent(in)         :: pver          ! Number of vertical levels
        1137             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
        1138             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
        1139             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1140             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
        1141             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
        1142             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
        1143             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
        1144             :     integer,            intent(inout)  :: tropLev(:)            ! tropopause level index
        1145             :     real(kind_phys), optional, intent(inout)  :: tropP(:)              ! tropopause pressure (Pa)
        1146             :     real(kind_phys), optional, intent(inout)  :: tropT(:)              ! tropopause temperature (K)
        1147             :     real(kind_phys), optional, intent(inout)  :: tropZ(:)              ! tropopause height (m)
        1148             :
        1149             :     ! Local Variables
        1150             :     real(kind_phys), parameter    :: ztrop_low   = 5000._kind_phys        ! lowest tropopause level allowed (m)
        1151             :     real(kind_phys), parameter    :: ztrop_high  = 25000._kind_phys       ! highest tropopause level allowed (m)
        1152             :
        1153             :     integer                 :: i
        1154             :     integer                 :: k, firstk, lastk
        1155             :     integer                 :: k2
        1156             :     real(kind_phys)                :: tZ                           ! tropopause height (m)
        1157             :     real(kind_phys)                :: tmin
        1158             :     real(kind_phys)                :: f0, f1, f2
        1159             :     real(kind_phys)                :: x0, x1, x2
        1160             :     real(kind_phys)                :: c0, c1, c2
        1161             :     real(kind_phys)                :: a, b, c
        1162             :
        1163             :     ! Iterate over all of the columns.
        1164    24865776 :     do i = 1, ncol
        1165             :
        1166    23376600 :       firstk = 0
        1167    23376600 :       lastk  = pver+1
        1168             :
        1169             :       ! Skip column in which the tropopause has already been found.
        1170    24865776 :       if (tropLev(i) == NOTFOUND) then
        1171    23376600 :         tmin = 1e6_kind_phys
        1172             :
        1173  1573256627 :         kloop: do k = pver-1, 2, -1
        1174             :
        1175             :           ! Skip levels below the minimum and stop if nothing is found
        1176             :           ! before the maximum.
        1177  1573256627 :           if (zm(i, k) < ztrop_low) then
        1178             :             firstk = k
        1179             :             cycle kloop
        1180  1048917338 :           else if (zm(i, k) > ztrop_high) then
        1181             :             lastk = k
        1182             :             exit kloop
        1183             :           end if
        1184             :
        1185             :           ! Find the coldest point
        1186  1025540738 :           if (t(i, k) < tmin) then
        1187   523032813 :             tropLev(i) = k
        1188   523032813 :             tmin = t(i,k)
        1189             :           end if
        1190             :         end do kloop
        1191             :
        1192             :
        1193             :         ! If the minimum is at the edge of the search range, then don't
        1194             :         ! consider this to be a minima
        1195    23376600 :         if ((tropLev(i) >= (firstk-1)) .or. (tropLev(i) <= (lastk+1))) then
        1196      268471 :           tropLev(i) = NOTFOUND
        1197             :         else
        1198             :
        1199             :           ! If returning P, Z, or T, then do a parabolic fit using the
        1200             :           ! cold point and it its 2 surrounding points to interpolate
        1201             :           ! between model levels.
        1202    23108129 :           if (present(tropP) .or. present(tropZ) .or. present(tropT)) then
        1203    23108129 :             f0 = t(i, tropLev(i)-1)
        1204    23108129 :             f1 = t(i, tropLev(i))
        1205    23108129 :             f2 = t(i, tropLev(i)+1)
        1206             :
        1207    23108129 :             x0 = zm(i, tropLev(i)-1)
        1208    23108129 :             x1 = zm(i, tropLev(i))
        1209    23108129 :             x2 = zm(i, tropLev(i)+1)
        1210             :
        1211    23108129 :             c0 = (x0-x1)*(x0-x2)
        1212    23108129 :             c1 = (x1-x0)*(x1-x2)
        1213    23108129 :             c2 = (x2-x0)*(x2-x1)
        1214             :
        1215             :             ! Determine the quadratic coefficients of:
        1216             :             !   T = a * z^2 - b*z + c
        1217    23108129 :             a = (f0/c0 + f1/c1 + f2/c2)
        1218    23108129 :             b = (f0/c0*(x1+x2) + f1/c1*(x0+x2) + f2/c2*(x0+x1))
        1219    23108129 :             c = f0/c0*x1*x2 + f1/c1*x0*x2 + f2/c2*x0*x1
        1220             :
        1221             :             ! Find the altitude of the minimum temperature
        1222    23108129 :             tZ = 0.5_kind_phys * b / a
        1223             :
        1224             :             ! The fit should be between the upper and lower points,
        1225             :             ! so skip the point if the fit fails.
        1226    23108129 :             if ((tZ >= x0) .or. (tZ <= x2)) then
        1227           0 :               tropLev(i) = NOTFOUND
        1228             :             else
        1229             :
        1230             :               ! Return the optional outputs
        1231    23108129 :               if (present(tropP)) then
        1232    23108129 :                 tropP(i) = tropopause_interpolateP(pver, pmid, zm, i, tropLev(i), tZ)
        1233             :               end if
        1234             :
        1235    23108129 :               if (present(tropT)) then
        1236    23108129 :                 tropT(i) = a * tZ*tZ - b*tZ + c
        1237             :               end if
        1238             :
        1239    23108129 :               if (present(tropZ)) then
        1240    23108129 :                 tropZ(i) = tZ
        1241             :               end if
        1242             :             end if
        1243             :           end if
        1244             :         end if
        1245             :       end if
        1246             :     end do
        1247             :
        1248     1489176 :   end subroutine tropopause_cpp
        1249             :
        1250             :   ! Searches all the columns and attempts to identify the "chemical"
        1251             :   ! tropopause. This is the lapse rate tropopause, backed up by the climatology
        1252             :   ! if the lapse rate fails to find the tropopause at pressures higher than a certain
        1253             :   ! threshold. This pressure threshold depends on latitude. Between 50S and 50N,
        1254             :   ! the climatology is used if the lapse rate tropopause is not found at P > 75 hPa.
        1255             :   ! At high latitude (poleward of 50), the threshold is increased to 125 hPa to
        1256             :   ! eliminate false events that are sometimes detected in the cold polar stratosphere.
        1257             :   !
        1258             :   ! NOTE: This routine was adapted from code in chemistry.F90 and mo_gasphase_chemdr.F90.
        1259             :   ! During the CCPP-ization, findChemTrop is now called from tropopause_find_run using method CHEMTROP
        1260             :   ! and now also returns the standard tropLev, tropP, tropT, tropZ outputs (optional).
        1261             :   ! The "backup" option is dropped as it is not used anywhere in current CAM.
        1262    10430424 :   subroutine tropopause_findChemTrop(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
        1263    10430424 :                                      calday, tropp_p_loc, tropp_days, &
        1264    10430424 :                                      tropLev, tropP, tropT, tropZ, errmsg, errflg)
        1265             :
        1266             :     integer,         intent(in)         :: ncol          ! Number of atmospheric columns
        1267             :     integer,         intent(in)         :: pver          ! Number of vertical levels
        1268             :     real(kind_phys), intent(in)         :: lat(:)        ! Latitudes (radians)
        1269             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
        1270             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1271             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
        1272             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
        1273             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
        1274             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
        1275             :
        1276             :     real(kind_phys), intent(in)         :: calday        ! Day of year including fraction of day
        1277             :
        1278             :     ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
        1279             :     real(kind_phys), intent(in)         :: tropp_p_loc(:,:)
        1280             :     real(kind_phys), intent(in)         :: tropp_days(:) ! Day-of-year for climo data, 12
        1281             :
        1282             :     integer,                   intent(out)     :: tropLev(:)            ! tropopause level index
        1283             :     real(kind_phys), optional, intent(inout)   :: tropP(:)              ! tropopause pressure (Pa)
        1284             :     real(kind_phys), optional, intent(inout)   :: tropT(:)              ! tropopause temperature (K)
        1285             :     real(kind_phys), optional, intent(inout)   :: tropZ(:)              ! tropopause height (m)
        1286             :
        1287             :     character(len=512), intent(out) :: errmsg
        1288             :     integer,            intent(out) :: errflg
        1289             :
        1290             :     ! Local Variable
        1291    10430424 :     real(kind_phys)     :: dlats(ncol)
        1292             :     integer             :: i
        1293             :
        1294    10430424 :     errmsg = ' '
        1295    10430424 :     errflg = 0
        1296             :
        1297             :     ! First use the lapse rate tropopause.
        1298    41721696 :     call tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
        1299             :
        1300             :     ! Now check high latitudes (poleward of 50) and set the level to the
        1301             :     ! climatology if the level was not found or is at P <= 125 hPa.
        1302   174163824 :     dlats(:ncol) = lat(:ncol) * cnst_rad2deg ! convert to degrees
        1303             :
        1304   174163824 :     do i = 1, ncol
        1305   174163824 :       if (abs(dlats(i)) > 50._kind_phys) then
        1306    36735576 :         if (tropLev(i) .ne. NOTFOUND) then
        1307    36735576 :           if (pmid(i, tropLev(i)) <= 12500._kind_phys) then
        1308          49 :             tropLev(i) = NOTFOUND
        1309             :           end if
        1310             :         end if
        1311             :       end if
        1312             :     end do
        1313             :
        1314             :     ! Now use the climatology backup
        1315   171598308 :     if (any(tropLev(:) == NOTFOUND)) then
        1316             :       call tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
        1317             :                               calday, tropp_p_loc, tropp_days, &
        1318      900996 :                               tropLev, tropP, tropT, tropZ)
        1319             :     end if
        1320             :
        1321    10430424 :   end subroutine tropopause_findChemTrop
        1322             :
        1323             :   ! This routine interpolates the pressures in the physics state to
        1324             :   ! find the pressure at the specified tropopause altitude.
        1325    23108129 :   function tropopause_interpolateP(pver, pmid, zm, icol, tropLev, tropZ)
        1326             :
        1327             :     implicit none
        1328             :
        1329             :     integer,         intent(in)         :: pver          ! Number of vertical levels
        1330             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1331             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
        1332             :     integer, intent(in)                 :: icol               ! column being processed
        1333             :     integer, intent(in)                 :: tropLev            ! tropopause level index
        1334             :     real(kind_phys), optional, intent(in)      :: tropZ              ! tropopause pressure (m)
        1335             :     real(kind_phys)                            :: tropopause_interpolateP
        1336             :
        1337             :     ! Local Variables
        1338             :     real(kind_phys)   :: tropP              ! tropopause pressure (Pa)
        1339             :     real(kind_phys)   :: dlogPdZ            ! dlog(p)/dZ
        1340             :
        1341             :     ! Interpolate the temperature linearly against log(P)
        1342             :
        1343             :     ! Is the tropopause at the midpoint?
        1344    23108129 :     if (tropZ == zm(icol, tropLev)) then
        1345           0 :       tropP = pmid(icol, tropLev)
        1346             :
        1347    23108129 :     else if (tropZ > zm(icol, tropLev)) then
        1348             :
        1349             :       ! It is above the midpoint? Make sure we aren't at the top.
        1350    11123538 :       if (tropLev > 1) then
        1351    22247076 :         dlogPdZ = (log(pmid(icol, tropLev)) - log(pmid(icol, tropLev - 1))) / &
        1352    22247076 :           (zm(icol, tropLev) - zm(icol, tropLev - 1))
        1353    11123538 :         tropP = pmid(icol, tropLev) + exp((tropZ - zm(icol, tropLev)) * dlogPdZ)
        1354             :       end if
        1355             :     else
        1356             :
        1357             :       ! It is below the midpoint. Make sure we aren't at the bottom.
        1358    11984591 :       if (tropLev < pver) then
        1359    11984591 :         dlogPdZ =  (log(pmid(icol, tropLev + 1)) - log(pmid(icol, tropLev))) / &
        1360    23969182 :           (zm(icol, tropLev + 1) - zm(icol, tropLev))
        1361    11984591 :         tropP = pmid(icol, tropLev) + exp((tropZ - zm(icol, tropLev)) * dlogPdZ)
        1362             :       end if
        1363             :     end if
        1364             :
        1365    23108129 :     tropopause_interpolateP = tropP
        1366    23108129 :   end function tropopause_interpolateP
        1367             :
        1368             :
        1369             :   ! This routine interpolates the temperatures in the physics state to
        1370             :   ! find the temperature at the specified tropopause pressure.
        1371    46599902 :   function tropopause_interpolateT(pver, pmid, t, icol, tropLev, tropP)
        1372             :
        1373             :     implicit none
        1374             :
        1375             :     integer,         intent(in)         :: pver          ! Number of vertical levels
        1376             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1377             :     real(kind_phys), intent(in)         :: t(:,:)        ! Temperature (K)
        1378             :     integer, intent(in)                 :: icol               ! column being processed
        1379             :     integer, intent(in)                 :: tropLev            ! tropopause level index
        1380             :     real(kind_phys), optional, intent(in)      :: tropP              ! tropopause pressure (Pa)
        1381             :     real(kind_phys)                            :: tropopause_interpolateT
        1382             :
        1383             :     ! Local Variables
        1384             :     real(kind_phys)   :: tropT              ! tropopause temperature (K)
        1385             :     real(kind_phys)   :: dTdlogP            ! dT/dlog(P)
        1386             :
        1387             :     ! Interpolate the temperature linearly against log(P)
        1388             :
        1389             :     ! Is the tropopause at the midpoint?
        1390    46599902 :     if (tropP == pmid(icol, tropLev)) then
        1391           0 :       tropT = t(icol, tropLev)
        1392             :
        1393    46599902 :     else if (tropP < pmid(icol, tropLev)) then
        1394             :
        1395             :       ! It is above the midpoint? Make sure we aren't at the top.
        1396    22334491 :       if (tropLev > 1) then
        1397    44668982 :         dTdlogP = (t(icol, tropLev) - t(icol, tropLev - 1)) / &
        1398    44668982 :           (log(pmid(icol, tropLev)) - log(pmid(icol, tropLev - 1)))
        1399    22334491 :         tropT = t(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dTdlogP
        1400             :       end if
        1401             :     else
        1402             :
        1403             :       ! It is below the midpoint. Make sure we aren't at the bottom.
        1404    24265411 :       if (tropLev < pver) then
        1405    24265411 :         dTdlogP = (t(icol, tropLev + 1) - t(icol, tropLev)) / &
        1406    48530822 :           (log(pmid(icol, tropLev + 1)) - log(pmid(icol, tropLev)))
        1407    24265411 :         tropT = t(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dTdlogP
        1408             :       end if
        1409             :     end if
        1410             :
        1411    46599902 :     tropopause_interpolateT = tropT
        1412    46599902 :   end function tropopause_interpolateT
        1413             :
        1414             :
        1415             :   ! This routine interpolates the geopotential height in the physics state to
        1416             :   ! find the geopotential height at the specified tropopause pressure.
        1417    46599902 :   function tropopause_interpolateZ(pint, pmid, zi, zm, phis, icol, tropLev, tropP)
        1418             :
        1419             :     real(kind_phys), intent(in)         :: pint(:,:)     ! Interface pressures (Pa)
        1420             :     real(kind_phys), intent(in)         :: pmid(:,:)     ! Midpoint pressures (Pa)
        1421             :     real(kind_phys), intent(in)         :: zi(:,:)       ! Geopotential height above surface at interfaces (m)
        1422             :     real(kind_phys), intent(in)         :: zm(:,:)       ! Geopotential height above surface at midpoints (m)
        1423             :     real(kind_phys), intent(in)         :: phis(:)       ! Surface geopotential (m2 s-2)
        1424             :     integer, intent(in)                 :: icol               ! column being processed
        1425             :     integer, intent(in)                 :: tropLev            ! tropopause level index
        1426             :     real(kind_phys), optional, intent(in)      :: tropP              ! tropopause pressure (Pa)
        1427             :     real(kind_phys)                            :: tropopause_interpolateZ
        1428             :
        1429             :     ! Local Variables
        1430             :     real(kind_phys)   :: tropZ              ! tropopause geopotential height (m)
        1431             :     real(kind_phys)   :: dZdlogP            ! dZ/dlog(P)
        1432             :
        1433             :     ! Interpolate the geopotential height linearly against log(P)
        1434             :
        1435             :     ! Is the tropopause at the midpoint?
        1436    46599902 :     if (tropP == pmid(icol, tropLev)) then
        1437           0 :       tropZ = zm(icol, tropLev)
        1438             :
        1439    46599902 :     else if (tropP < pmid(icol, tropLev)) then
        1440             :
        1441             :       ! It is above the midpoint? Make sure we aren't at the top.
        1442    22334491 :       dZdlogP = (zm(icol, tropLev) - zi(icol, tropLev)) / &
        1443    22334491 :         (log(pmid(icol, tropLev)) - log(pint(icol, tropLev)))
        1444    22334491 :       tropZ = zm(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dZdlogP
        1445             :     else
        1446             :
        1447             :       ! It is below the midpoint. Make sure we aren't at the bottom.
        1448    48530822 :       dZdlogP = (zm(icol, tropLev) - zi(icol, tropLev+1)) / &
        1449    48530822 :         (log(pmid(icol, tropLev)) - log(pint(icol, tropLev+1)))
        1450    24265411 :       tropZ = zm(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dZdlogP
        1451             :     end if
        1452             :
        1453    46599902 :     tropopause_interpolateZ = tropZ + phis(icol)*cnst_rga
        1454    46599902 :   end function tropopause_interpolateZ
        1455             : end module tropopause_find

Generated by: LCOV version 1.14