chem/module_tropopause.F

   1 ! This module is used to diagnose the location of the tropopause.
   2 ! Two routines are specifed, a primary routine (twmo) which is tried first
   3 ! and a backup routine (climate) which will be tried only if the first routine fails.
   4 ! If the tropopause can not be identified by either routine, then a NOTFOUND is returned.
   5 !
   6 ! These routines are based upon code in the WACCM chemistry module
   7 ! including mo_tropoause.F90 and llnl_set_chem_trop.F90. The code
   8 ! for the Reichler et al. [2003] algorithm is from:
   9 !
  10 !   http://www.gfdl.noaa.gov/~tjr/TROPO/tropocode.htm
  11 !
  12 ! Author: Jeff Lee
  13 ! Created: Feb., 2011
  14
  15  module module_tropopause
  16
  17       implicit none
  18
  19       private
  20
  21       public  :: tropopause_init
  22       public  :: tropopause_driver
  23
  24       save
  25
  26       integer,  parameter :: r8    = selected_real_kind(12)
  27
  28       real(r8), parameter :: pi    =3.14159265358979323846_r8
  29       real(r8), parameter :: mb2Pa = 100._r8
  30       real(r8), parameter :: d2r   = pi/180._r8
  31       real(r8), parameter :: zero  = 0._r8
  32       real(r8), parameter :: twopi = pi*2._r8
  33
  34       integer,  parameter :: NOTFOUND = -1
  35
  36       integer :: iend
  37       integer :: jend
  38
  39       integer :: tropo_month_n      ! # of month in tropopause file
  40       integer :: tropo_lon_n        ! # of lon   in tropopause file
  41       integer :: tropo_lat_n        ! # of lat   in tropopause file
  42
  43       type tropo_type
  44           real(r8), pointer :: tropo_p_loc(:,:,:)  ! climatological tropopause pressures(Pa)
  45                                                    ! on wrf model grids
  46           logical           :: is_allocated
  47       end type tropo_type
  48
  49       type(tropo_type), allocatable :: tropo_bc(:)
  50
  51       ! physical constants
  52       ! These constants are set in module variables rather than as parameters
  53       ! to support the aquaplanet mode in which the constants have values determined
  54       ! by the experiment protocol
  55
  56       real(r8),parameter :: CONST_BOLTZ  = 1.38065e-23_r8  ! Boltzmann's constant ~ J/K/molecule
  57       real(r8),parameter :: CONST_AVOGAD = 6.02214e26_r8   ! Avogadro's number ~ molecules/kmole
  58       real(r8),parameter :: CONST_RGAS   = CONST_AVOGAD*CONST_BOLTZ ! Universal gas constant ~ J/K/kmole
  59       real(r8),parameter :: CONST_MWDAIR = 28.966_r8       ! molecular weight dry air ~ kg/kmole
  60       real(r8),parameter :: CONST_RDAIR  = CONST_RGAS/CONST_MWDAIR ! Dry air gas constant ~ J/K/kg
  61       real(r8),parameter :: CONST_CPDAIR  = 1.00464e3_r8    ! specific heat of dry air   ~ J/kg/K
  62
  63       real(r8),parameter :: gravit       = 9.80616_r8
  64       real(r8),parameter :: rair         = CONST_RDAIR
  65       real(r8),parameter :: cappa        = CONST_RDAIR/CONST_CPDAIR   !R/CP
  66
  67       real(r8),parameter :: cnst_kap    = cappa
  68       real(r8),parameter :: cnst_faktor = -gravit/rair
  69       real(r8),parameter :: cnst_ka1    = cnst_kap - 1._r8
  70
  71       contains
  72
  73       subroutine tropopause_driver( id, dtstep, current_date_char,  &
  74                                     t_phy, p_phy, p8w, zmid, z8w,   &
  75                                     tropo_lev, tropo_p, tropo_z,    &
  76                                     ids, ide, jds, jde, kds, kde,   &
  77                                     ims, ime, jms, jme, kms, kme,   &
  78                                     its, ite, jts, jte, kts, kte    )
  79
  80       implicit none
  81
  82 !----------------------------------------------------
  83 !     input arguments
  84 !----------------------------------------------------
  85       integer, intent(in   ) :: id,                             &
  86                                 ids,ide, jds,jde, kds,kde,      &
  87                                 ims,ime, jms,jme, kms,kme,      &
  88                                 its,ite, jts,jte, kts,kte
  89
  90       real, intent(in  ) :: dtstep
  91
  92       real, dimension( ims:ime, kms:kme, jms:jme ),             &
  93                intent(in   ) :: t_phy,  & ! t at mid-level
  94                                 p_phy,  & ! p at mid_level (Pa)
  95                                 zmid,   & ! z at mid_level (meters)
  96                                 z8w,    & ! z at interface (meters)
  97                                 p8w       ! p at interface (Pa)
  98
  99       real, dimension( ims:ime, jms:jme ),             &
 100                intent(inout) :: tropo_p,  & ! tropopause pressure (Pa)
 101                                 tropo_z     ! tropopause height   (meters)
 102       integer, dimension( ims:ime, jms:jme ),          &
 103                intent(inout) :: tropo_lev   ! tropopause level
 104
 105       CHARACTER (LEN=256),intent(in) :: current_date_char
 106 !----------------------------------------------------
 107 !     local scalars
 108 !----------------------------------------------------
 109
 110       integer :: i, j, k
 111
 112 !----------------------------------------------------
 113 !----------------------------------------------------
 114 ! ... tile dimensions, needed for nest domains
 115 !----------------------------------------------------
 116       iend = min(ite,ide-1)
 117       jend = min(jte,jde-1)
 118
 119       tropo_lev(its:iend,jts:jend) = NOTFOUND
 120
 121 ! --  This is the primary routine
 122
 123       call tropopause_twmo (id, t_phy, p_phy, p8w, zmid, z8w, &
 124                             tropo_lev, tropo_p, tropo_z,      &
 125                             ids,ide, jds,jde, kds,kde,        &
 126                             ims,ime, jms,jme, kms,kme,        &
 127                             its,ite, jts,jte, kts,kte         )
 128
 129 ! --  This is the backup routine
 130
 131       if ( any(tropo_lev(its:iend,jts:jend) == NOTFOUND) ) then
 132
 133          call tropopause_climate (id, current_date_char,       &
 134                                   p_phy, p8w, zmid, z8w,       &
 135                                   tropo_lev, tropo_p, tropo_z, &
 136                                   ids,ide, jds,jde, kds,kde,   &
 137                                   ims,ime, jms,jme, kms,kme,   &
 138                                   its,ite, jts,jte, kts,kte    )
 139       end if
 140
 141    end subroutine tropopause_driver
 142
 143 !===========================================================================
 144
 145       subroutine tropopause_init (id, xlat, xlon, config_flags, &
 146                                   ids,ide, jds,jde, kds,kde,    &
 147                                   ims,ime, jms,jme, kms,kme,    &
 148                                   its,ite, jts,jte, kts,kte     )
 149 !---------------------------------------------------------------------
 150 !       ... new initialization routine for tropopause
 151 !---------------------------------------------------------------------
 152       use module_interpolate,  only : lininterp_init, lininterp, interp_type, lininterp_finish
 153       use module_configure,    only : grid_config_rec_type
 154
 155       implicit none
 156
 157 !---------------------------------------------------------------------
 158 ! ... dummy arguments
 159 !---------------------------------------------------------------------
 160       integer, intent(in) :: id,                          &
 161                              ids,ide, jds,jde, kds,kde,   &
 162                              ims,ime, jms,jme, kms,kme,   &
 163                              its,ite, jts,jte, kts,kte
 164
 165       real,    intent(in) :: xlat(ims:ime, jms:jme)
 166       real,    intent(in) :: xlon(ims:ime, jms:jme)
 167       type(grid_config_rec_type), intent(in) :: config_flags
 168
 169 !---------------------------------------------------------------------
 170 ! ... local variables
 171 !---------------------------------------------------------------------
 172
 173       type(interp_type) :: lon_wgts, lat_wgts
 174
 175       integer :: max_dom
 176       integer :: astat
 177       integer :: ncid
 178       integer :: varid
 179       integer :: dimid(3)
 180       integer :: start(3)
 181       integer :: count(3)
 182       integer :: dimid_lon
 183       integer :: dimid_lat
 184       integer :: dimid_month
 185       integer :: ndims
 186
 187       character(len=128) :: err_msg
 188       character(len=64)  :: filename
 189       character(len=3)   :: id_num
 190
 191       character(len=80) :: attribute
 192
 193       real(r8), allocatable :: tropo_p_in(:,:,:)  ! values of pressure levels (Pa) in tropopause file
 194       real(r8), allocatable :: tropo_lat(:)       ! values of lat (-90~90)in tropopause file
 195       real(r8), allocatable :: tropo_lon(:)       ! values of lon (0~360) in tropopause file
 196
 197       real(r8) :: wrf_lon(1)       ! to match kind, values of lon (0~360)
 198       real(r8) :: wrf_lat(1)       ! input to lininterp_init needs to be an array, not scalar
 199       real(r8) :: tmp_out(1)
 200
 201       integer  :: i, j, m
 202       CHARACTER(LEN=132) :: message
 203
 204       LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
 205
 206 include 'netcdf.inc'
 207 !---------------------------------------------------------------------
 208 ! ... tile dimensions
 209 !---------------------------------------------------------------------
 210       iend = min(ite,ide-1)
 211       jend = min(jte,jde-1)
 212
 213 !---------------------------------------------------------------------
 214 ! ... allocate tropo_bc
 215 !---------------------------------------------------------------------
 216    if( id == 1 .and. .not. allocated(tropo_bc) ) then
 217        CALL nl_get_max_dom( 1,max_dom )
 218
 219        allocate(tropo_bc(max_dom),stat=astat )
 220        if( astat /= 0 ) then
 221           CALL wrf_message( 'tropopause_init: failed to allocate tropo_bc' )
 222           CALL wrf_abort
 223        end if
 224
 225        tropo_bc(:)%is_allocated = .false.
 226   endif
 227
 228 tropo_bc_allocated : &
 229   if( .not. tropo_bc(id)%is_allocated ) then
 230
 231 !---------------------------------------------------------------------
 232 !... allocate tropo_bc type
 233 !--------------------------------------------------------------------
 234
 235 master_proc : &
 236    IF( wrf_dm_on_monitor() ) THEN
 237       write(id_num,'(i3)') 100+id
 238       write(message,*) 'tropopause_init: intializing domain ' // id_num(2:3)
 239       call wrf_message( trim(message) )
 240
 241 !---------------------------------------------------------------------
 242 ! ... open climate tropopause netcdf file
 243 !---------------------------------------------------------------------
 244 !     filename = 'clim_p_trop.nc'
 245       filename = config_flags%trop_lev_inname
 246       if( filename == ' ' ) then
 247         call wrf_message( 'tropopause_init: input filename not specified in namelist' )
 248         call wrf_abort
 249       endif
 250
 251       err_msg = 'tropopause_init: failed to open file ' // trim(filename)
 252       call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
 253       write(message,*) 'tropopause_init: open filename= ',filename
 254       call wrf_message( trim(message) )
 255 !---------------------------------------------------------------------
 256 ! ... get dimensions
 257 !---------------------------------------------------------------------
 258       err_msg = 'tropopause_init: failed to get time id'
 259       call handle_ncerr( nf_inq_dimid( ncid, 'time', dimid_month ), trim(err_msg) )
 260       err_msg = 'tropopause_init: failed to get time'
 261       call handle_ncerr( nf_inq_dimlen( ncid, dimid_month, tropo_month_n ), trim(err_msg) )
 262       if( tropo_month_n /= 12 )then
 263        write(message,*) 'tropopause_init: number of months = ',tropo_month_n,'; expecting 12'
 264        call wrf_message( trim(message) )
 265        call wrf_abort
 266       end if
 267
 268       err_msg = 'tropopause_init: failed to get lat id'
 269       call handle_ncerr( nf_inq_dimid( ncid, 'lat', dimid_lat ), trim(err_msg) )
 270       err_msg = 'tropopause_init: failed to get lat'
 271       call handle_ncerr( nf_inq_dimlen( ncid, dimid_lat, tropo_lat_n ), trim(err_msg) )
 272
 273       err_msg = 'tropopause_init: failed to get lon id'
 274       call handle_ncerr( nf_inq_dimid( ncid, 'lon', dimid_lon ), trim(err_msg) )
 275       err_msg = 'tropopause_init: failed to get lon'
 276       call handle_ncerr( nf_inq_dimlen( ncid, dimid_lon, tropo_lon_n ), trim(err_msg) )
 277
 278    END IF master_proc
 279 !---------------------------------------------------------------------
 280 ! ... bcast the dimensions
 281 !---------------------------------------------------------------------
 282 #ifdef DM_PARALLEL
 283       CALL wrf_dm_bcast_integer ( tropo_month_n , 1 )
 284       CALL wrf_dm_bcast_integer ( tropo_lat_n   , 1 )
 285       CALL wrf_dm_bcast_integer ( tropo_lon_n   , 1 )
 286 #endif
 287
 288 !---------------------------------------------------------------------
 289 ! ... allocate local arrays
 290 !---------------------------------------------------------------------
 291
 292       allocate( tropo_lat(tropo_lat_n), stat=astat )
 293       if( astat /= 0 ) then
 294          call wrf_message( 'tropopause_init: failed to allocate tropo_lat' )
 295          call wrf_abort
 296       end if
 297
 298       allocate( tropo_lon(tropo_lon_n), stat=astat )
 299       if( astat /= 0 ) then
 300          call wrf_message( 'tropopause_init: failed to allocate tropo_lon' )
 301          call wrf_abort
 302       end if
 303
 304       allocate( tropo_p_in(tropo_lon_n,tropo_lat_n,tropo_month_n), stat=astat )
 305       if( astat /= 0 ) then
 306          call wrf_message( 'tropopause_init: failed to allocate tropo_p_in' )
 307          call wrf_abort
 308       end if
 309
 310 !---------------------------------------------------------------------
 311 ! ... allocate tropo_bc(id) component arrays
 312 !---------------------------------------------------------------------
 313
 314       allocate( tropo_bc(id)%tropo_p_loc(its:iend,jts:jend,tropo_month_n), stat=astat )
 315       if( astat /= 0 ) then
 316          call wrf_message( 'tropopause_init: failed to allocate tropo_bc(id)%tropo_p_loc' )
 317          call wrf_abort
 318       end if
 319
 320       tropo_bc(id)%is_allocated = .true.
 321
 322 !---------------------------------------------------------------------
 323 ! ... read arrays
 324 !---------------------------------------------------------------------
 325 master_proc_a : &
 326    IF ( wrf_dm_on_monitor() ) THEN
 327 !---------------------------
 328 !lat
 329       err_msg = 'tropopause_init: failed to get lat variable id'
 330       call handle_ncerr( nf_inq_varid( ncid, 'lat', varid ), trim(err_msg) )
 331       err_msg = 'tropopause_init: failed to read lat variable'
 332       call handle_ncerr( nf_get_var_double( ncid, varid, tropo_lat ), trim(err_msg) )
 333
 334       !-------- check unit
 335      tropo_lat(:) = tropo_lat(:) * d2r
 336
 337 !---------------------------
 338 !lon
 339       err_msg = 'tropopause_init: failed to get lon variable id'
 340       call handle_ncerr( nf_inq_varid( ncid, 'lon', varid ), trim(err_msg) )
 341       err_msg = 'tropopause_init: failed to read lon variable'
 342       call handle_ncerr( nf_get_var_double( ncid, varid, tropo_lon ), trim(err_msg) )
 343
 344       !-------- check unit and note: 0-360 degree
 345       tropo_lon(:) = tropo_lon(:) * d2r
 346 !---------------------------
 347 !tropo_p_in
 348       err_msg = 'tropopause_init: failed to get trop_p variable id'
 349       call handle_ncerr( nf_inq_varid( ncid, 'trop_p', varid ), trim(err_msg) )
 350
 351       ! check dimensions
 352
 353       err_msg = 'tropopause_init: failed to get ndims of trop_p variable'
 354       call handle_ncerr( nf_inq_varndims( ncid, varid, ndims ), trim(err_msg) )
 355
 356       if( ndims /= 3 ) then
 357          write(message,*) 'tropopause_init: error! variable trop_p has ndims = ',ndims,', expecting 3'
 358          call wrf_message( trim(message) )
 359          call wrf_abort
 360       end if
 361
 362       err_msg = 'tropopause_init: failed to get dimid of vmr variable'
 363       call handle_ncerr( nf_inq_vardimid( ncid, varid, dimid ), trim(err_msg) )
 364
 365       if( dimid(1) /= dimid_lon   .or. dimid(2) /= dimid_lat .or. &
 366           dimid(3) /= dimid_month )then
 367           write(message,*) 'tropopause_init: error! dimensions in wrong order for variable trop_p,'// &
 368                'expecting (lon,lat,month)'
 369           call wrf_message( trim(message) )
 370           call wrf_abort
 371       end if
 372
 373       !------------------------------------------------------------------
 374       ! ... read in the tropo_p_in values
 375       !------------------------------------------------------------------
 376       err_msg = 'tropopause_init: failed to read trop_p variable'
 377       call handle_ncerr( nf_get_var_double( ncid, varid, tropo_p_in ), trim(err_msg) )
 378
 379 !---------------------------------------------------------------------
 380 !... close input netcdf file
 381 !---------------------------------------------------------------------
 382       err_msg = 'tropopause_init: failed to close file ' // trim(filename)
 383       call handle_ncerr( nf_close( ncid ), trim(err_msg) )
 384
 385    END IF master_proc_a
 386
 387 !---------------------------------------------------------------------
 388 ! ... bcast the variables
 389 !---------------------------------------------------------------------
 390 #ifdef DM_PARALLEL
 391       CALL wrf_dm_bcast_double ( tropo_lat, size(tropo_lat) )
 392       CALL wrf_dm_bcast_double ( tropo_lon, size(tropo_lon) )
 393       CALL wrf_dm_bcast_double ( tropo_p_in, size(tropo_p_in) )
 394 #endif
 395
 396 !------------------------------------------------------------------
 397 ! ... linear interpolation from tropo_p_in to tropo_p_loc
 398 !------------------------------------------------------------------
 399       !-------------------------------------
 400       !... get wrf_lat_1d and wrf_lon_1d
 401       !... note: grid%XLAT(ide,:)  =0.
 402       !... note: grid%XLONG(:,jde) =0.
 403       !...       => iend = min(ite,ide-1)
 404       !             jend = min(jte,jde-1)
 405       !-------------------------------------
 406
 407       !-------------------------------------------
 408       !... for every point in the tile
 409       !-------------------------------------------
 410       do i = its,iend
 411         do j = jts,jend   ! input to lininterp_init needs to be an array, not scalar
 412          !-------------------------------------------
 413          !...  from degrees to radians
 414          !-------------------------------------------
 415           wrf_lat(1) = xlat(i,j)*d2r
 416
 417           wrf_lon(1) = xlon(i,j)
 418           if( wrf_lon(1)  < 0.0_r8 ) wrf_lon(1) = wrf_lon(1) + 360.0_r8
 419           wrf_lon(1) = wrf_lon(1)*d2r
 420
 421          !-------------------------------------
 422          !...initialization interp routine
 423          !-------------------------------------
 424           call lininterp_init( tropo_lon, tropo_lon_n, wrf_lon, 1, 2, lon_wgts, zero, twopi )
 425           call lininterp_init( tropo_lat, tropo_lat_n, wrf_lat, 1, 1, lat_wgts )
 426
 427          !-------------------------------------
 428          !... linear interpolation
 429          !-------------------------------------
 430           do m = 1,tropo_month_n
 431             call lininterp( tropo_p_in(:,:,m), tropo_lon_n, tropo_lat_n, &
 432                             tmp_out, 1 , lon_wgts, lat_wgts)
 433             tropo_bc(id)%tropo_p_loc(i,j,m) = tmp_out(1)
 434           end do
 435
 436         end do
 437       end do
 438
 439       call lininterp_finish( lon_wgts )
 440       call lininterp_finish( lat_wgts )
 441
 442       deallocate(tropo_lon)
 443       deallocate(tropo_lat)
 444       deallocate(tropo_p_in)
 445
 446       call wrf_message( ' ' )
 447       write(message,*) 'tropopause_init: DONE intialized domain ',id
 448       call wrf_message( trim(message) )
 449       call wrf_message( ' ' )
 450
 451 endif tropo_bc_allocated
 452
 453       end subroutine tropopause_init
 454
 455 !===========================================================================
 456
 457       subroutine handle_ncerr( ret, mes )
 458 !---------------------------------------------------------------------
 459 !       ... netcdf error handling routine
 460 !---------------------------------------------------------------------
 461
 462       implicit none
 463
 464 !---------------------------------------------------------------------
 465 !       ... dummy arguments
 466 !---------------------------------------------------------------------
 467       integer, intent(in) :: ret
 468       character(len=*), intent(in) :: mes
 469
 470 include 'netcdf.inc'
 471
 472       if( ret /= nf_noerr ) then
 473          call wrf_message( trim(mes) )
 474          call wrf_message( trim(nf_strerror(ret)) )
 475          call wrf_abort
 476       end if
 477
 478       end subroutine handle_ncerr
 479
 480 !===========================================================================
 481
 482   ! This routine uses an implementation of Reichler et al. [2003] done by
 483   ! Reichler and downloaded from his web site. This is similar to the WMO
 484   ! routines, but is designed for GCMs with a coarse vertical grid.
 485
 486   subroutine tropopause_twmo (id, t_phy, p_phy, p8w, zmid, z8w, &
 487                               tropo_lev, tropo_p,tropo_z,       &
 488                               ids,ide, jds,jde, kds,kde,        &
 489                               ims,ime, jms,jme, kms,kme,        &
 490                               its,ite, jts,jte, kts,kte         )
 491       implicit none
 492
 493 !----------------------------------------------------
 494 !     input arguments
 495 !----------------------------------------------------
 496       integer, intent(in) :: id,                         &
 497                              ids,ide, jds,jde, kds,kde,  &
 498                              ims,ime, jms,jme, kms,kme,  &
 499                              its,ite, jts,jte, kts,kte
 500
 501       real, dimension( ims:ime , kms:kme , jms:jme ),    &
 502             intent(in)    :: t_phy,  & ! t at mid-level
 503                              p_phy,  & ! p at mid_level (Pa)
 504                              zmid,   & ! z at mid_level (meters)
 505                              z8w,    & ! z at interface (meters)
 506                              p8w       ! p at interface (Pa)
 507
 508       real, dimension( ims:ime, jms:jme ),             &
 509                intent(inout) :: tropo_p,  & ! tropopause pressure (Pa)
 510                                 tropo_z     ! tropopause height   (meters)
 511       integer, dimension( ims:ime, jms:jme ),          &
 512                intent(inout) :: tropo_lev   ! tropopause level
 513 !----------------------------------------------------
 514 !     local scalars
 515 !----------------------------------------------------
 516     real, dimension( kts:kte) :: t_temp    ! store t from top to bottom
 517     real, dimension( kts:kte) :: p_temp    ! store p from top to bottom
 518
 519     real(r8), parameter     :: gamma    = -0.002_r8    ! K/m
 520     real(r8), parameter     :: plimu    = 45000._r8    ! Pa
 521     real(r8), parameter     :: pliml    = 7500._r8     ! Pa
 522
 523     integer                 :: i
 524     integer                 :: j
 525     integer                 :: k
 526     integer                 :: kk
 527     integer                 :: pLev
 528
 529     real(r8)                :: tP                       ! tropopause pressure (Pa)
 530     real(r8)                :: dZdlogP
 531 !----------------------------------------------------
 532
 533     pLev = kte - kts + 1
 534
 535     ! Iterate over all of the grids.
 536     do i = its,iend
 537     do j = jts,jend
 538
 539        ! subroutine twmo expects variables from top to bottom
 540        ! t_phy and p_phy are from bottom to top
 541
 542        do k = kts,kte
 543           kk = pLev - k + 1
 544           t_temp(kk) = t_phy(i,k,j)
 545           p_temp(kk) = p_phy(i,k,j)
 546        end do
 547
 548        ! Use the routine from Reichler.
 549
 550        call twmo(pLev, t_temp, p_temp, plimu, pliml, gamma, tP)
 551
 552        ! if successful, store the results and find the level and temperature.
 553        if (tP > 0) then
 554           ! Find the associated level, from bottom to top
 555           do k = kts,kte-1
 556              if (tP >= p8w(i,k,j)) then
 557                 tropo_lev(i,j) = k - 1      ! needed for wrf, from bottom to top
 558                 tropo_p  (i,j) = tP
 559                 exit
 560              end if
 561           end do
 562           !----------------------------------------------------------
 563           ! get tropopause height
 564           ! Intrepolate the geopotential height linearly against log(P)
 565
 566           k = k - 1     ! needed for wrf, from bottom to top
 567
 568           ! Is the tropoause at the midpoint?
 569           if (tropo_p(i,j) == p_phy(i,k,j)) then
 570              tropo_z(i,j)= zmid(i,k,j)
 571           else if (tropo_p(i,j) < p_phy(i,k,j)) then
 572              ! It is below the midpoint. Make sure we aren't at the bottom.
 573              if ( k > kts ) then
 574                 dZdlogP = (zmid(i,k,j) - z8w(i,k-1,j)) / &
 575                           (log(p_phy(i,k,j)) - log(p8w(i,k-1,j)) )
 576                 tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP
 577              end if
 578           else
 579              ! It is above the midpoint? Make sure we aren't at the top.
 580              if ( k < kte ) then
 581                 dZdlogP = (zmid(i,k,j) - z8w(i,k,j)) / &
 582                           (log(p_phy(i,k,j)) - log(p8w(i,k,j)) )
 583                 tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP
 584              end if
 585           end if
 586           !----------------------------------------------------------
 587        end if
 588
 589 !      if ( tropo_lev(i,j) == NOTFOUND ) then
 590 !         write (*,*) 'tropopause_twmo: NOTFOUND at id, i, j = ', id, i, j
 591 !      end if
 592
 593     end do
 594     end do
 595
 596   end subroutine tropopause_twmo
 597
 598 !===========================================================================
 599
 600   ! This routine is an implementation of Reichler et al. [2003] done by
 601   ! Reichler and downloaded from his web site. Minimal modifications were
 602   ! made to have the routine work within the CAM framework (i.e. using
 603   ! CAM constants and types).
 604   !
 605   ! NOTE: I am not a big fan of the goto's and multiple returns in this
 606   ! code, but for the moment I have left them to preserve as much of the
 607   ! original and presumably well tested code as possible.
 608   !
 609   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 610   !
 611   ! determination of tropopause height from gridded temperature data
 612   !
 613   ! reference: Reichler, T., M. Dameris, and R. Sausen (2003)
 614   !
 615   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 616
 617   subroutine twmo(level, t, p, plimu, pliml, gamma, trp)
 618
 619     integer, intent(in)                     :: level
 620 !jeff
 621 !   real(r8), intent(in), dimension(level)  :: t, p
 622     real,     intent(in), dimension(level)  :: t, p
 623
 624     real(r8), intent(in)                    :: plimu, pliml, gamma
 625     real(r8), intent(out)                   :: trp
 626
 627     real(r8), parameter                     :: deltaz = 2000.0_r8
 628
 629     real(r8)                                :: pmk, pm, a, b, tm, dtdp, dtdz
 630     real(r8)                                :: ag, bg, ptph
 631     real(r8)                                :: pm0, pmk0, dtdz0
 632     real(r8)                                :: p2km, asum, aquer
 633     real(r8)                                :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
 634     integer                                 :: icount, jj
 635     integer                                 :: j
 636
 637
 638     trp=-99.0_r8                           ! negative means not valid
 639
 640 lev_loop : &
 641     do j=level,2,-1
 642
 643       ! dt/dz
 644       pmk= .5_r8 * (p(j-1)**cnst_kap+p(j)**cnst_kap)
 645       pm = pmk**(1/cnst_kap)
 646       a = (t(j-1)-t(j))/(p(j-1)**cnst_kap-p(j)**cnst_kap)
 647       b = t(j)-(a*p(j)**cnst_kap)
 648       tm = a * pmk + b
 649       dtdp = a * cnst_kap * (pm**cnst_ka1)
 650       dtdz = cnst_faktor*dtdp*pm/tm
 651       ! dt/dz valid?
 652 !     if (j.eq.level)    go to 999     ! no, start level, initialize first
 653 !     if (dtdz.le.gamma) go to 999     ! no, dt/dz < -2 K/km
 654 !     if (pm.gt.plimu)   go to 999     ! no, too low
 655       if( j == level .or. dtdz <= gamma .or. pm > plimu ) then
 656         pm0   = pm
 657         pmk0  = pmk
 658         dtdz0 = dtdz
 659         cycle lev_loop
 660       endif
 661
 662       ! dtdz is valid, calculate tropopause pressure
 663       if (dtdz0 < gamma) then
 664         ag = (dtdz-dtdz0) / (pmk-pmk0)
 665         bg = dtdz0 - (ag * pmk0)
 666         ptph = exp(log((gamma-bg)/ag)/cnst_kap)
 667       else
 668         ptph = pm
 669       endif
 670
 671 !     if (ptph.lt.pliml) go to 999
 672 !     if (ptph.gt.plimu) go to 999
 673       if( ptph < pliml .or. ptph > plimu ) then
 674         pm0   = pm
 675         pmk0  = pmk
 676         dtdz0 = dtdz
 677         cycle lev_loop
 678       endif
 679
 680       ! 2nd test: dtdz above 2 km must not exceed gamma
 681       p2km = ptph + deltaz*(pm/tm)*cnst_faktor     ! p at ptph + 2km
 682       asum = 0.0_r8                                ! dtdz above
 683       icount = 0                                   ! number of levels above
 684
 685       ! test until apm < p2km
 686 lev_loop_a : &
 687       do jj=j,2,-1
 688
 689         pmk2 = .5_r8 * (p(jj-1)**cnst_kap+p(jj)**cnst_kap) ! p mean ^kappa
 690         pm2 = pmk2**(1/cnst_kap)                           ! p mean
 691         if( pm2 > ptph ) then                     ! doesn't happen
 692           cycle lev_loop_a
 693         endif
 694         if( pm2 < p2km ) then                     ! ptropo is valid
 695           trp = ptph
 696           exit lev_loop
 697         endif
 698
 699         a2 = (t(jj-1)-t(jj))                     ! a
 700         a2 = a2/(p(jj-1)**cnst_kap-p(jj)**cnst_kap)
 701         b2 = t(jj)-(a2*p(jj)**cnst_kap)          ! b
 702         tm2 = a2 * pmk2 + b2                     ! T mean
 703         dtdp2 = a2 * cnst_kap * (pm2**(cnst_kap-1))  ! dt/dp
 704         dtdz2 = cnst_faktor*dtdp2*pm2/tm2
 705         asum = asum+dtdz2
 706         icount = icount+1
 707         aquer = asum/float(icount)               ! dt/dz mean
 708
 709         ! discard ptropo ?
 710         if( aquer <= gamma ) then                ! dt/dz above < gamma
 711           pm0   = pm
 712           pmk0  = pmk
 713           dtdz0 = dtdz
 714           exit lev_loop_a
 715         endif
 716
 717       enddo lev_loop_a
 718
 719     enddo lev_loop
 720
 721   end subroutine twmo
 722
 723 !===========================================================================
 724
 725   ! Read the tropopause pressure in from a file containging a climatology. The
 726   ! data is interpolated to the current dat of year and latitude.
 727   !
 728   ! NOTE: The data is read in during tropopause_init and stored in the module
 729   ! variable tropo_bc(id)%tropo_p_loc
 730
 731   subroutine tropopause_climate (id, current_date_char,        &
 732                                  p_phy, p8w, zmid, z8w,        &
 733                                  tropo_lev, tropo_p,  tropo_z, &
 734                                  ids,ide, jds,jde, kds,kde,    &
 735                                  ims,ime, jms,jme, kms,kme,    &
 736                                  its,ite, jts,jte, kts,kte     )
 737       implicit none
 738
 739 !----------------------------------------------------
 740 !     input arguments
 741 !----------------------------------------------------
 742       integer, intent(in) :: id,                         &
 743                              ids,ide, jds,jde, kds,kde,  &
 744                              ims,ime, jms,jme, kms,kme,  &
 745                              its,ite, jts,jte, kts,kte
 746
 747       real, dimension( ims:ime , kms:kme , jms:jme ),    &
 748             intent(in)    :: p_phy,  & ! p at mid-level (Pa)
 749                              zmid,   & ! z at mid_level (meters)
 750                              z8w,    & ! z at interface (meters)
 751                              p8w       ! p at interface (Pa)
 752
 753       real, dimension( ims:ime, jms:jme ),             &
 754                intent(inout) :: tropo_p,  & ! tropopause pressure (Pa)
 755                                 tropo_z     ! tropopause height   (meters)
 756       integer, dimension( ims:ime, jms:jme ),          &
 757                intent(inout) :: tropo_lev   ! tropopause level
 758
 759       CHARACTER (LEN=256),intent(in) :: current_date_char
 760 !----------------------------------------------------
 761
 762     ! Local Variables
 763
 764     real      :: del_time
 765     integer   :: last
 766     integer   :: next
 767
 768     real(r8)  :: tP                       ! tropopause pressure (Pa)
 769     real(r8)  :: dZdlogP
 770
 771     integer   :: i
 772     integer   :: j
 773     integer   :: k
 774     CHARACTER(LEN=132) :: message
 775 !-----------------------------------------------------------------------
 776
 777     if (any(tropo_lev(its:iend,jts:jend) == NOTFOUND)) then
 778
 779        !------------------------------------------------------------------------
 780        ! setup time interpolation
 781        !------------------------------------------------------------------------
 782
 783        call get_time_interp_factors( current_date_char, last, next, del_time )
 784
 785        ! Iterate over all of the grids.
 786
 787        do i = its,iend
 788        do j = jts,jend
 789
 790           ! Skip grids in which the tropopause has already been found.
 791
 792           if (tropo_lev(i,j) == NOTFOUND) then
 793
 794              !--------------------------------------------------------
 795              ! ... get tropopause level from climatology
 796              !--------------------------------------------------------
 797              ! Interpolate the tropopause pressure.
 798              tP =  tropo_bc(id)%tropo_p_loc(i,j,last) &
 799                 + (tropo_bc(id)%tropo_p_loc(i,j,next) &
 800                 -  tropo_bc(id)%tropo_p_loc(i,j,last))*del_time
 801
 802              if (tP > 0) then
 803                 ! Find the associated level, from bottom to top
 804                 do k = kts,kte-1
 805                    if (tP >= p8w(i,k,j)) then
 806                       tropo_lev(i,j) = k - 1   ! needed for wrf, from bottom to top
 807                       tropo_p  (i,j) = tP
 808                       exit
 809                    end if
 810                 end do
 811
 812                 !----------------------------------------------------------
 813                 ! get tropopause height
 814                 ! Intrepolate the geopotential height linearly against log(P)
 815
 816                 k = k - 1     ! needed for wrf
 817
 818                 ! Is the tropoause at the midpoint?
 819                 if (tropo_p(i,j) == p_phy(i,k,j)) then
 820                    tropo_z(i,j)= zmid(i,k,j)
 821
 822                 else if (tropo_p(i,j) < p_phy(i,k,j)) then
 823                    ! It is below the midpoint. Make sure we aren't at the bottom.
 824                    if ( k > kts ) then
 825                       dZdlogP = (zmid(i,k,j) - z8w(i,k-1,j)) / &
 826                                 (log(p_phy(i,k,j)) - log(p8w(i,k-1,j)) )
 827                       tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP
 828                    end if
 829                 else
 830                    ! It is above the midpoint? Make sure we aren't at the top.
 831                    if ( k < kte ) then
 832                       dZdlogP = (zmid(i,k,j) - z8w(i,k,j)) / &
 833                                 (log(p_phy(i,k,j)) - log(p8w(i,k,j)) )
 834                       tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP
 835                    end if
 836                 end if
 837                 !----------------------------------------------------------
 838              end if
 839           end if
 840        end do
 841        end do
 842     end if
 843
 844     if (any(tropo_lev(its:iend,jts:jend) == NOTFOUND)) then
 845        write(message,*) 'tropopause_climate: Warning: some tropopause levels still NOTFOUND'
 846     else
 847        write(message,*) 'tropopause_climate: Warning: Done finding tropopause'
 848     end if
 849     call wrf_message( trim(message) )
 850
 851   end subroutine tropopause_climate
 852
 853 !===========================================================================
 854
 855       subroutine get_time_interp_factors(current_date_char , last, next, del_time )
 856 !-----------------------------------------------------------------------------
 857 ! ... setup the time interpolation
 858 !-----------------------------------------------------------------------------
 859
 860       use module_date_time, only : get_julgmt
 861
 862       implicit none
 863
 864 !-----------------------------------------------------------------------------
 865 ! ... dummy arguments
 866 !-----------------------------------------------------------------------------
 867       CHARACTER (LEN=256),intent(in) :: current_date_char
 868
 869       integer, intent(out) :: next, last
 870       real,    intent(out) :: del_time
 871 !-----------------------------------------------------------------------------
 872 ! ... local variables
 873 !-----------------------------------------------------------------------------
 874
 875       integer, parameter :: day_of_year(12) = (/ 16, 45, 75, 105, 136, 166, 197, &
 876                                                  228, 258, 289, 319, 350        /)
 877       integer :: julyr , julday
 878       real    :: gmt
 879       real    :: calday
 880
 881       integer  :: m
 882 !---------------------------------------------------------------------------------
 883
 884       call get_julgmt ( current_date_char, julyr, julday, gmt )
 885
 886       calday = real(julday) + gmt
 887
 888       if( calday < day_of_year(1) ) then
 889          next = 1
 890          last = 12
 891          del_time = (365. + calday - day_of_year(12)) &
 892                 / (365. + day_of_year(1) - day_of_year(12))
 893       else if( calday >= day_of_year(12) ) then
 894          next = 1
 895          last = 12
 896          del_time = (calday - day_of_year(12)) &
 897                 / (365. + day_of_year(1) - day_of_year(12))
 898       else
 899          do m = 11,1,-1
 900             if( calday >= day_of_year(m) ) then
 901                exit
 902             end if
 903          end do
 904          last = m
 905          next = m + 1
 906          del_time = (calday - day_of_year(m)) / (day_of_year(m+1) - day_of_year(m))
 907       end if
 908
 909       del_time = max( min( 1.,del_time ),0. )
 910
 911       end subroutine get_time_interp_factors
 912
 913 !===========================================================================
 914       end module module_tropopause
 915