LCOV - code coverage report
Current view: | top level - physics/cam - tropopause.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 165 | 195 | 84.6 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 6 | 6 | 100.0 %

              Line data    Source code
           1             : ! This is the CAM interface to the CCPP-ized tropopause_find scheme.
           2             : ! Full compatibility, bit-for-bit, to old CAM approach is achieved through
           3             : ! this module, however this module will not be necessary in CAM-SIMA.
           4             : !
           5             : ! For science description of the underlying algorithms, refer to
           6             : ! atmospheric_physics/tropopause_find/tropopause_find.F90.
           7             : ! (hplin, 8/20/24)
           8             :
           9             : module tropopause
          10             :   !---------------------------------------------------------------
          11             :   ! ... variables for the tropopause module
          12             :   !---------------------------------------------------------------
          13             :
          14             :   use shr_kind_mod,         only : r8 => shr_kind_r8
          15             :   use shr_const_mod,        only : pi => shr_const_pi
          16             :   use ppgrid,               only : pcols, pver, pverp, begchunk, endchunk
          17             :   use cam_abortutils,       only : endrun
          18             :   use cam_logfile,          only : iulog
          19             :   use cam_history_support,  only : fillvalue
          20             :   use physics_types,        only : physics_state
          21             :   use spmd_utils,           only : masterproc
          22             :
          23             :   implicit none
          24             :
          25             :   private
          26             :
          27             :   public  :: tropopause_readnl, tropopause_init, tropopause_find_cam, tropopause_output
          28             :   public  :: tropopause_findChemTrop
          29             :   public  :: TROP_ALG_NONE, TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
          30             :   public  :: TROP_ALG_STOBIE, TROP_ALG_HYBSTOB, TROP_ALG_TWMO, TROP_ALG_WMO
          31             :   public  :: TROP_ALG_CPP
          32             :   public  :: NOTFOUND
          33             :
          34             :   save
          35             :
          36             :   ! These parameters define and enumeration to be used to define the primary
          37             :   ! and backup algorithms to be used with the tropopause_find() method. The
          38             :   ! backup algorithm is meant to provide a solution when the primary algorithm
          39             :   ! fail. The algorithms that can't fail are: TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
          40             :   ! and TROP_ALG_STOBIE.
          41             :   integer, parameter    :: TROP_ALG_NONE      = 1    ! Don't evaluate
          42             :   integer, parameter    :: TROP_ALG_ANALYTIC  = 2    ! Analytic Expression
          43             :   integer, parameter    :: TROP_ALG_CLIMATE   = 3    ! Climatology
          44             :   integer, parameter    :: TROP_ALG_STOBIE    = 4    ! Stobie Algorithm
          45             :   integer, parameter    :: TROP_ALG_TWMO      = 5    ! WMO Definition, Reichler et al. [2003]
          46             :   integer, parameter    :: TROP_ALG_WMO       = 6    ! WMO Definition
          47             :   integer, parameter    :: TROP_ALG_HYBSTOB   = 7    ! Hybrid Stobie Algorithm
          48             :   integer, parameter    :: TROP_ALG_CPP       = 8    ! Cold Point Parabolic
          49             :   integer, parameter    :: TROP_ALG_CHEMTROP  = 9    ! Chemical tropopause
          50             :
          51             :   ! Note: exclude CHEMTROP here as it is a new flag added in CCPP-ized routines to unify the chemTrop routine. (hplin, 8/20/24)
          52             :   integer, parameter    :: TROP_NALG          = 8    ! Number of Algorithms
          53             :   character,parameter   :: TROP_LETTER(TROP_NALG) = (/ ' ', 'A', 'C', 'S', 'T', 'W', 'H', 'F' /)
          54             :                                                      ! unique identifier for output, don't use P
          55             :
          56             :   ! These variables should probably be controlled by namelist entries.
          57             :   logical ,parameter    :: output_all         = .False.              ! output tropopause info from all algorithms
          58             :   integer ,parameter    :: default_primary    = TROP_ALG_TWMO        ! default primary algorithm
          59             :   integer ,parameter    :: default_backup     = TROP_ALG_CLIMATE     ! default backup algorithm
          60             :
          61             :   ! Namelist variables
          62             :   character(len=256)    :: tropopause_climo_file = 'trop_climo'      ! absolute filepath of climatology file
          63             :
          64             :   ! These variables are used to store the climatology data.
          65             :   real(r8)              :: days(12)                                  ! days in the climatology
          66             :   real(r8), pointer     :: tropp_p_loc(:,:,:)                        ! climatological tropopause pressures
          67             :
          68             :   integer, parameter :: NOTFOUND = -1
          69             :
          70             :   ! physical constants
          71             :   ! These constants are set in module variables rather than as parameters
          72             :   ! to support the aquaplanet mode in which the constants have values determined
          73             :   ! by the experiment protocol
          74             :   real(r8) :: cnst_kap     ! = cappa
          75             :   real(r8) :: cnst_faktor  ! = -gravit/rair
          76             :   real(r8) :: cnst_ka1     ! = cnst_kap - 1._r8
          77             :
          78             : !================================================================================================
          79             : contains
          80             : !================================================================================================
          81             :
          82             :    ! Read namelist variables.
          83     1490712 :    subroutine tropopause_readnl(nlfile)
          84             :
          85             :       use namelist_utils,  only: find_group_name
          86             :       use units,           only: getunit, freeunit
          87             :       use mpishorthand
          88             :
          89             :       character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          90             :
          91             :       ! Local variables
          92             :       integer :: unitn, ierr
          93             :       character(len=*), parameter :: subname = 'tropopause_readnl'
          94             :
          95             :       namelist /tropopause_nl/ tropopause_climo_file
          96             :       !-----------------------------------------------------------------------------
          97             :
          98        1536 :       if (masterproc) then
          99           2 :          unitn = getunit()
         100           2 :          open( unitn, file=trim(nlfile), status='old' )
         101           2 :          call find_group_name(unitn, 'tropopause_nl', status=ierr)
         102           2 :          if (ierr == 0) then
         103           2 :             read(unitn, tropopause_nl, iostat=ierr)
         104           2 :             if (ierr /= 0) then
         105           0 :                call endrun(subname // ':: ERROR reading namelist')
         106             :             end if
         107             :          end if
         108           2 :          close(unitn)
         109           2 :          call freeunit(unitn)
         110             :       end if
         111             :
         112             : #ifdef SPMD
         113             :       ! Broadcast namelist variables
         114        1536 :       call mpibcast(tropopause_climo_file, len(tropopause_climo_file), mpichar, 0, mpicom)
         115             : #endif
         116             :
         117        1536 :    end subroutine tropopause_readnl
         118             :
         119             :
         120             :   ! This routine is called during intialization and must be called before the
         121             :   ! other methods in this module can be used. Its main tasks are to read in the
         122             :   ! climatology from a file and to define the output fields. Much of this code
         123             :   ! is taken from mo_tropopause.
         124        1536 :   subroutine tropopause_init()
         125             :
         126             :     use cam_history,     only: addfld, horiz_only
         127             :     use tropopause_find, only: tropopause_find_init
         128             :     use physconst,       only: cappa, rair, gravit, pi
         129             :
         130             :     character(len=512) :: errmsg
         131             :     integer            :: errflg
         132             :
         133             :     ! Call underlying CCPP-initialization routine.
         134        1536 :     call tropopause_find_init(cappa, rair, gravit, pi, errmsg, errflg)
         135             :
         136             :     ! Define the output fields.
         137        1536 :     call addfld('TROP_P',          horiz_only,  'A', 'Pa',          'Tropopause Pressure',              flag_xyfill=.True.)
         138        1536 :     call addfld('TROP_T',          horiz_only,  'A', 'K',           'Tropopause Temperature',           flag_xyfill=.True.)
         139        1536 :     call addfld('TROP_Z',          horiz_only,  'A', 'm',           'Tropopause Height',                flag_xyfill=.True.)
         140        3072 :     call addfld('TROP_DZ',         (/ 'lev' /), 'A', 'm',           'Relative Tropopause Height')
         141        3072 :     call addfld('TROP_PD',         (/ 'lev' /), 'A', 'probability', 'Tropopause Probabilty')
         142        1536 :     call addfld('TROP_FD',         horiz_only,  'A', 'probability', 'Tropopause Found')
         143             :
         144        1536 :     call addfld('TROPP_P',         horiz_only,  'A', 'Pa',          'Tropopause Pressure (primary)',    flag_xyfill=.True.)
         145        1536 :     call addfld('TROPP_T',         horiz_only,  'A', 'K',           'Tropopause Temperature (primary)', flag_xyfill=.True.)
         146        1536 :     call addfld('TROPP_Z',         horiz_only,  'A', 'm',           'Tropopause Height (primary)',      flag_xyfill=.True.)
         147        3072 :     call addfld('TROPP_DZ',        (/ 'lev' /), 'A', 'm',           'Relative Tropopause Height (primary)')
         148        3072 :     call addfld('TROPP_PD',        (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (primary)')
         149        1536 :     call addfld('TROPP_FD',        horiz_only,  'A', 'probability', 'Tropopause Found (primary)')
         150             :
         151        1536 :     call addfld('TROPF_P',         horiz_only,  'A',  'Pa',         'Tropopause Pressure (cold point)',    flag_xyfill=.True.)
         152        1536 :     call addfld('TROPF_T',         horiz_only,  'A',  'K',          'Tropopause Temperature (cold point)', flag_xyfill=.True.)
         153        1536 :     call addfld('TROPF_Z',         horiz_only,  'A',  'm',          'Tropopause Height (cold point)',      flag_xyfill=.True.)
         154        3072 :     call addfld('TROPF_DZ',        (/ 'lev' /),  'A', 'm',          'Relative Tropopause Height (cold point)', flag_xyfill=.True.)
         155        3072 :     call addfld('TROPF_PD',        (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (cold point)')
         156        1536 :     call addfld('TROPF_FD',        horiz_only,  'A', 'probability', 'Tropopause Found (cold point)')
         157             :
         158        3072 :     call addfld( 'hstobie_trop',   (/ 'lev' /), 'I',  'fraction of model time', 'Lowest level with stratospheric chemsitry')
         159        3072 :     call addfld( 'hstobie_linoz',  (/ 'lev' /), 'I',  'fraction of model time', 'Lowest possible Linoz level')
         160             :     call addfld( 'hstobie_tropop', (/ 'lev' /), 'I', 'fraction of model time', &
         161        3072 :          'Troposphere boundary calculated in chemistry' )
         162             :
         163             :     ! If requested, be prepared to output results from all of the methods.
         164             :     if (output_all) then
         165             :       call addfld('TROPA_P',  horiz_only,  'A',  'Pa',          'Tropopause Pressure (analytic)',        flag_xyfill=.True.)
         166             :       call addfld('TROPA_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (analytic)',      flag_xyfill=.True.)
         167             :       call addfld('TROPA_Z',  horiz_only,  'A',  'm',          'Tropopause Height (analytic)',           flag_xyfill=.True.)
         168             :       call addfld('TROPA_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (analytic)')
         169             :       call addfld('TROPA_FD', horiz_only,  'A', 'probability', 'Tropopause Found (analytic)')
         170             :
         171             :       call addfld('TROPC_P',  horiz_only,  'A',  'Pa',         'Tropopause Pressure (climatology)',      flag_xyfill=.True.)
         172             :       call addfld('TROPC_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (climatology)',   flag_xyfill=.True.)
         173             :       call addfld('TROPC_Z',  horiz_only,  'A',  'm',          'Tropopause Height (climatology)',        flag_xyfill=.True.)
         174             :       call addfld('TROPC_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (climatology)')
         175             :       call addfld('TROPC_FD', horiz_only,  'A', 'probability', 'Tropopause Found (climatology)')
         176             :
         177             :       call addfld('TROPS_P',  horiz_only,  'A',  'Pa',         'Tropopause Pressure (stobie)',           flag_xyfill=.True.)
         178             :       call addfld('TROPS_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (stobie)',        flag_xyfill=.True.)
         179             :       call addfld('TROPS_Z',  horiz_only,  'A',  'm',          'Tropopause Height (stobie)',             flag_xyfill=.True.)
         180             :       call addfld('TROPS_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (stobie)')
         181             :       call addfld('TROPS_FD', horiz_only,  'A', 'probability', 'Tropopause Found (stobie)')
         182             :
         183             :       call addfld('TROPT_P',  horiz_only,  'A',  'Pa',         'Tropopause Pressure (twmo)',             flag_xyfill=.True.)
         184             :       call addfld('TROPT_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (twmo)',          flag_xyfill=.True.)
         185             :       call addfld('TROPT_Z',  horiz_only,  'A',  'm',          'Tropopause Height (twmo)',               flag_xyfill=.True.)
         186             :       call addfld('TROPT_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (twmo)')
         187             :       call addfld('TROPT_FD', horiz_only,  'A', 'probability', 'Tropopause Found (twmo)')
         188             :
         189             :       call addfld('TROPW_P',  horiz_only,  'A',  'Pa',         'Tropopause Pressure (WMO)',              flag_xyfill=.True.)
         190             :       call addfld('TROPW_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (WMO)',           flag_xyfill=.True.)
         191             :       call addfld('TROPW_Z',  horiz_only,  'A',  'm',          'Tropopause Height (WMO)',                flag_xyfill=.True.)
         192             :       call addfld('TROPW_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (WMO)')
         193             :       call addfld('TROPW_FD', horiz_only,  'A', 'probability', 'Tropopause Found (WMO)')
         194             :
         195             :       call addfld('TROPH_P',  horiz_only,  'A',  'Pa',         'Tropopause Pressure (Hybrid Stobie)',    flag_xyfill=.True.)
         196             :       call addfld('TROPH_T',  horiz_only,  'A',  'K',          'Tropopause Temperature (Hybrid Stobie)', flag_xyfill=.True.)
         197             :       call addfld('TROPH_Z',  horiz_only,  'A',  'm',          'Tropopause Height (Hybrid Stobie)',      flag_xyfill=.True.)
         198             :       call addfld('TROPH_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (Hybrid Stobie)')
         199             :       call addfld('TROPH_FD', horiz_only,  'A', 'probability', 'Tropopause Found (Hybrid Stobie)')
         200             :     end if
         201             :
         202             :
         203        1536 :     call tropopause_read_file()
         204             :
         205             :
         206        1536 :   end subroutine tropopause_init
         207             :
         208             :
         209        1536 :   subroutine tropopause_read_file
         210             :     !------------------------------------------------------------------
         211             :     ! ... initialize upper boundary values
         212             :     !------------------------------------------------------------------
         213        1536 :     use interpolate_data,  only : lininterp_init, lininterp, interp_type, lininterp_finish
         214             :     use dyn_grid,     only : get_dyn_grid_parm
         215             :     use phys_grid,    only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
         216             :     use ioFileMod,    only : getfil
         217             :     use time_manager, only : get_calday
         218             :     use physconst,    only : pi
         219             :     use cam_pio_utils, only: cam_pio_openfile
         220             :     use pio,          only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, &
         221             :          pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite
         222             :
         223             :     !------------------------------------------------------------------
         224             :     ! ... local variables
         225             :     !------------------------------------------------------------------
         226             :     integer :: i, j, n
         227             :     integer :: ierr
         228             :     type(file_desc_t) :: pio_id
         229             :     integer :: dimid
         230             :     type(var_desc_t) :: vid
         231             :     integer :: nlon, nlat, ntimes
         232             :     integer :: start(3)
         233             :     integer :: count(3)
         234             :     integer, parameter :: dates(12) = (/ 116, 214, 316, 415,  516,  615, &
         235             :          716, 816, 915, 1016, 1115, 1216 /)
         236             :     integer :: plon, plat
         237             :     type(interp_type) :: lon_wgts, lat_wgts
         238        1536 :     real(r8), allocatable :: tropp_p_in(:,:,:)
         239        1536 :     real(r8), allocatable :: lat(:)
         240        1536 :     real(r8), allocatable :: lon(:)
         241             :     real(r8) :: to_lats(pcols), to_lons(pcols)
         242             :     real(r8), parameter :: d2r=pi/180._r8, zero=0._r8, twopi=pi*2._r8
         243             :     character(len=256) :: locfn
         244             :     integer  :: c, ncols
         245             :
         246             :
         247             :     plon = get_dyn_grid_parm('plon')
         248             :     plat = get_dyn_grid_parm('plat')
         249             :
         250             :
         251             :     !-----------------------------------------------------------------------
         252             :     !       ... open netcdf file
         253             :     !-----------------------------------------------------------------------
         254        1536 :     call getfil (tropopause_climo_file, locfn, 0)
         255        1536 :     call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE)
         256             :
         257             :     !-----------------------------------------------------------------------
         258             :     !       ... get time dimension
         259             :     !-----------------------------------------------------------------------
         260        1536 :     ierr = pio_inq_dimid( pio_id, 'time', dimid )
         261        1536 :     ierr = pio_inq_dimlen( pio_id, dimid, ntimes )
         262        1536 :     if( ntimes /= 12 )then
         263           0 :        write(iulog,*) 'tropopause_init: number of months = ',ntimes,'; expecting 12'
         264           0 :        call endrun
         265             :     end if
         266             :     !-----------------------------------------------------------------------
         267             :     !       ... get latitudes
         268             :     !-----------------------------------------------------------------------
         269        1536 :     ierr = pio_inq_dimid( pio_id, 'lat', dimid )
         270        1536 :     ierr = pio_inq_dimlen( pio_id, dimid, nlat )
         271        4608 :     allocate( lat(nlat), stat=ierr )
         272        1536 :     if( ierr /= 0 ) then
         273           0 :        write(iulog,*) 'tropopause_init: lat allocation error = ',ierr
         274           0 :        call endrun
         275             :     end if
         276        1536 :     ierr = pio_inq_varid( pio_id, 'lat', vid )
         277        1536 :     ierr = pio_get_var( pio_id, vid, lat )
         278      115200 :     lat(:nlat) = lat(:nlat) * d2r
         279             :     !-----------------------------------------------------------------------
         280             :     !       ... get longitudes
         281             :     !-----------------------------------------------------------------------
         282        1536 :     ierr = pio_inq_dimid( pio_id, 'lon', dimid )
         283        1536 :     ierr = pio_inq_dimlen( pio_id, dimid, nlon )
         284        4608 :     allocate( lon(nlon), stat=ierr )
         285        1536 :     if( ierr /= 0 ) then
         286           0 :        write(iulog,*) 'tropopause_init: lon allocation error = ',ierr
         287           0 :        call endrun
         288             :     end if
         289        1536 :     ierr = pio_inq_varid( pio_id, 'lon', vid )
         290        1536 :     ierr = pio_get_var( pio_id, vid, lon )
         291      224256 :     lon(:nlon) = lon(:nlon) * d2r
         292             :
         293             :     !------------------------------------------------------------------
         294             :     !  ... allocate arrays
         295             :     !------------------------------------------------------------------
         296        7680 :     allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr )
         297        1536 :     if( ierr /= 0 ) then
         298           0 :        write(iulog,*) 'tropopause_init: tropp_p_in allocation error = ',ierr
         299           0 :        call endrun
         300             :     end if
         301             :     !------------------------------------------------------------------
         302             :     !  ... read in the tropopause pressure
         303             :     !------------------------------------------------------------------
         304        1536 :     ierr = pio_inq_varid( pio_id, 'trop_p', vid )
         305        1536 :     start = (/ 1, 1, 1 /)
         306        6144 :     count = (/ nlon, nlat, ntimes /)
         307        1536 :     ierr = pio_get_var( pio_id, vid, start, count, tropp_p_in )
         308             :
         309             :     !------------------------------------------------------------------
         310             :     !  ... close the netcdf file
         311             :     !------------------------------------------------------------------
         312        1536 :     call pio_closefile( pio_id )
         313             :
         314             :     !--------------------------------------------------------------------
         315             :     !  ... regrid
         316             :     !--------------------------------------------------------------------
         317             :
         318        6144 :     allocate( tropp_p_loc(pcols,begchunk:endchunk,ntimes), stat=ierr )
         319             :
         320        1536 :     if( ierr /= 0 ) then
         321           0 :       write(iulog,*) 'tropopause_init: tropp_p_loc allocation error = ',ierr
         322           0 :       call endrun
         323             :     end if
         324             :
         325        7728 :     do c=begchunk,endchunk
         326        6192 :        ncols = get_ncols_p(c)
         327        6192 :        call get_rlat_all_p(c, pcols, to_lats)
         328        6192 :        call get_rlon_all_p(c, pcols, to_lons)
         329        6192 :        call lininterp_init(lon, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
         330        6192 :        call lininterp_init(lat, nlat, to_lats, ncols, 1, lat_wgts)
         331       80496 :        do n=1,ntimes
         332       80496 :           call lininterp(tropp_p_in(:,:,n), nlon, nlat, tropp_p_loc(1:ncols,c,n), ncols, lon_wgts, lat_wgts)
         333             :        end do
         334        6192 :        call lininterp_finish(lon_wgts)
         335        7728 :        call lininterp_finish(lat_wgts)
         336             :     end do
         337        1536 :     deallocate(lon)
         338        1536 :     deallocate(lat)
         339        1536 :     deallocate(tropp_p_in)
         340             :
         341             :     !--------------------------------------------------------
         342             :     ! ... initialize the monthly day of year times
         343             :     !--------------------------------------------------------
         344             :
         345       19968 :     do n = 1,12
         346       19968 :        days(n) = get_calday( dates(n), 0 )
         347             :     end do
         348        1536 :     if (masterproc) then
         349           2 :        write(iulog,*) 'tropopause_init : days'
         350           2 :        write(iulog,'(1p,5g15.8)') days(:)
         351             :     endif
         352             :
         353        3072 :   end subroutine tropopause_read_file
         354             :
         355             :   ! Searches all the columns in the chunk and attempts to identify the tropopause.
         356             :   ! Two routines can be specifed, a primary routine which is tried first and a
         357             :   ! backup routine which will be tried only if the first routine fails. If the
         358             :   ! tropopause can not be identified by either routine, then a NOTFOUND is returned
         359             :   ! for the tropopause level, temperature and pressure.
         360     6702840 :   subroutine tropopause_find_cam(pstate, tropLev, tropP, tropT, tropZ, primary, backup)
         361             :
         362        1536 :     use tropopause_find, only: tropopause_findWithBackup
         363             :
         364             :     use cam_history,     only: outfld
         365             :     use time_manager,    only: get_curr_calday
         366             :
         367             :     implicit none
         368             :
         369             :     type(physics_state), intent(in)     :: pstate
         370             :     integer, optional, intent(in)       :: primary                   ! primary detection algorithm
         371             :     integer, optional, intent(in)       :: backup                    ! backup detection algorithm
         372             :     integer,            intent(out)     :: tropLev(:)                ! tropopause level index
         373             :     real(r8), optional, intent(out)     :: tropP(:)                  ! tropopause pressure (Pa)
         374             :     real(r8), optional, intent(out)     :: tropT(:)                  ! tropopause temperature (K)
         375             :     real(r8), optional, intent(out)     :: tropZ(:)                  ! tropopause height (m)
         376             :
         377             :     ! Local Variable
         378             :     integer       :: primAlg            ! Primary algorithm
         379             :     integer       :: backAlg            ! Backup algorithm
         380             :
         381             :     real(r8)      :: calday
         382             :     integer       :: ncol
         383             :
         384             :     real(r8)      :: hstobie_trop  (pcols, pver)
         385             :     real(r8)      :: hstobie_linoz (pcols, pver)
         386             :     real(r8)      :: hstobie_tropop(pcols, pver)
         387             :
         388             :     character(len=512) :: errmsg
         389             :     integer            :: errflg
         390             :
         391             :     ! Get compatibility variables for CCPP-ized routine
         392     6702840 :     ncol   = pstate%ncol
         393     6702840 :     calday = get_curr_calday()
         394             :
         395             :     ! Initialize the results to a missing value, so that the algorithms will
         396             :     ! attempt to find the tropopause for all of them. Only do this for the active columns.
         397   118624680 :     tropLev(:ncol) = NOTFOUND
         398    81300168 :     if (present(tropP)) tropP(:ncol) = fillvalue
         399    81300168 :     if (present(tropT)) tropT(:ncol) = fillvalue
         400    81300168 :     if (present(tropZ)) tropZ(:ncol) = fillvalue
         401             :
         402             :     ! Set the algorithms to be used, either the ones provided or the defaults.
         403     6702840 :     if (present(primary)) then
         404     1489176 :       primAlg = primary
         405             :     else
         406     5213664 :       primAlg = default_primary
         407             :     end if
         408             :
         409     6702840 :     if (present(backup)) then
         410     2978352 :       backAlg = backup
         411             :     else
         412     3724488 :       backAlg = default_backup
         413             :     end if
         414             :
         415             :     ! This does not call the tropopause_find_run routine directly, because it
         416             :     ! computes multiple needed tropopauses simultaneously. Instead, here we
         417             :     ! specify the algorithm needed directly to the algorithm driver routine.
         418             :     call tropopause_findWithBackup( &
         419             :          ncol           = ncol, &
         420             :          pver           = pver, &
         421             :          fillvalue      = fillvalue, &
         422           0 :          lat            = pstate%lat(:ncol), &
         423           0 :          pint           = pstate%pint(:ncol, :pverp), &
         424           0 :          pmid           = pstate%pmid(:ncol, :pver), &
         425           0 :          t              = pstate%t(:ncol, :pver), &
         426           0 :          zi             = pstate%zi(:ncol, :pverp), &
         427           0 :          zm             = pstate%zm(:ncol, :pver), &
         428           0 :          phis           = pstate%phis(:ncol), &
         429             :          calday         = calday, &
         430           0 :          tropp_p_loc    = tropp_p_loc(:ncol,pstate%lchnk,:), &  ! Subset into chunk as the underlying routines are no longer chunkized.
         431             :          tropp_days     = days, &
         432     6702840 :          tropLev        = tropLev(:ncol), &
         433             :          tropP          = tropP, &
         434             :          tropT          = tropT, &
         435             :          tropZ          = tropZ, &
         436             :          primary        = primAlg, &
         437             :          backup         = backAlg, &
         438     6702840 :          hstobie_trop   = hstobie_trop(:ncol, :pver), &    ! Only used if TROP_ALG_HYBSTOB
         439             :          hstobie_linoz  = hstobie_linoz(:ncol, :pver), &   ! Only used if TROP_ALG_HYBSTOB
         440             :          hstobie_tropop = hstobie_tropop(:ncol, :pver), &  ! Only used if TROP_ALG_HYBSTOB
         441             :          errmsg         = errmsg, &
         442             :          errflg         = errflg &
         443    24579144 :     )
         444             :
         445             :     ! Output hybridstobie specific fields
         446     6702840 :     if(primAlg == TROP_ALG_HYBSTOB) then
         447           0 :        call outfld('hstobie_trop',   hstobie_trop(:ncol,:),    ncol, pstate%lchnk )
         448           0 :        call outfld('hstobie_linoz',  hstobie_linoz(:ncol,:),   ncol, pstate%lchnk )
         449           0 :        call outfld('hstobie_tropop', hstobie_tropop(:ncol,:),  ncol, pstate%lchnk )
         450             :     endif
         451     6702840 :   end subroutine tropopause_find_cam
         452             :
         453             :   ! Searches all the columns in the chunk and attempts to identify the "chemical"
         454             :   ! tropopause. This is the lapse rate tropopause, backed up by the climatology
         455             :   ! if the lapse rate fails to find the tropopause at pressures higher than a certain
         456             :   ! threshold. This pressure threshold depends on latitude. Between 50S and 50N,
         457             :   ! the climatology is used if the lapse rate tropopause is not found at P > 75 hPa.
         458             :   ! At high latitude (poleward of 50), the threshold is increased to 125 hPa to
         459             :   ! eliminate false events that are sometimes detected in the cold polar stratosphere.
         460             :   !
         461             :   ! NOTE: This routine was adapted from code in chemistry.F90 and mo_gasphase_chemdr.F90.
         462    12665736 :   subroutine tropopause_findChemTrop(pstate, tropLev)
         463             :
         464     6702840 :     use tropopause_find, only: tropopause_findWithBackup
         465             :
         466             :     use time_manager,    only: get_curr_calday
         467             :
         468             :     implicit none
         469             :
         470             :     type(physics_state), intent(in)     :: pstate
         471             :     integer,             intent(out)    :: tropLev(:)            ! tropopause level index
         472             :
         473             :     ! Local Variable
         474             :     real(r8)            :: calday
         475             :     integer             :: i
         476             :     integer             :: ncol
         477             :
         478             :     character(len=512) :: errmsg
         479             :     integer            :: errflg
         480             :
         481             :     ! Get compatibility variables for CCPP-ized routine
         482    12665736 :     ncol   = pstate%ncol
         483    12665736 :     calday = get_curr_calday()
         484             :
         485             :     ! Now call the unified routine with the CHEMTROP option, which has automatic
         486             :     ! backup fall to climatology.
         487             :     call tropopause_findWithBackup( &
         488             :          ncol           = ncol, &
         489             :          pver           = pver, &
         490             :          fillvalue      = fillvalue, &
         491           0 :          lat            = pstate%lat(:ncol), &
         492           0 :          pint           = pstate%pint(:ncol, :pverp), &
         493           0 :          pmid           = pstate%pmid(:ncol, :pver), &
         494           0 :          t              = pstate%t(:ncol, :pver), &
         495           0 :          zi             = pstate%zi(:ncol, :pverp), &
         496           0 :          zm             = pstate%zm(:ncol, :pver), &
         497           0 :          phis           = pstate%phis(:ncol), &
         498             :          calday         = calday, &
         499           0 :          tropp_p_loc    = tropp_p_loc(:ncol,pstate%lchnk,:), &  ! Subset into chunk as the underlying routines are no longer chunkized.
         500             :          tropp_days     = days, &
         501    12665736 :          tropLev        = tropLev(1:ncol), &
         502             :          primary        = TROP_ALG_CHEMTROP, &
         503             :          backup         = TROP_ALG_CLIMATE, &
         504             :          errmsg         = errmsg, &
         505             :          errflg         = errflg &
         506    25331472 :     )
         507    12665736 :   end subroutine tropopause_findChemTrop
         508             :
         509             :   ! Output the tropopause pressure and temperature to the history files. Two sets
         510             :   ! of output will be generated, one for the default algorithm and another one
         511             :   ! using the default routine, but backed by a climatology when the default
         512             :   ! algorithm fails.
         513     1489176 :   subroutine tropopause_output(pstate)
         514    12665736 :     use cam_history,  only : outfld
         515             :
         516             :     implicit none
         517             :
         518             :     type(physics_state), intent(in)     :: pstate
         519             :
         520             :     ! Local Variables
         521             :     integer       :: i
         522             :     integer       :: alg
         523             :     integer       :: ncol                     ! number of cloumns in the chunk
         524             :     integer       :: lchnk                    ! chunk identifier
         525             :     integer       :: tropLev(pcols)           ! tropopause level index
         526             :     real(r8)      :: tropP(pcols)             ! tropopause pressure (Pa)
         527             :     real(r8)      :: tropT(pcols)             ! tropopause temperature (K)
         528             :     real(r8)      :: tropZ(pcols)             ! tropopause height (m)
         529             :     real(r8)      :: tropFound(pcols)         ! tropopause found
         530             :     real(r8)      :: tropDZ(pcols, pver)      ! relative tropopause height (m)
         531             :     real(r8)      :: tropPdf(pcols, pver)     ! tropopause probability distribution
         532             :
         533             :     ! Information about the chunk.
         534     1489176 :     lchnk = pstate%lchnk
         535     1489176 :     ncol  = pstate%ncol
         536             :
         537             :     ! Find the tropopause using the default algorithm backed by the climatology.
         538     1489176 :     call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ)
         539             :
         540     1489176 :     tropPdf(:,:) = 0._r8
         541     1489176 :     tropFound(:) = 0._r8
         542  2355876432 :     tropDZ(:,:) = fillvalue
         543    24865776 :     do i = 1, ncol
         544    24865776 :       if (tropLev(i) /= NOTFOUND) then
         545    23376600 :         tropPdf(i, tropLev(i)) = 1._r8
         546    23376600 :         tropFound(i) = 1._r8
         547  2197400400 :         tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
         548             :       end if
         549             :     end do
         550             :
         551     1489176 :     call outfld('TROP_P',   tropP(:ncol),      ncol, lchnk)
         552     1489176 :     call outfld('TROP_T',   tropT(:ncol),      ncol, lchnk)
         553     1489176 :     call outfld('TROP_Z',   tropZ(:ncol),      ncol, lchnk)
         554   543653136 :     call outfld('TROP_DZ',  tropDZ(:ncol, :), ncol, lchnk)
         555   543653136 :     call outfld('TROP_PD',  tropPdf(:ncol, :), ncol, lchnk)
         556     1489176 :     call outfld('TROP_FD',  tropFound(:ncol),  ncol, lchnk)
         557             :
         558             :
         559             :     ! Find the tropopause using just the primary algorithm.
         560     1489176 :     call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, backup=TROP_ALG_NONE)
         561             :
         562     1489176 :     tropPdf(:,:) = 0._r8
         563     1489176 :     tropFound(:) = 0._r8
         564  2355876432 :     tropDZ(:,:) = fillvalue
         565             :
         566    24865776 :     do i = 1, ncol
         567    24865776 :       if (tropLev(i) /= NOTFOUND) then
         568    23205946 :         tropPdf(i, tropLev(i)) = 1._r8
         569    23205946 :         tropFound(i) = 1._r8
         570  2181358924 :         tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
         571             :       end if
         572             :     end do
         573             :
         574     1489176 :     call outfld('TROPP_P',   tropP(:ncol),      ncol, lchnk)
         575     1489176 :     call outfld('TROPP_T',   tropT(:ncol),      ncol, lchnk)
         576     1489176 :     call outfld('TROPP_Z',   tropZ(:ncol),      ncol, lchnk)
         577   543653136 :     call outfld('TROPP_DZ',  tropDZ(:ncol, :), ncol, lchnk)
         578   543653136 :     call outfld('TROPP_PD',  tropPdf(:ncol, :), ncol, lchnk)
         579     1489176 :     call outfld('TROPP_FD',  tropFound(:ncol),  ncol, lchnk)
         580             :
         581             :
         582             :     ! Find the tropopause using just the cold point algorithm.
         583     1489176 :     call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=TROP_ALG_CPP, backup=TROP_ALG_NONE)
         584             :
         585     1489176 :     tropPdf(:,:) = 0._r8
         586     1489176 :     tropFound(:) = 0._r8
         587  2355876432 :     tropDZ(:,:) = fillvalue
         588             :
         589    24865776 :     do i = 1, ncol
         590    24865776 :       if (tropLev(i) /= NOTFOUND) then
         591    23101516 :         tropPdf(i, tropLev(i)) = 1._r8
         592    23101516 :         tropFound(i) = 1._r8
         593  2171542504 :         tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
         594             :       end if
         595             :     end do
         596             :
         597     1489176 :     call outfld('TROPF_P',   tropP(:ncol),      ncol, lchnk)
         598     1489176 :     call outfld('TROPF_T',   tropT(:ncol),      ncol, lchnk)
         599     1489176 :     call outfld('TROPF_Z',   tropZ(:ncol),      ncol, lchnk)
         600   543653136 :     call outfld('TROPF_DZ',  tropDZ(:ncol, :), ncol, lchnk)
         601   543653136 :     call outfld('TROPF_PD',  tropPdf(:ncol, :), ncol, lchnk)
         602     1489176 :     call outfld('TROPF_FD',  tropFound(:ncol),  ncol, lchnk)
         603             :
         604             :
         605             :     ! If requested, do all of the algorithms.
         606             :     if (output_all) then
         607             :
         608             :       do alg = 2, TROP_NALG
         609             :
         610             :         ! Find the tropopause using just the analytic algorithm.
         611             :         call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=alg, backup=TROP_ALG_NONE)
         612             :
         613             :         tropPdf(:,:) = 0._r8
         614             :         tropFound(:) = 0._r8
         615             :
         616             :         do i = 1, ncol
         617             :           if (tropLev(i) /= NOTFOUND) then
         618             :             tropPdf(i, tropLev(i)) = 1._r8
         619             :             tropFound(i) = 1._r8
         620             :           end if
         621             :         end do
         622             :
         623             :         call outfld('TROP' // TROP_LETTER(alg) // '_P',   tropP(:ncol),      ncol, lchnk)
         624             :         call outfld('TROP' // TROP_LETTER(alg) // '_T',   tropT(:ncol),      ncol, lchnk)
         625             :         call outfld('TROP' // TROP_LETTER(alg) // '_Z',   tropZ(:ncol),      ncol, lchnk)
         626             :         call outfld('TROP' // TROP_LETTER(alg) // '_PD',  tropPdf(:ncol, :), ncol, lchnk)
         627             :         call outfld('TROP' // TROP_LETTER(alg) // '_FD',  tropFound(:ncol),  ncol, lchnk)
         628             :       end do
         629             :     end if
         630     1489176 :   end subroutine tropopause_output
         631             : end module tropopause

Generated by: LCOV version 1.14