! This module is used to diagnose the location of the tropopause. ! Two routines are specifed, a primary routine (twmo) which is tried first ! and a backup routine (climate) which will be tried only if the first routine fails. ! If the tropopause can not be identified by either routine, then a NOTFOUND is returned. ! ! These routines are based upon code in the WACCM chemistry module ! including mo_tropoause.F90 and llnl_set_chem_trop.F90. The code ! for the Reichler et al. [2003] algorithm is from: ! ! http://www.gfdl.noaa.gov/~tjr/TROPO/tropocode.htm ! ! Author: Jeff Lee ! Created: Feb., 2011 module module_tropopause implicit none private public :: tropopause_init public :: tropopause_driver save integer, parameter :: r8 = selected_real_kind(12) real(r8), parameter :: pi =3.14159265358979323846_r8 real(r8), parameter :: mb2Pa = 100._r8 real(r8), parameter :: d2r = pi/180._r8 real(r8), parameter :: zero = 0._r8 real(r8), parameter :: twopi = pi*2._r8 integer, parameter :: NOTFOUND = -1 integer :: iend integer :: jend integer :: tropo_month_n ! # of month in tropopause file integer :: tropo_lon_n ! # of lon in tropopause file integer :: tropo_lat_n ! # of lat in tropopause file type tropo_type real(r8), pointer :: tropo_p_loc(:,:,:) ! climatological tropopause pressures(Pa) ! on wrf model grids logical :: is_allocated end type tropo_type type(tropo_type), allocatable :: tropo_bc(:) ! physical constants ! These constants are set in module variables rather than as parameters ! to support the aquaplanet mode in which the constants have values determined ! by the experiment protocol real(r8),parameter :: CONST_BOLTZ = 1.38065e-23_r8 ! Boltzmann's constant ~ J/K/molecule real(r8),parameter :: CONST_AVOGAD = 6.02214e26_r8 ! Avogadro's number ~ molecules/kmole real(r8),parameter :: CONST_RGAS = CONST_AVOGAD*CONST_BOLTZ ! Universal gas constant ~ J/K/kmole real(r8),parameter :: CONST_MWDAIR = 28.966_r8 ! molecular weight dry air ~ kg/kmole real(r8),parameter :: CONST_RDAIR = CONST_RGAS/CONST_MWDAIR ! Dry air gas constant ~ J/K/kg real(r8),parameter :: CONST_CPDAIR = 1.00464e3_r8 ! specific heat of dry air ~ J/kg/K real(r8),parameter :: gravit = 9.80616_r8 real(r8),parameter :: rair = CONST_RDAIR real(r8),parameter :: cappa = CONST_RDAIR/CONST_CPDAIR !R/CP real(r8),parameter :: cnst_kap = cappa real(r8),parameter :: cnst_faktor = -gravit/rair real(r8),parameter :: cnst_ka1 = cnst_kap - 1._r8 contains subroutine tropopause_driver( id, dtstep, current_date_char, & t_phy, p_phy, p8w, zmid, z8w, & tropo_lev, tropo_p, tropo_z, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) implicit none !---------------------------------------------------- ! input arguments !---------------------------------------------------- integer, intent(in ) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real, intent(in ) :: dtstep real, dimension( ims:ime, kms:kme, jms:jme ), & intent(in ) :: t_phy, & ! t at mid-level p_phy, & ! p at mid_level (Pa) zmid, & ! z at mid_level (meters) z8w, & ! z at interface (meters) p8w ! p at interface (Pa) real, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_p, & ! tropopause pressure (Pa) tropo_z ! tropopause height (meters) integer, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_lev ! tropopause level CHARACTER (LEN=256),intent(in) :: current_date_char !---------------------------------------------------- ! local scalars !---------------------------------------------------- integer :: i, j, k !---------------------------------------------------- !---------------------------------------------------- ! ... tile dimensions, needed for nest domains !---------------------------------------------------- iend = min(ite,ide-1) jend = min(jte,jde-1) tropo_lev(its:iend,jts:jend) = NOTFOUND ! -- This is the primary routine call tropopause_twmo (id, t_phy, p_phy, p8w, zmid, z8w, & tropo_lev, tropo_p, tropo_z, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) ! -- This is the backup routine if ( any(tropo_lev(its:iend,jts:jend) == NOTFOUND) ) then call tropopause_climate (id, current_date_char, & p_phy, p8w, zmid, z8w, & tropo_lev, tropo_p, tropo_z, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) end if end subroutine tropopause_driver !=========================================================================== subroutine tropopause_init (id, xlat, xlon, config_flags, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !--------------------------------------------------------------------- ! ... new initialization routine for tropopause !--------------------------------------------------------------------- use module_interpolate, only : lininterp_init, lininterp, interp_type, lininterp_finish use module_configure, only : grid_config_rec_type implicit none !--------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------- integer, intent(in) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real, intent(in) :: xlat(ims:ime, jms:jme) real, intent(in) :: xlon(ims:ime, jms:jme) type(grid_config_rec_type), intent(in) :: config_flags !--------------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------------- type(interp_type) :: lon_wgts, lat_wgts integer :: max_dom integer :: astat integer :: ncid integer :: varid integer :: dimid(3) integer :: start(3) integer :: count(3) integer :: dimid_lon integer :: dimid_lat integer :: dimid_month integer :: ndims character(len=128) :: err_msg character(len=64) :: filename character(len=3) :: id_num character(len=80) :: attribute real(r8), allocatable :: tropo_p_in(:,:,:) ! values of pressure levels (Pa) in tropopause file real(r8), allocatable :: tropo_lat(:) ! values of lat (-90~90)in tropopause file real(r8), allocatable :: tropo_lon(:) ! values of lon (0~360) in tropopause file real(r8) :: wrf_lon(1) ! to match kind, values of lon (0~360) real(r8) :: wrf_lat(1) ! input to lininterp_init needs to be an array, not scalar real(r8) :: tmp_out(1) integer :: i, j, m CHARACTER(LEN=132) :: message LOGICAL , EXTERNAL :: wrf_dm_on_monitor include 'netcdf.inc' !--------------------------------------------------------------------- ! ... tile dimensions !--------------------------------------------------------------------- iend = min(ite,ide-1) jend = min(jte,jde-1) !--------------------------------------------------------------------- ! ... allocate tropo_bc !--------------------------------------------------------------------- if( id == 1 .and. .not. allocated(tropo_bc) ) then CALL nl_get_max_dom( 1,max_dom ) allocate(tropo_bc(max_dom),stat=astat ) if( astat /= 0 ) then CALL wrf_message( 'tropopause_init: failed to allocate tropo_bc' ) CALL wrf_abort end if tropo_bc(:)%is_allocated = .false. endif tropo_bc_allocated : & if( .not. tropo_bc(id)%is_allocated ) then !--------------------------------------------------------------------- !... allocate tropo_bc type !-------------------------------------------------------------------- master_proc : & IF( wrf_dm_on_monitor() ) THEN write(id_num,'(i3)') 100+id write(message,*) 'tropopause_init: intializing domain ' // id_num(2:3) call wrf_message( trim(message) ) !--------------------------------------------------------------------- ! ... open climate tropopause netcdf file !--------------------------------------------------------------------- ! filename = 'clim_p_trop.nc' filename = config_flags%trop_lev_inname if( filename == ' ' ) then call wrf_message( 'tropopause_init: input filename not specified in namelist' ) call wrf_abort endif err_msg = 'tropopause_init: failed to open file ' // trim(filename) call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) ) write(message,*) 'tropopause_init: open filename= ',filename call wrf_message( trim(message) ) !--------------------------------------------------------------------- ! ... get dimensions !--------------------------------------------------------------------- err_msg = 'tropopause_init: failed to get time id' call handle_ncerr( nf_inq_dimid( ncid, 'time', dimid_month ), trim(err_msg) ) err_msg = 'tropopause_init: failed to get time' call handle_ncerr( nf_inq_dimlen( ncid, dimid_month, tropo_month_n ), trim(err_msg) ) if( tropo_month_n /= 12 )then write(message,*) 'tropopause_init: number of months = ',tropo_month_n,'; expecting 12' call wrf_message( trim(message) ) call wrf_abort end if err_msg = 'tropopause_init: failed to get lat id' call handle_ncerr( nf_inq_dimid( ncid, 'lat', dimid_lat ), trim(err_msg) ) err_msg = 'tropopause_init: failed to get lat' call handle_ncerr( nf_inq_dimlen( ncid, dimid_lat, tropo_lat_n ), trim(err_msg) ) err_msg = 'tropopause_init: failed to get lon id' call handle_ncerr( nf_inq_dimid( ncid, 'lon', dimid_lon ), trim(err_msg) ) err_msg = 'tropopause_init: failed to get lon' call handle_ncerr( nf_inq_dimlen( ncid, dimid_lon, tropo_lon_n ), trim(err_msg) ) END IF master_proc !--------------------------------------------------------------------- ! ... bcast the dimensions !--------------------------------------------------------------------- #ifdef DM_PARALLEL CALL wrf_dm_bcast_integer ( tropo_month_n , 1 ) CALL wrf_dm_bcast_integer ( tropo_lat_n , 1 ) CALL wrf_dm_bcast_integer ( tropo_lon_n , 1 ) #endif !--------------------------------------------------------------------- ! ... allocate local arrays !--------------------------------------------------------------------- allocate( tropo_lat(tropo_lat_n), stat=astat ) if( astat /= 0 ) then call wrf_message( 'tropopause_init: failed to allocate tropo_lat' ) call wrf_abort end if allocate( tropo_lon(tropo_lon_n), stat=astat ) if( astat /= 0 ) then call wrf_message( 'tropopause_init: failed to allocate tropo_lon' ) call wrf_abort end if allocate( tropo_p_in(tropo_lon_n,tropo_lat_n,tropo_month_n), stat=astat ) if( astat /= 0 ) then call wrf_message( 'tropopause_init: failed to allocate tropo_p_in' ) call wrf_abort end if !--------------------------------------------------------------------- ! ... allocate tropo_bc(id) component arrays !--------------------------------------------------------------------- allocate( tropo_bc(id)%tropo_p_loc(its:iend,jts:jend,tropo_month_n), stat=astat ) if( astat /= 0 ) then call wrf_message( 'tropopause_init: failed to allocate tropo_bc(id)%tropo_p_loc' ) call wrf_abort end if tropo_bc(id)%is_allocated = .true. !--------------------------------------------------------------------- ! ... read arrays !--------------------------------------------------------------------- master_proc_a : & IF ( wrf_dm_on_monitor() ) THEN !--------------------------- !lat err_msg = 'tropopause_init: failed to get lat variable id' call handle_ncerr( nf_inq_varid( ncid, 'lat', varid ), trim(err_msg) ) err_msg = 'tropopause_init: failed to read lat variable' call handle_ncerr( nf_get_var_double( ncid, varid, tropo_lat ), trim(err_msg) ) !-------- check unit tropo_lat(:) = tropo_lat(:) * d2r !--------------------------- !lon err_msg = 'tropopause_init: failed to get lon variable id' call handle_ncerr( nf_inq_varid( ncid, 'lon', varid ), trim(err_msg) ) err_msg = 'tropopause_init: failed to read lon variable' call handle_ncerr( nf_get_var_double( ncid, varid, tropo_lon ), trim(err_msg) ) !-------- check unit and note: 0-360 degree tropo_lon(:) = tropo_lon(:) * d2r !--------------------------- !tropo_p_in err_msg = 'tropopause_init: failed to get trop_p variable id' call handle_ncerr( nf_inq_varid( ncid, 'trop_p', varid ), trim(err_msg) ) ! check dimensions err_msg = 'tropopause_init: failed to get ndims of trop_p variable' call handle_ncerr( nf_inq_varndims( ncid, varid, ndims ), trim(err_msg) ) if( ndims /= 3 ) then write(message,*) 'tropopause_init: error! variable trop_p has ndims = ',ndims,', expecting 3' call wrf_message( trim(message) ) call wrf_abort end if err_msg = 'tropopause_init: failed to get dimid of vmr variable' call handle_ncerr( nf_inq_vardimid( ncid, varid, dimid ), trim(err_msg) ) if( dimid(1) /= dimid_lon .or. dimid(2) /= dimid_lat .or. & dimid(3) /= dimid_month )then write(message,*) 'tropopause_init: error! dimensions in wrong order for variable trop_p,'// & 'expecting (lon,lat,month)' call wrf_message( trim(message) ) call wrf_abort end if !------------------------------------------------------------------ ! ... read in the tropo_p_in values !------------------------------------------------------------------ err_msg = 'tropopause_init: failed to read trop_p variable' call handle_ncerr( nf_get_var_double( ncid, varid, tropo_p_in ), trim(err_msg) ) !--------------------------------------------------------------------- !... close input netcdf file !--------------------------------------------------------------------- err_msg = 'tropopause_init: failed to close file ' // trim(filename) call handle_ncerr( nf_close( ncid ), trim(err_msg) ) END IF master_proc_a !--------------------------------------------------------------------- ! ... bcast the variables !--------------------------------------------------------------------- #ifdef DM_PARALLEL CALL wrf_dm_bcast_double ( tropo_lat, size(tropo_lat) ) CALL wrf_dm_bcast_double ( tropo_lon, size(tropo_lon) ) CALL wrf_dm_bcast_double ( tropo_p_in, size(tropo_p_in) ) #endif !------------------------------------------------------------------ ! ... linear interpolation from tropo_p_in to tropo_p_loc !------------------------------------------------------------------ !------------------------------------- !... get wrf_lat_1d and wrf_lon_1d !... note: grid%XLAT(ide,:) =0. !... note: grid%XLONG(:,jde) =0. !... => iend = min(ite,ide-1) ! jend = min(jte,jde-1) !------------------------------------- !------------------------------------------- !... for every point in the tile !------------------------------------------- do i = its,iend do j = jts,jend ! input to lininterp_init needs to be an array, not scalar !------------------------------------------- !... from degrees to radians !------------------------------------------- wrf_lat(1) = xlat(i,j)*d2r wrf_lon(1) = xlon(i,j) if( wrf_lon(1) < 0.0_r8 ) wrf_lon(1) = wrf_lon(1) + 360.0_r8 wrf_lon(1) = wrf_lon(1)*d2r !------------------------------------- !...initialization interp routine !------------------------------------- call lininterp_init( tropo_lon, tropo_lon_n, wrf_lon, 1, 2, lon_wgts, zero, twopi ) call lininterp_init( tropo_lat, tropo_lat_n, wrf_lat, 1, 1, lat_wgts ) !------------------------------------- !... linear interpolation !------------------------------------- do m = 1,tropo_month_n call lininterp( tropo_p_in(:,:,m), tropo_lon_n, tropo_lat_n, & tmp_out, 1 , lon_wgts, lat_wgts) tropo_bc(id)%tropo_p_loc(i,j,m) = tmp_out(1) end do end do end do call lininterp_finish( lon_wgts ) call lininterp_finish( lat_wgts ) deallocate(tropo_lon) deallocate(tropo_lat) deallocate(tropo_p_in) call wrf_message( ' ' ) write(message,*) 'tropopause_init: DONE intialized domain ',id call wrf_message( trim(message) ) call wrf_message( ' ' ) endif tropo_bc_allocated end subroutine tropopause_init !=========================================================================== subroutine handle_ncerr( ret, mes ) !--------------------------------------------------------------------- ! ... netcdf error handling routine !--------------------------------------------------------------------- implicit none !--------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------- integer, intent(in) :: ret character(len=*), intent(in) :: mes include 'netcdf.inc' if( ret /= nf_noerr ) then call wrf_message( trim(mes) ) call wrf_message( trim(nf_strerror(ret)) ) call wrf_abort end if end subroutine handle_ncerr !=========================================================================== ! This routine uses an implementation of Reichler et al. [2003] done by ! Reichler and downloaded from his web site. This is similar to the WMO ! routines, but is designed for GCMs with a coarse vertical grid. subroutine tropopause_twmo (id, t_phy, p_phy, p8w, zmid, z8w, & tropo_lev, tropo_p,tropo_z, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) implicit none !---------------------------------------------------- ! input arguments !---------------------------------------------------- integer, intent(in) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real, dimension( ims:ime , kms:kme , jms:jme ), & intent(in) :: t_phy, & ! t at mid-level p_phy, & ! p at mid_level (Pa) zmid, & ! z at mid_level (meters) z8w, & ! z at interface (meters) p8w ! p at interface (Pa) real, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_p, & ! tropopause pressure (Pa) tropo_z ! tropopause height (meters) integer, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_lev ! tropopause level !---------------------------------------------------- ! local scalars !---------------------------------------------------- real, dimension( kts:kte) :: t_temp ! store t from top to bottom real, dimension( kts:kte) :: p_temp ! store p from top to bottom real(r8), parameter :: gamma = -0.002_r8 ! K/m real(r8), parameter :: plimu = 45000._r8 ! Pa real(r8), parameter :: pliml = 7500._r8 ! Pa integer :: i integer :: j integer :: k integer :: kk integer :: pLev real(r8) :: tP ! tropopause pressure (Pa) real(r8) :: dZdlogP !---------------------------------------------------- pLev = kte - kts + 1 ! Iterate over all of the grids. do i = its,iend do j = jts,jend ! subroutine twmo expects variables from top to bottom ! t_phy and p_phy are from bottom to top do k = kts,kte kk = pLev - k + 1 t_temp(kk) = t_phy(i,k,j) p_temp(kk) = p_phy(i,k,j) end do ! Use the routine from Reichler. call twmo(pLev, t_temp, p_temp, plimu, pliml, gamma, tP) ! if successful, store the results and find the level and temperature. if (tP > 0) then ! Find the associated level, from bottom to top do k = kts,kte-1 if (tP >= p8w(i,k,j)) then tropo_lev(i,j) = k - 1 ! needed for wrf, from bottom to top tropo_p (i,j) = tP exit end if end do !---------------------------------------------------------- ! get tropopause height ! Intrepolate the geopotential height linearly against log(P) k = k - 1 ! needed for wrf, from bottom to top ! Is the tropoause at the midpoint? if (tropo_p(i,j) == p_phy(i,k,j)) then tropo_z(i,j)= zmid(i,k,j) else if (tropo_p(i,j) < p_phy(i,k,j)) then ! It is below the midpoint. Make sure we aren't at the bottom. if ( k > kts ) then dZdlogP = (zmid(i,k,j) - z8w(i,k-1,j)) / & (log(p_phy(i,k,j)) - log(p8w(i,k-1,j)) ) tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP end if else ! It is above the midpoint? Make sure we aren't at the top. if ( k < kte ) then dZdlogP = (zmid(i,k,j) - z8w(i,k,j)) / & (log(p_phy(i,k,j)) - log(p8w(i,k,j)) ) tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP end if end if !---------------------------------------------------------- end if ! if ( tropo_lev(i,j) == NOTFOUND ) then ! write (*,*) 'tropopause_twmo: NOTFOUND at id, i, j = ', id, i, j ! end if end do end do end subroutine tropopause_twmo !=========================================================================== ! This routine is an implementation of Reichler et al. [2003] done by ! Reichler and downloaded from his web site. Minimal modifications were ! made to have the routine work within the CAM framework (i.e. using ! CAM constants and types). ! ! NOTE: I am not a big fan of the goto's and multiple returns in this ! code, but for the moment I have left them to preserve as much of the ! original and presumably well tested code as possible. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! determination of tropopause height from gridded temperature data ! ! reference: Reichler, T., M. Dameris, and R. Sausen (2003) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine twmo(level, t, p, plimu, pliml, gamma, trp) integer, intent(in) :: level !jeff ! real(r8), intent(in), dimension(level) :: t, p real, intent(in), dimension(level) :: t, p real(r8), intent(in) :: plimu, pliml, gamma real(r8), intent(out) :: trp real(r8), parameter :: deltaz = 2000.0_r8 real(r8) :: pmk, pm, a, b, tm, dtdp, dtdz real(r8) :: ag, bg, ptph real(r8) :: pm0, pmk0, dtdz0 real(r8) :: p2km, asum, aquer real(r8) :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2 integer :: icount, jj integer :: j trp=-99.0_r8 ! negative means not valid lev_loop : & do j=level,2,-1 ! dt/dz pmk= .5_r8 * (p(j-1)**cnst_kap+p(j)**cnst_kap) pm = pmk**(1/cnst_kap) a = (t(j-1)-t(j))/(p(j-1)**cnst_kap-p(j)**cnst_kap) b = t(j)-(a*p(j)**cnst_kap) tm = a * pmk + b dtdp = a * cnst_kap * (pm**cnst_ka1) dtdz = cnst_faktor*dtdp*pm/tm ! dt/dz valid? ! if (j.eq.level) go to 999 ! no, start level, initialize first ! if (dtdz.le.gamma) go to 999 ! no, dt/dz < -2 K/km ! if (pm.gt.plimu) go to 999 ! no, too low if( j == level .or. dtdz <= gamma .or. pm > plimu ) then pm0 = pm pmk0 = pmk dtdz0 = dtdz cycle lev_loop endif ! dtdz is valid, calculate tropopause pressure if (dtdz0 < gamma) then ag = (dtdz-dtdz0) / (pmk-pmk0) bg = dtdz0 - (ag * pmk0) ptph = exp(log((gamma-bg)/ag)/cnst_kap) else ptph = pm endif ! if (ptph.lt.pliml) go to 999 ! if (ptph.gt.plimu) go to 999 if( ptph < pliml .or. ptph > plimu ) then pm0 = pm pmk0 = pmk dtdz0 = dtdz cycle lev_loop endif ! 2nd test: dtdz above 2 km must not exceed gamma p2km = ptph + deltaz*(pm/tm)*cnst_faktor ! p at ptph + 2km asum = 0.0_r8 ! dtdz above icount = 0 ! number of levels above ! test until apm < p2km lev_loop_a : & do jj=j,2,-1 pmk2 = .5_r8 * (p(jj-1)**cnst_kap+p(jj)**cnst_kap) ! p mean ^kappa pm2 = pmk2**(1/cnst_kap) ! p mean if( pm2 > ptph ) then ! doesn't happen cycle lev_loop_a endif if( pm2 < p2km ) then ! ptropo is valid trp = ptph exit lev_loop endif a2 = (t(jj-1)-t(jj)) ! a a2 = a2/(p(jj-1)**cnst_kap-p(jj)**cnst_kap) b2 = t(jj)-(a2*p(jj)**cnst_kap) ! b tm2 = a2 * pmk2 + b2 ! T mean dtdp2 = a2 * cnst_kap * (pm2**(cnst_kap-1)) ! dt/dp dtdz2 = cnst_faktor*dtdp2*pm2/tm2 asum = asum+dtdz2 icount = icount+1 aquer = asum/float(icount) ! dt/dz mean ! discard ptropo ? if( aquer <= gamma ) then ! dt/dz above < gamma pm0 = pm pmk0 = pmk dtdz0 = dtdz exit lev_loop_a endif enddo lev_loop_a enddo lev_loop end subroutine twmo !=========================================================================== ! Read the tropopause pressure in from a file containging a climatology. The ! data is interpolated to the current dat of year and latitude. ! ! NOTE: The data is read in during tropopause_init and stored in the module ! variable tropo_bc(id)%tropo_p_loc subroutine tropopause_climate (id, current_date_char, & p_phy, p8w, zmid, z8w, & tropo_lev, tropo_p, tropo_z, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) implicit none !---------------------------------------------------- ! input arguments !---------------------------------------------------- integer, intent(in) :: id, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte real, dimension( ims:ime , kms:kme , jms:jme ), & intent(in) :: p_phy, & ! p at mid-level (Pa) zmid, & ! z at mid_level (meters) z8w, & ! z at interface (meters) p8w ! p at interface (Pa) real, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_p, & ! tropopause pressure (Pa) tropo_z ! tropopause height (meters) integer, dimension( ims:ime, jms:jme ), & intent(inout) :: tropo_lev ! tropopause level CHARACTER (LEN=256),intent(in) :: current_date_char !---------------------------------------------------- ! Local Variables real :: del_time integer :: last integer :: next real(r8) :: tP ! tropopause pressure (Pa) real(r8) :: dZdlogP integer :: i integer :: j integer :: k CHARACTER(LEN=132) :: message !----------------------------------------------------------------------- if (any(tropo_lev(its:iend,jts:jend) == NOTFOUND)) then !------------------------------------------------------------------------ ! setup time interpolation !------------------------------------------------------------------------ call get_time_interp_factors( current_date_char, last, next, del_time ) ! Iterate over all of the grids. do i = its,iend do j = jts,jend ! Skip grids in which the tropopause has already been found. if (tropo_lev(i,j) == NOTFOUND) then !-------------------------------------------------------- ! ... get tropopause level from climatology !-------------------------------------------------------- ! Interpolate the tropopause pressure. tP = tropo_bc(id)%tropo_p_loc(i,j,last) & + (tropo_bc(id)%tropo_p_loc(i,j,next) & - tropo_bc(id)%tropo_p_loc(i,j,last))*del_time if (tP > 0) then ! Find the associated level, from bottom to top do k = kts,kte-1 if (tP >= p8w(i,k,j)) then tropo_lev(i,j) = k - 1 ! needed for wrf, from bottom to top tropo_p (i,j) = tP exit end if end do !---------------------------------------------------------- ! get tropopause height ! Intrepolate the geopotential height linearly against log(P) k = k - 1 ! needed for wrf ! Is the tropoause at the midpoint? if (tropo_p(i,j) == p_phy(i,k,j)) then tropo_z(i,j)= zmid(i,k,j) else if (tropo_p(i,j) < p_phy(i,k,j)) then ! It is below the midpoint. Make sure we aren't at the bottom. if ( k > kts ) then dZdlogP = (zmid(i,k,j) - z8w(i,k-1,j)) / & (log(p_phy(i,k,j)) - log(p8w(i,k-1,j)) ) tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP end if else ! It is above the midpoint? Make sure we aren't at the top. if ( k < kte ) then dZdlogP = (zmid(i,k,j) - z8w(i,k,j)) / & (log(p_phy(i,k,j)) - log(p8w(i,k,j)) ) tropo_z(i,j)= zmid(i,k,j) + (log(tropo_p(i,j)) - log(p_phy(i,k,j))) * dZdlogP end if end if !---------------------------------------------------------- end if end if end do end do end if if (any(tropo_lev(its:iend,jts:jend) == NOTFOUND)) then write(message,*) 'tropopause_climate: Warning: some tropopause levels still NOTFOUND' else write(message,*) 'tropopause_climate: Warning: Done finding tropopause' end if call wrf_message( trim(message) ) end subroutine tropopause_climate !=========================================================================== subroutine get_time_interp_factors(current_date_char , last, next, del_time ) !----------------------------------------------------------------------------- ! ... setup the time interpolation !----------------------------------------------------------------------------- use module_date_time, only : get_julgmt implicit none !----------------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------------- CHARACTER (LEN=256),intent(in) :: current_date_char integer, intent(out) :: next, last real, intent(out) :: del_time !----------------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------------- integer, parameter :: day_of_year(12) = (/ 16, 45, 75, 105, 136, 166, 197, & 228, 258, 289, 319, 350 /) integer :: julyr , julday real :: gmt real :: calday integer :: m !--------------------------------------------------------------------------------- call get_julgmt ( current_date_char, julyr, julday, gmt ) calday = real(julday) + gmt if( calday < day_of_year(1) ) then next = 1 last = 12 del_time = (365. + calday - day_of_year(12)) & / (365. + day_of_year(1) - day_of_year(12)) else if( calday >= day_of_year(12) ) then next = 1 last = 12 del_time = (calday - day_of_year(12)) & / (365. + day_of_year(1) - day_of_year(12)) else do m = 11,1,-1 if( calday >= day_of_year(m) ) then exit end if end do last = m next = m + 1 del_time = (calday - day_of_year(m)) / (day_of_year(m+1) - day_of_year(m)) end if del_time = max( min( 1.,del_time ),0. ) end subroutine get_time_interp_factors !=========================================================================== end module module_tropopause