module edt_mod !======================================================================= ! ! ! ! EDT (Entrainment and Diagnostic Turbulence) MODULE ! ! ! February 2002 ! Contact person: Steve Klein ! ! ! These routines calculate the diffusivity coefficients for ! momentum and temperature-moisture-scalars using the moist ! thermodynamcs modules based on: ! ! H. Grenier and C. Bretherton, 2001: A moist PBL parameterization ! for large-scale models and its application to subtropical ! cloud-topped marine boundary layers. Mon. Wea. Rev., 129, ! 357-377. ! ! The actual routine is not described in this paper but is ! a simplified extension of the parameterization discussed ! here. The original code, given to Steve Klein from ! Chris Bretherton in May 2001, was tested in the NCAR ! atmospheric model, formerly known as CCM. The code has ! been adapted for the FMS system by Steve Klein and Paul ! Kushner. ! ! ! To quote the Bretherton and Grenier description: ! ! Driver routine to compute eddy diffusion coefficients for ! momentum, moisture, trace constituents and static energy. Uses ! first order closure for stable turbulent layers. For convective ! layers, an entrainment closure is used, coupled to a diagnosis ! of layer-average TKE from the instantaneous thermodynamic and ! velocity profiles. Convective layers are diagnosed by extending ! layers of moist static instability into adjacent weakly stably ! stratified interfaces, stopping if the stability is too strong. ! This allows a realistic depiction of dry convective boundary ! layers with a downgradient approach." ! ! Authors: Herve Grenier, 06/2000, Chris Bretherton 09/2000 ! !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! ! outside modules used ! use constants_mod, only: grav,vonkarm,cp_air,rdgas,rvgas,hlv,hls, & tfreeze,radian use mpp_mod, only: input_nml_file use fms_mod, only: file_exist, open_namelist_file, error_mesg, FATAL,& NOTE, mpp_pe, mpp_root_pe, close_file, read_data, & write_data, write_version_number, stdlog, & open_file, check_nml_error use fms_io_mod, only: register_restart_field, restart_file_type, & save_restart, restore_state use diag_manager_mod, only: register_diag_field, send_data use time_manager_mod, only: time_type, get_date, month_name use monin_obukhov_mod, only: mo_diff use sat_vapor_pres_mod, only: compute_qs implicit none private !----------------------------------------------------------------------- ! ! public interfaces public edt, edt_init, edt_end, edt_on, qaturb, qcturb,tblyrtau !----------------------------------------------------------------------- ! ! global storage variable ! real, allocatable, dimension(:,:,:) :: qaturb ! cloud fraction diagnosed ! from turbulence model ! (fraction) real, allocatable, dimension(:,:,:) :: qcturb ! cloud condensate ! diagnosed from turb. ! model (kg liq/kg air) real, allocatable, dimension(:,:,:) :: tblyrtau ! turbulent layer ! time scale ! (seconds) real, allocatable, dimension(:,:,:) :: sigmas ! standard deviation of ! water perturbation ! (kg water/kg air) type(restart_file_type), pointer, save :: edt_restart => NULL() !----------------------------------------------------------------------- ! ! set default values to namelist parameters ! character :: sftype = "z" ! method for calculating sat frac real :: qcminfrac = 1.e-3 ! Min condensate counted as cloud ! as a fraction of qsat. logical :: use_qcmin = .true. ! Use qcminfrac as indicator of ! of cloud top. logical :: use_extrapolated_ql = .false. ! should the layer top ! liquid water be used to ! estimate the evaporative ! enhancement integer :: n_print_levels = 14 ! how many of the lowest levels ! should be printed out integer, dimension(2) :: edt_pts = 0 ! the global indices for i,j ! at which diagnostics will ! print out logical :: do_print = .false. ! should selected variables ! be sent to logfile logical :: column_match = .false. ! should this column be printed ! out? integer :: dpu = 0 ! unit # for do_print output real :: min_adj_time = 0.25 ! minimum adjustment time of ! turbulent layer for thermo- ! dynamics, as a fraction of ! the physics time step. ! ! the following quantities only deal with the gaussian cloud model ! logical :: do_gaussian_cloud = .false. real :: kappa = 0.5 ! absolute value of the correlation ! coefficient between sli and qt ! fluctions real :: mesovar = 0.02 ! amplitude to mesoscale fluctuations ! in sigma-s as a fraction of the ! saturation specific humidity integer, parameter :: MAX_PTS = 20 integer, dimension (MAX_PTS) :: i_edtprt_gl=0, j_edtprt_gl=0 real, dimension(MAX_PTS) :: lat_edtprt=999., lon_edtprt=999. integer :: num_pts_ij = 0 integer :: num_pts_latlon = 0 namelist /edt_nml/sftype, qcminfrac, use_qcmin,kappa, mesovar,edt_pts, & do_gaussian_cloud, n_print_levels, use_extrapolated_ql, & min_adj_time, i_edtprt_gl, j_edtprt_gl, & num_pts_ij, num_pts_latlon, lat_edtprt, lon_edtprt integer :: num_pts ! total number of columns in which ! diagnostics are desired !--------------------------------------------------------------------- ! deglon1 and deglat1 are the longitude and latitude of the columns ! at which diagnostics will be calculated (degrees). !--------------------------------------------------------------------- real, dimension(:), allocatable :: deglon1, deglat1 !--------------------------------------------------------------------- ! iradprt and jradprt are the processor-based i and j coordinates ! of the desired diagnostics columns. !--------------------------------------------------------------------- integer, dimension(:), allocatable :: j_edtprt, i_edtprt !--------------------------------------------------------------------- ! do_raddg is an array of logicals indicating which latitude rows ! belonging to the processor contain diagnostics columns. !--------------------------------------------------------------------- logical, dimension(:), allocatable :: do_edt_dg !----------------------------------------------------------------------- ! ! diagnostic fields ! character(len=10) :: mod_name = 'edt' real :: missing_value = -999. integer :: id_fq_cv_int, id_fq_cv_top, id_fq_cv_bot, & id_fq_st, id_n2, id_s2, id_ri, & id_leng, id_bprod, id_sprod, id_trans, & id_diss, id_qaedt, id_qcedt, id_tauinv, & id_eddytau, id_fq_turb, id_sigmas, id_radf, & id_sh, id_sm, id_gh, id_evhc !----------------------------------------------------------------------- ! ! set default values to parameters ! logical :: edt_on = .false. logical :: init = .false. real, parameter :: small = 1.e-8 real, parameter :: fpi = 3.14159 ! pi real, parameter :: d608 = (rvgas-rdgas)/rdgas real, parameter :: d622 = rdgas/rvgas real, parameter :: d378 = 1. - d622 real, parameter :: frac_sfclyr = 0.1 ! height of surface layer top as a ! fraction of the pbl height !----------------------------------------------------------------------- ! ! the following parameters are those defined only in the routines ! provided by Chris Bretherton and Herve Grenier ! real, parameter :: ntzero = 1.e-10 ! not zero, used to set min value ! to s2, shear-squared. real, parameter :: zvir = d608 real, parameter :: b1 = 5.8 ! TKE dissipation = e^3/(b1*leng) real :: b123 = 3.2281 ! b1**(2/3) real, parameter :: tunl = 0.085 ! Asympt leng = ! tunl*(turb lyr depth) real, parameter :: alph1 = 0.5562 ! Galperin stability fn params real, parameter :: alph2 = -4.3640 real, parameter :: alph3 = -34.6764 real, parameter :: alph4 = -6.1272 real, parameter :: alph5 = 0.6986 real, parameter :: ricrit = 0.19 ! Critical Richardson # for turb. real, parameter :: mu = 70. ! used in finding e/ebar in ! convective layers (CLs) real, parameter :: rinc = -0.5 ! Min W/ for incorp into CL ! params governing entr. effic. A=a1l*evhc, evhc=1+a2l*a3l*L*ql/jt2slv ! where ql is cloud-top liq. water and jt2slv is the jump in slv across ! the cloud-top entrainment zone. real, parameter :: a1l = 0.10 ! Dry entrainment efficiency param ! Herve set to 0.05 due to excess ! TKE but a1l = 0.10 = 0.2*tunl* ! erat^-1.5 should be the "real" ! value, where erat = /wstar^2 ! for dry CBL = 0.3 real, parameter :: a2l = 15. ! Moist entrainment enhancement ! param Herve's SCCM value was 15 real, parameter :: a3l = 0.8 ! real, parameter :: jbumin = 0.001 ! Min buoyancy jump at an entrain- ! ment interface (m/s2) ! (~ jump in K/30) real, parameter :: evhcmax = 10. ! Max entrainment efficiency real, parameter :: rimaxentr = 0. ! Limiting Ri for entraining turb ! layer ! parameters affecting TKE real, parameter :: rmin = 0.1 ! Min allowable e/ in a CL real, parameter :: rmax = 2.0 ! Max allowable e/ in a CL real, parameter :: tkemax = 20. ! tke capped at tkemax (m2/s2) real, parameter :: tkemin = 1.e-6 ! tke minimum (m2/s2) !----------------------------------------------------------------------- ! ! declare version number ! character(len=128) :: Version = '$Id: edt.F90,v 19.0 2012/01/06 20:09:18 fms Exp $' character(len=128) :: Tagname = '$Name: tikal $' logical :: module_is_initialized = .false. !----------------------------------------------------------------------- ! ! fms module subroutines include: ! ! edt main driver program of the module ! ! edt_init initialization routine ! ! edt_tend adds in the longwave heating rate to the ! global storage variable ! ! edt_end ending routine ! ! ! Grenier-Bretherton subroutines are described after the ! subroutines listed above. ! integer, dimension(1) :: restart_versions = (/ 1 /) contains !======================================================================= ! ! subroutine edt_init ! ! ! this subroutine reads the namelist file and restart data ! and initializes some constants. ! subroutine edt_init(lonb, latb, axes,time,idim,jdim,kdim) !----------------------------------------------------------------------- ! ! variables ! ! ----- ! input ! ----- ! ! idim,jdim,kdim size of the first 3 dimensions ! axes, time variables needed for netcdf diagnostics ! latb, lonb latitudes and longitudes at grid box corners ! ! ! -------- ! internal ! -------- ! ! unit unit number for namelist and restart file ! io internal variable for reading of namelist file ! full indices for full level axes coordinates ! half indices for half level axes coordinates ! !----------------------------------------------------------------------- integer, intent(in) :: idim,jdim,kdim,axes(4) type(time_type), intent(in) :: time real, dimension(:,:), intent(in) :: lonb, latb integer :: unit, io, logunit, ierr, vers, vers2, nn, i, j, id_restart real :: dellat, dellon character(len=4) :: chvers integer, dimension(3) :: full = (/1,2,3/), half = (/1,2,4/) !----------------------------------------------------------------------- ! ! namelist functions #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=edt_nml, iostat=io) ierr = check_nml_error(io,'edt_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file () ierr=1 ; Do While (ierr .ne. 0) Read (unit, nml=edt_nml, iostat=io, End=10) ierr = check_nml_error(io,'edt_nml') enddo 10 Call Close_File (unit) endif #endif !------- write version number and namelist --------- if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) logunit = stdlog() write (logunit, nml=edt_nml) endif !--------------------------------------------------------------------- ! allocate and initialize a flag array which indicates the latitudes ! containing columns where radiation diagnostics are desired. !--------------------------------------------------------------------- allocate (do_edt_dg (size(latb,2)-1) ) do_edt_dg(:) = .false. !------------------------------------------------------------------- ! define the total number of points at which diagnostics are desired. ! points may be specified either by lat-lon pairs or by global index ! pairs. !------------------------------------------------------------------- num_pts = num_pts_latlon + num_pts_ij !------------------------------------------------------------------- ! continue on only if diagnostics are desired in at least one column. !------------------------------------------------------------------- if (num_pts > 0) then !------------------------------------------------------------------- ! if more points are desired than space has been reserved for, print ! a message. !------------------------------------------------------------------- if (num_pts > MAX_PTS) then call error_mesg ( 'edt_mod', & 'must reset MAX_PTS or reduce number of diagnostics points', & FATAL) endif !------------------------------------------------------------------- ! allocate space for arrays which will contain the lat and lon and ! processor-local i and j indices. !------------------------------------------------------------------- allocate ( deglon1 (num_pts)) allocate ( deglat1 (num_pts)) allocate ( j_edtprt (num_pts)) allocate ( i_edtprt (num_pts)) !--------------------------------------------------------------------- ! if any points for diagnostics are specified by (i,j) global ! indices, determine their lat-lon coordinates. assumption is made ! that the deltas of latitude and longitude are uniform over ! the globe. !--------------------------------------------------------------------- do nn=1,num_pts_ij dellat = latb(1,2) - latb(1,1) dellon = lonb(2,1) - lonb(1,1) lat_edtprt(nn + num_pts_latlon) = & (-0.5*acos(-1.0) + (j_edtprt_gl(nn) - 0.5)* & dellat) * radian lon_edtprt(nn + num_pts_latlon) = & (i_edtprt_gl(nn) - 0.5)*dellon*radian end do !-------------------------------------------------------------------- ! determine if the lat/lon values are within the global grid, ! latitude between -90 and 90 degrees and longitude between 0 and ! 360 degrees. !-------------------------------------------------------------------- do nn=1,num_pts j_edtprt(nn) = 0 i_edtprt(nn) = 0 deglat1(nn) = 0.0 deglon1(nn) = 0.0 if (lat_edtprt(nn) .ge. -90. .and. & lat_edtprt(nn) .le. 90.) then else call error_mesg ('edt_mod', & ' invalid latitude for edt diagnostics ', FATAL) endif if (lon_edtprt(nn) .ge. 0. .and. & lon_edtprt(nn) .le. 360.) then else call error_mesg ('edt_mod', & ' invalid longitude for edt diagnostics ', FATAL) endif !-------------------------------------------------------------------- ! determine if the diagnostics column is within the current ! processor's domain. if so, set a logical flag indicating the ! presence of a diagnostic column on the particular row, define the ! i and j processor-coordinates and the latitude and longitude of ! the diagnostics column. !-------------------------------------------------------------------- do j=1,size(latb,2) - 1 if (lat_edtprt(nn) .ge. latb(1,j)*radian .and. & lat_edtprt(nn) .lt. latb(1,j+1)*radian) then do i=1,size(lonb,1) - 1 if (lon_edtprt(nn) .ge. lonb(i,1)*radian & .and.& lon_edtprt(nn) .lt. lonb(i+1,1)*radian) & then do_edt_dg(j) = .true. j_edtprt(nn) = j i_edtprt(nn) = i deglon1(nn) = 0.5*(lonb(i,1) + lonb(i+1,1))*radian deglat1(nn) = 0.5*(latb(1,j) + latb(1,j+1))*radian exit endif end do exit endif end do end do !---------------------------------------------------------------------- ! open a unit for the radiation diagnostics output. !--------------------------------------------------------------------- dpu = open_file ('edt.out', action='write', & threading='multi', form='formatted') do_print = .true. if ( mpp_pe() == mpp_root_pe() ) then call write_version_number(version, tagname) write (dpu ,nml=edt_nml) endif endif ! (num_pts > 0) !----------------------------------------------------------------------- ! ! initialize edt_on edt_on = .TRUE. module_is_initialized = .true. !----------------------------------------------------------------------- ! ! initialize b123 = b1**(2/3) b123 = b1**(2./3.) !----------------------------------------------------------------------- ! ! handle global storage if (allocated(qaturb)) deallocate (qaturb) allocate(qaturb(idim,jdim,kdim)) if (allocated(qcturb)) deallocate (qcturb) allocate(qcturb(idim,jdim,kdim)) if (allocated(tblyrtau)) deallocate (tblyrtau) allocate(tblyrtau(idim,jdim,kdim)) if (allocated(sigmas)) deallocate (sigmas) allocate(sigmas(idim,jdim,kdim+1)) allocate(edt_restart) id_restart = register_restart_field(edt_restart, 'edt.res.nc', 'qaturb' , qaturb ) id_restart = register_restart_field(edt_restart, 'edt.res.nc', 'qcturb' , qcturb ) id_restart = register_restart_field(edt_restart, 'edt.res.nc', 'tblyrtau', tblyrtau ) id_restart = register_restart_field(edt_restart, 'edt.res.nc', 'sigmas' , sigmas ) if (File_Exist('INPUT/edt.res.nc')) then if (mpp_pe() == mpp_root_pe() ) then call error_mesg ('edt_mod', 'Reading netCDF formatted restart file: INPUT/edt.res.nc', & NOTE) endif call restore_state(edt_restart) elseif (File_Exist('INPUT/edt.res')) then call error_mesg ('edt_mod', 'Native format restart file read no longer supported.',& FATAL) ! unit = Open_restart_File (FILE='INPUT/edt.res', ACTION='read') ! read (unit, iostat=io, err=142) vers, vers2 !142 continue ! !!-------------------------------------------------------------------- !! if eor is not encountered, then the file includes tdtlw as the !! first record (which this read statement read). that data is not !! needed; note this and continue by reading next record. !!-------------------------------------------------------------------- ! if (io == 0) then ! call error_mesg ('edt_mod', & ! 'reading pre-version number edt.res file, '//& ! 'ignoring tdtlw', NOTE) ! !!-------------------------------------------------------------------- !! if the first record was only one word long, then the file is a !! newer one, and that record was the version number, read into vers. !! if it is not a valid version, stop execution with a message. !!-------------------------------------------------------------------- ! else ! if (.not. any(vers == restart_versions) ) then ! write (chvers, '(i4)') vers ! call error_mesg ('edt_mod', & ! 'restart version ' // chvers//' cannot be read '//& ! 'by this version of edt_mod.', FATAL) ! endif ! endif ! !!--------------------------------------------------------------------- ! call read_data (unit, qaturb) ! call read_data (unit, qcturb) ! call read_data (unit, tblyrtau) ! call read_data (unit, sigmas) ! call Close_File (unit) else qaturb (:,:,:) = 0. qcturb (:,:,:) = 0. tblyrtau(:,:,:) = 0. sigmas (:,:,:) = 0. endif !----------------------------------------------------------------------- ! ! register diagnostic fields id_fq_cv_int = register_diag_field (mod_name, 'fq_cv_int', axes(half), & time, 'Frequency that the interface is in '// & 'the interior of a convective layer from EDT', & 'none', missing_value=missing_value ) id_fq_cv_top = register_diag_field (mod_name, 'fq_cv_top', axes(half), & time, 'Frequency that the interface is at '// & 'the top of a convective layer from EDT', & 'none', missing_value=missing_value ) id_fq_cv_bot = register_diag_field (mod_name, 'fq_cv_bot', axes(half), & time, 'Frequency that the interface is at '// & 'the bottom of a convective layer from EDT', & 'none', missing_value=missing_value ) id_fq_st = register_diag_field (mod_name, 'fq_st', axes(half), & time, 'Frequency of stable turbulence from EDT', & 'none', missing_value=missing_value ) id_fq_turb = register_diag_field (mod_name, 'fq_turb', axes(full), & time, 'Frequency that the layer is fully '// & 'turbulent from EDT','none', missing_value=missing_value ) id_n2 = register_diag_field (mod_name, 'n2', axes(half), & time,'Moist Vaisala Frequency from EDT', & '(1/sec)**2', missing_value=missing_value ) id_s2 = register_diag_field (mod_name, 's2', axes(half), & time,'Shear vector magnitude squared from EDT', & '(1/sec)**2', missing_value=missing_value ) id_ri = register_diag_field (mod_name, 'ri', axes(half), & time, 'Moist Richardson number from EDT', & 'none', missing_value=missing_value ) id_leng = register_diag_field (mod_name, 'leng', axes(half), & time, 'Turbulent Length Scale from EDT', & 'meters', missing_value=missing_value ) id_bprod = register_diag_field (mod_name, 'bprod', axes(half), & time, 'Buoyancy production of TKE from EDT', & '(meters**2)/(sec**3)', missing_value=missing_value ) id_sprod = register_diag_field (mod_name, 'sprod', axes(half), & time, 'Shear production of TKE from EDT', & '(meters**2)/(sec**3)', missing_value=missing_value ) id_trans = register_diag_field (mod_name, 'trans', axes(half), & time, 'TKE transport from EDT', & '(meters**2)/(sec**3)', missing_value=missing_value ) id_diss = register_diag_field (mod_name, 'diss', axes(half), & time, 'TKE dissipation from EDT', & '(meters**2)/(sec**3)', missing_value=missing_value ) id_radf = register_diag_field (mod_name, 'radf', axes(half), & time, 'TKE radiative forcing from EDT', & '(meters**2)/(sec**3)', missing_value=missing_value ) id_sh = register_diag_field (mod_name, 'sh', axes(half), & time, 'Galperin heat stability coefficient from EDT', & 'none', missing_value=missing_value ) id_sm = register_diag_field (mod_name, 'sm', axes(half), & time, 'Galperin momentum stability coefficient from EDT', & 'none', missing_value=missing_value ) id_gh = register_diag_field (mod_name, 'gh', axes(half), & time, 'Galperin stability ratio from EDT', & 'none', missing_value=missing_value ) id_evhc = register_diag_field (mod_name, 'evhc', axes(half), & time, 'Evaporative cooling enhancement factor from EDT', & 'none', missing_value=missing_value ) id_sigmas = register_diag_field (mod_name, 'sigmas', axes(half), & time, 'Std. dev. of water perturbation from EDT', & '(kg water)/(kg air)', missing_value=missing_value ) id_qaedt = register_diag_field (mod_name, 'qaedt', axes(full), & time, 'statistical cloud fraction from EDT', & 'fraction', missing_value=missing_value ) id_qcedt = register_diag_field (mod_name, 'qcedt', axes(full), & time, 'statistical cloud condensate from EDT', & 'kg condensate/kg air', missing_value=missing_value ) id_tauinv = register_diag_field (mod_name, 'tauinv', axes(half),& time, 'inverse large-eddy turnover time from EDT', & '1/second', missing_value=missing_value ) id_eddytau = register_diag_field (mod_name, 'eddytau', axes(full), & time, 'large-eddy turnover time from EDT', & 'seconds', missing_value=missing_value ) !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine edt_init ! !======================================================================= !======================================================================= ! ! subroutine edt ! ! ! this subroutine is the main driver program to the routines ! provided by Chris Bretherton ! subroutine edt(is,ie,js,je,dt,time,tdtlw_in, u_star,b_star,q_star,t,qv,ql,qi,qa, & u,v,z_full,p_full,z_half,p_half,stbltop,k_m,k_t,pblh,kbot,tke) !----------------------------------------------------------------------- ! ! variables ! ! ----- ! input ! ----- ! ! is,ie,js,je i,j indices marking the slab of model working on ! dt physics time step (seconds) ! time variable needed for netcdf diagnostics ! u_star friction velocity (m/s) ! b_star buoyancy scale (m/(s**2)) ! q_star moisture scale (kg vapor/kg air) ! ! three dimensional fields on model full levels, reals dimensioned ! (:,:,pressure), third index running from top of atmosphere to ! bottom ! ! t temperature (K) ! qv water vapor specific humidity (kg vapor/kg air) ! ql liquid water specific humidity (kg cond/kg air) ! qi ice water specific humidity (kg cond/kg air) ! qa cloud fraction ! u zonal wind (m/s) ! v meridional wind (m/s) ! z_full height of full levels (m) ! p_full pressure (Pa) ! ! the following two fields are on the model half levels, with ! size(z_half,3) = size(t,3) +1, z_half(:,:,size(z_half,3)) ! must be height of surface (if you are not using eta-model) ! ! z_half height at half levels (m) ! p_half pressure at half levels (Pa) ! ! ------ ! output ! ------ ! ! stbltop maximum altitude the very stable boundary layer ! is permitted to operate ! ! The following variables are defined at half levels and are ! dimensions 1:nlev+1. ! ! k_m diffusivity for momentum (m**2/s) ! k_t diffusivity for temperature and scalars (m**2/s) ! ! k_m and k_t are defined at half-levels so that size(k_m,3) ! should be at least as large as size(t,3). Note, however, that ! only the returned values at levels 2 to size(t,3) are ! meaningful; other values will be returned as zero. ! ! -------------- ! optional input ! -------------- ! ! kbot integer indicating the lowest true layer of atmosphere ! ! --------------- ! optional output ! --------------- ! ! pblh depth of planetary boundary layer (m) ! tke turbulent kinetic energy (m*m)/(s*s) ! ! -------- ! internal ! -------- ! ! z_surf height of surface (m) ! z_full_ag height of full model levels above the surface (m) ! qt total water specific humidity (kg water/kg air) ! qc cloud condensate spec. hum. (kg condensate/kg air) ! qsl saturation specific humidity at the midpoint pressure ! and liquid-ice water temperature (kg water/kg air) ! dqsldtl temperature derivative of qsl (kg H20/kg air/Kelvin) ! hleff effective latent heat of vaporization (J/kg) ! sli ice-liq water static energy (J/kg) ! sliv ice-liq water virtual static energy (J/kg) ! khfs surface kinematic heat flux (K*(m/s)) ! kqfs surface kinematic vapor flux (kg/kg)*(m/s) ! slislope sli slope wrt pressure in thermo layer (J/kg/Pa) ! qtslope qt slope wrt pressure in thermo layer (kg/kg/Pa) ! qxtop saturation excess at the top of the layer ! (kg wat/kg air) ! qxbot saturation excess at the bottom of the layer ! (kg wat/kg air) ! sfuh saturated fraction in upper half-layer ! sflh sflh saturated fraction in lower half-layer ! sfclyr_h height of the surface layer top (m) ! ql_new new value of cloud liquid (kg cond/kg air) ! qi_new new value of cloud ice (kg cond/kg air) ! qa_new new value of cloud fraction ! mask real array indicating the point is above the surface & ! if equal to 1.0 and indicating the point is below the & ! surface if equal to 0. ! ! mineddytau minimum value to adjustment time (1/sec) ! ! ! the following variables are defined on model half levels ! (1:kdim+1) ! ! z_half_ag height of half model levels above the surface (m) ! chu heat var. coef for dry states (1/m) ! chs heat var. coef for sat states (1/m) ! cmu moisture var. coef for dry states (kg/kg)*(m/s*s) ! cms moisture var. coef for sat states (kg/kg)*(m/s*s) ! n2 moist squared buoyancy freq (1/s*s) ! s2 squared deformation, or shear vector mag. (1/s*s) ! ri gradient Richardson number ! ! formulas for selected internal variables ! ---------------------------------------- ! ! qt = qv + ql + qi qc = ql + qi ! sli = cp*T + g*z - hleff*qc sliv = sli*(1.+d608*qt) ! khfs = mean of (w'T') kqfs = mean of (w'q') ! !----------------------------------------------------------------------- integer, intent(in) :: is,ie,js,je real, intent(in) :: dt type(time_type), intent(in) :: time real, intent(in), dimension(:,:) :: u_star, b_star, q_star real, intent(in), dimension(:,:,:) :: tdtlw_in, & t,qv,ql,qi,qa, & u, v, & z_full, p_full, & z_half, p_half real, intent(out), dimension(:,:) :: stbltop, pblh real, intent(out), dimension(:,:,:) :: k_m,k_t integer, intent(in), dimension(:,:), optional:: kbot real, intent(out), dimension(:,:,:),optional:: tke integer :: i,j,k,kk,ibot integer :: ipt,jpt integer :: nlev,nlat,nlon integer, dimension(4,size(t,1),size(t,2),size(t,3)+1):: turbtype real :: khfs,kqfs real :: sfclyr_h_max real :: mineddytau real, dimension(size(t,1),size(t,2)) :: z_surf, sfclyr_h real, dimension(size(t,1),size(t,2),size(t,3)) :: z_full_ag real, dimension(size(t,1),size(t,2),size(t,3)) :: qt,sli,sliv real, dimension(size(t,1),size(t,2),size(t,3)) :: esl, qsl, dqsldtl real, dimension(size(t,1),size(t,2),size(t,3)) :: hleff,density real, dimension(size(t,1),size(t,2),size(t,3)) :: slislope, qtslope real, dimension(size(t,1),size(t,2),size(t,3)) :: qxtop, qxbot real, dimension(size(t,1),size(t,2),size(t,3)) :: qa_new, qc_new real, dimension(size(t,1),size(t,2),size(t,3)) :: mask, isturb real, dimension(size(t,1),size(t,2),size(t,3)) :: eddytau,tauinvtmp real, dimension(size(t,1),size(t,2),size(t,3)+1):: z_half_ag real, dimension(size(t,1),size(t,2),size(t,3)+1):: n2, s2, ri, tauinv real, dimension(size(t,1),size(t,2),size(t,3)+1):: leng, bprod, sprod real, dimension(size(t,1),size(t,2),size(t,3)+1):: trans, diss, evhc real, dimension(size(t,1),size(t,2),size(t,3)+1):: radf, sh, sm, gh real, dimension(size(t,1),size(t,2),size(t,3)+1):: mask3, tmpdat real, dimension(size(t,1),size(t,2),size(t,3)+1):: k_t_mo, k_m_mo integer, dimension(4,size(t,3)+1) :: gb_turbtype real :: gb_pblh, gb_u_star real, dimension(size(t,3)) :: gb_t, gb_u, gb_v real, dimension(size(t,3)) :: gb_qv real, dimension(size(t,3)) :: gb_dqsldtl, gb_qa_new, gb_qc_new real, dimension(size(t,3)) :: gb_sli, gb_qt real, dimension(size(t,3)) :: gb_qc, gb_sliv, gb_sflh real, dimension(size(t,3)) :: gb_slisl, gb_qtsl, gb_sfuh real, dimension(size(t,3)) :: gb_tdtlw, gb_hleff real, dimension(size(t,3)) :: gb_qsl, gb_esl, gb_qxtop real, dimension(size(t,3)) :: gb_qxbot, gb_density, gb_z_full real, dimension(size(t,3)) :: gb_p_full, gb_isturb, gb_qltop real, dimension(size(t,3)+1):: gb_cmu, gb_chu real, dimension(size(t,3)+1):: gb_chs, gb_cms real, dimension(size(t,3)+1):: gb_k_m, gb_k_t, gb_tke real, dimension(size(t,3)+1):: gb_n2, gb_s2, gb_ri real, dimension(size(t,3)+1):: gb_leng, gb_bprod, gb_sprod real, dimension(size(t,3)+1):: gb_diss, gb_trans, gb_tauinv real, dimension(size(t,3)+1):: gb_z_half, gb_sigmas, gb_p_half real, dimension(size(t,3)+1):: gb_radf, gb_sh, gb_sm real, dimension(size(t,3)+1):: gb_gh, gb_evhc logical :: used, topfound integer, dimension(MAX_PTS) :: nsave integer :: iloc(MAX_PTS), jloc(MAX_PTS), nn, npts, nnsave integer :: year, month, day, hour, minute, second character(len=16) :: mon !----------------------------------------------------------------------- ! ! initialize variables pblh = 0.0 k_t = 0.0 k_m = 0.0 tke = 0.0 turbtype = 0 isturb = 0.0 n2 = 0.0 s2 = 0.0 ri = 0.0 tauinv = 0.0 leng = 0.0 bprod = 0.0 sprod = 0.0 trans = 0.0 diss = 0.0 radf = 0.0 sh = 0.0 sm = 0.0 evhc = 0.0 gh = 0.0 slislope = 0.0 qtslope = 0.0 qxbot = 0.0 qxtop = 0.0 qc_new = 0.0 qa_new = 0.0 mineddytau = min_adj_time * dt !----------------------------------------------------------------------- ! ! compute height above surface ! nlev = size(t,3) nlat = size(t,2) nlon = size(t,1) mask = 1.0 if (present(kbot)) then do j=1,nlat do i=1,nlon z_surf(i,j) = z_half(i,j,kbot(i,j)+1) enddo enddo else z_surf(:,:) = z_half(:,:,nlev+1) end if do k = 1, nlev z_full_ag(:,:,k) = z_full(:,:,k) - z_surf(:,:) z_half_ag(:,:,k) = z_half(:,:,k) - z_surf(:,:) end do z_half_ag(:,:,nlev+1) = z_half(:,:,nlev+1) - z_surf(:,:) if (present(kbot)) then where (z_full_ag < 0) mask = 0.0 endwhere end if !----------------------------------------------------------------------- ! ! Calculate saturation specific humidity and its temperature ! derivative, the effective latent heat ! ! These are calculated according to the formulas: ! ! qsl = d622*esl/ [p_full - (1.-d622)*esl] ! ! dqsldtl = d622*p_full*(desat/dT)/[p_full-(1.-d622)*esl]**2. ! ! ! where d622 = rdgas/rvgas; esl = saturation vapor pressure; ! and desat/dT is the temperature derivative of esl. Note that: ! ! { hlv for t > tfreeze } ! hleff = { 0.05*(t-tfreeze+20.)*hlv + 0.05*(tfreeze-t)*hls } ! { for tfreeze-20.< t < tfreeze} ! { hls for t < tfreeze-20. } ! ! This linear form is chosen because at Tfreeze-20. es = esi, and ! at Tfreeze, es = esl, with linear interpolation in between. ! ! !calculate effective latent heat hleff = (min(1.,max(0.,0.05*(t -tfreeze+20.)))*hlv + & min(1.,max(0.,0.05*(tfreeze -t )))*hls) !calculate qsl and dqsldtl. return es for diagnostic use. call compute_qs (t-(hleff*(ql+qi)/cp_air), p_full, qsl, & esat = esl, dqsdT=dqsldtl) !----------------------------------------------------------------------- ! ! set up specific humidities and static energies ! compute airdensity qt = qv + ql + qi sli = cp_air*t + grav*z_full_ag - hleff*(ql + qi) sliv = sli*(1+zvir*qt) density = p_full/rdgas/(t *(1.+d608*qv-ql-qi)) !----------------------------------------------------------------------- ! ! big loop over points ! ibot = nlev do j=1,nlat npts = 0 if (do_edt_dg(j+js-1) ) then do nn=1,num_pts if ( js == j_edtprt(nn) .and. & i_edtprt(nn) >= is .and. i_edtprt(nn) <= ie) then iloc(npts+1) = i_edtprt(nn) - is + 1 jloc(npts+1) = j_edtprt(nn) - js + 1 nsave(npts+1) = nn npts = npts + 1 endif end do ! (num_points) else ipt = 0 jpt = 0 column_match = .false. endif do i=1,nlon if (npts > 0) then do nn=1,npts ipt = iloc(nn) jpt = jloc(nn) if (i == ipt ) then column_match = .true. nnsave = nsave(nn) exit else column_match = .false. endif end do nn = nnsave else column_match = .false. nn = 0 endif !----------------------------------------------------------- ! should diagnostics be printed out for this column ! !----------------------------------------------------------- ! ! extract column data to pass to caleddy ! if (present(kbot)) ibot = kbot(i,j) gb_t (1:ibot ) = t (i,j,1:ibot ) gb_qv (1:ibot ) = max(0., qv(i,j,1:ibot ) ) gb_qc (1:ibot ) = max(0., ql(i,j,1:ibot ) + & qi(i,j,1:ibot ) ) gb_qt (1:ibot ) = max(0., qt(i,j,1:ibot ) ) gb_sli (1:ibot ) = sli (i,j,1:ibot ) gb_sliv (1:ibot ) = sliv (i,j,1:ibot ) gb_u (1:ibot ) = u (i,j,1:ibot ) gb_v (1:ibot ) = v (i,j,1:ibot ) gb_z_half (1:ibot+1) = z_half_ag (i,j,1:ibot+1) gb_z_full (1:ibot ) = z_full_ag (i,j,1:ibot ) gb_p_half (1:ibot+1) = p_half (i,j,1:ibot+1) gb_p_full (1:ibot ) = p_full (i,j,1:ibot ) gb_qsl (1:ibot ) = qsl (i,j,1:ibot ) gb_esl (1:ibot ) = esl (i,j,1:ibot ) gb_hleff (1:ibot ) = hleff (i,j,1:ibot ) gb_density(1:ibot ) = density (i,j,1:ibot ) gb_dqsldtl(1:ibot ) = dqsldtl (i,j,1:ibot ) gb_tdtlw (1:ibot ) = tdtlw_in (i,j,1:ibot) gb_sigmas (1:ibot+1) = sigmas(is-1+i,js-1+j,1:ibot+1) kqfs = u_star(i,j)*q_star(i,j) khfs = u_star(i,j)*b_star(i,j)*gb_t(ibot)/ & grav gb_u_star = u_star(i,j) !compute sigmas at surface gb_sigmas(ibot+1) = ((mesovar*gb_qsl(ibot))**2.0) + & ((q_star(i,j)-gb_dqsldtl(ibot)*b_star(i,j)*gb_t(ibot)/ & grav)**2.0) gb_sigmas(ibot+1) = sqrt( gb_sigmas(ibot+1) ) / ( 1. + & gb_hleff(ibot)*gb_dqsldtl(ibot)/cp_air ) if (column_match) then call get_date(Time, year, month, day, hour, minute, second) mon = month_name (month) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') '===================================='//& '==================' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ENTERING EDT ' write (dpu,'(a)') ' ' write (dpu,'(a, i6,a,i4,i4,i4,i4)') ' time stamp:', & year, trim(mon), day, & hour, minute, second write (dpu,'(a)') ' DIAGNOSTIC POINT COORDINATES :' write (dpu,'(a)') ' ' write (dpu,'(a,f8.3,a,f8.3)') ' longitude = ', deglon1(nn),& ' latitude = ', deglat1(nn) write (dpu,'(a,i6,a,i6)') ' global i =', i_edtprt_gl(nn), & ' global j = ', j_edtprt_gl(nn) write (dpu,'(a,i6,a,i6)') ' processor i =', i_edtprt(nn), & ' processor j = ',j_edtprt(nn) write (dpu,'(a,i6,a,i6)') ' window i =', ipt, & ' window j = ',jpt write (dpu,'(a)') ' ' write (dpu,'(a)') ' sigmas at the surface .... ' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ',1000.* & gb_sigmas(ibot+1), ' g/kg' write (dpu,'(a,f14.7) ') ' acoef = ', 1./( 1.+ & gb_hleff(ibot)* gb_dqsldtl(ibot)/cp_air ) write (dpu,'(a,f14.7,a)') ' sigmas/a = ', 1000.* & gb_sigmas(ibot+1)*( 1. + gb_hleff(ibot)* & gb_dqsldtl(ibot)/cp_air ), ' g/kg' write (dpu,'(a,f14.7)' ) ' mesovar = ', mesovar write (dpu,'(a,f14.7,a)') ' mesovar*qsl = ',mesovar* & gb_qsl(ibot)*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' turb.fluct = ',1000.* & (q_star(i,j)-gb_dqsldtl(ibot)*b_star(i,j)*gb_t(ibot)/ & grav),' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k T u v '//& ' qv qt ' write (dpu,'(a)') ' (K) (m/s) (m/s) '//& '(g/kg) (g/kg)' write (dpu,'(a)') '------------------------------------'//& '-----------------' write (dpu,'(a)') ' ' do kk = nlev-n_print_levels,nlev write(dpu,18) kk,gb_t(kk),gb_u(kk),gb_v(kk),1000.* & gb_qv(kk), 1000.*gb_qt(kk) end do 18 format(1X,i2,1X,5(f9.4,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' k qa qc sli/cp_air '//& 'sliv/cp_air tdtlw' write (dpu,'(a)') ' (g/kg) (K) '//& ' (K) (K/day)' write (dpu,'(a)') '------------------------------------'//& '-----------------' write (dpu,'(a)') ' ' do kk = nlev-n_print_levels,nlev write(dpu,19) kk,qa(i,j,kk),1000.*gb_qc(kk), & gb_sli(kk)/cp_air,gb_sliv(kk)/cp_air,gb_tdtlw(kk)*86400. enddo 19 format(1X,i2,1X,5(f9.4,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' k z_full z_half p_full p'//& '_half sigmas' write (dpu,'(a)') ' (m) (m) (mb) '//& '(mb) (g/kg)' write (dpu,'(a)') '------------------------------------'//& '-----------------' write (dpu,'(a)') ' ' do kk = nlev-n_print_levels,nlev write(dpu,19) kk,gb_z_full(kk),gb_z_half(kk+1), & gb_p_full(kk)/100.,gb_p_half(kk+1)/100.,1000.* & gb_sigmas(kk+1) enddo write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k esl qsl dqsldtl '//& ' hleff ' write (dpu,'(a)') ' (mb) (g/kg) (g/kg/K) '//& '(MJ/kg)' write (dpu,'(a)') '------------------------------------'//& '-------' write (dpu,'(a)') ' ' do kk = nlev-n_print_levels,nlev write(dpu,19) kk, gb_esl(kk)/100.,gb_qsl(kk)*1000., & gb_dqsldtl(kk)*1000.,gb_hleff(kk)/1.0e+06 enddo write (dpu,'(a)') ' ' end if ! call the initialization routine for variables needed ! by caleddy call trbintd(gb_t, gb_qv, gb_qt, gb_qc, & gb_sli, gb_sliv, gb_u, gb_v, & gb_z_full, gb_z_half, gb_p_full, gb_p_half, & gb_sigmas, gb_qsl, gb_esl, gb_hleff, & gb_dqsldtl, gb_slisl, gb_qtsl, gb_qxtop, & gb_qxbot, gb_qltop, gb_sfuh, gb_sflh, & gb_qc_new, gb_qa_new, gb_chu, gb_chs, & gb_cmu, gb_cms, gb_n2, gb_s2, & gb_ri ) !----------------------------------------------------------- ! ! call caleddy ! call caleddy(gb_u_star, kqfs, khfs, gb_sli, & gb_qt, gb_qc_new, gb_qa_new, gb_sliv, & gb_u, gb_v, gb_p_full, gb_z_full, & gb_sfuh, gb_sflh, gb_slisl, gb_qtsl, & gb_tdtlw, gb_hleff, gb_density,gb_qsl, & gb_dqsldtl,gb_qltop, gb_p_half, gb_z_half, & gb_chu, gb_chs, gb_cmu, gb_cms, & gb_n2, gb_s2, gb_ri, gb_pblh, & gb_turbtype,gb_k_t, gb_k_m, gb_tke, & gb_leng, gb_bprod, gb_sprod, gb_trans, & gb_diss, gb_isturb, gb_tauinv, gb_sigmas, & gb_radf, gb_sh, gb_sm, gb_gh, & gb_evhc) !----------------------------------------------------------- ! ! paste back outputs ! pblh(i,j) = gb_pblh k_t(i,j,1:ibot) = gb_k_t(1:ibot) k_m(i,j,1:ibot) = gb_k_m(1:ibot) tke(i,j,1:ibot+1) = gb_tke(1:ibot+1) turbtype(:,i,j,1:ibot+1) = gb_turbtype(:,1:ibot+1) isturb(i,j,1:ibot) = gb_isturb(1:ibot) n2(i,j,1:ibot+1) = gb_n2(1:ibot+1) s2(i,j,1:ibot+1) = gb_s2(1:ibot+1) ri(i,j,1:ibot+1) = gb_ri(1:ibot+1) tauinv(i,j,1:ibot+1) = gb_tauinv(1:ibot+1) leng(i,j,1:ibot+1) = gb_leng(1:ibot+1) bprod(i,j,1:ibot+1) = gb_bprod(1:ibot+1) sprod(i,j,1:ibot+1) = gb_sprod(1:ibot+1) trans(i,j,1:ibot+1) = gb_trans(1:ibot+1) diss(i,j,1:ibot+1) = gb_diss(1:ibot+1) radf(i,j,1:ibot+1) = gb_radf(1:ibot+1) sh(i,j,1:ibot+1) = gb_sh(1:ibot+1) sm(i,j,1:ibot+1) = gb_sm(1:ibot+1) gh(i,j,1:ibot+1) = gb_gh(1:ibot+1) evhc(i,j,1:ibot+1) = gb_evhc(1:ibot+1) slislope(i,j,1:ibot ) = gb_slisl(1:ibot ) qtslope (i,j,1:ibot ) = gb_qtsl (1:ibot ) qxtop (i,j,1:ibot ) = gb_qxtop(1:ibot ) qxbot (i,j,1:ibot ) = gb_qxbot(1:ibot ) qc_new (i,j,1:ibot ) = gb_qc_new(1:ibot ) qa_new (i,j,1:ibot ) = gb_qa_new(1:ibot ) sigmas(is-1+i,js-1+j,1:ibot+1) = gb_sigmas (1:ibot+1) !determine maximum altitude that the very stable pbl !is permitted to operate. this is set to the half !level altitude just beneath the first turbulent !level stbltop(i,j) = 0. topfound = .false. kk = ibot + 1 do while (.not.topfound.and.kk.gt.1) kk = kk - 1 if (gb_isturb(kk).gt.0.5) topfound = .true. enddo stbltop(i,j) = gb_z_half(kk+1) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' stbltop = ',stbltop(i,j), ' meters' write (dpu,'(a)') ' ' end if enddo enddo !----------------------------------------------------------------------- ! ! diagnose cloud fraction and condensate tendencies ! !---------------------------------------------------------------- ! If the interface at the top or the bottom of the layer is part ! of a turbulent layer then set the eddytau to the minimum ! timescale of the two. ! do k = 1, nlev tauinvtmp(:,:,k) = max(tauinv(:,:,k),tauinv(:,:,k+1)) enddo where (tauinvtmp .gt. 1.e-10 .and. isturb .gt. 0.5) eddytau = max ( 1./tauinvtmp, mineddytau ) elsewhere eddytau = missing_value endwhere qaturb(is:ie,js:je,:) = qa_new qcturb(is:ie,js:je,:) = qc_new tblyrtau(is:ie,js:je,:) = eddytau !----------------------------------------------------------------------- ! ! blend in monin-obukhov similarity theory mixing coefficients at ! interface levels which are inside the surface layer ! sfclyr_h = frac_sfclyr*pblh sfclyr_h_max = maxval(sfclyr_h) kk = nlev do k = 2, nlev if (minval(z_half_ag(:,:,k)) < sfclyr_h_max) then kk = k exit end if end do k_m_mo = 0.0 k_t_mo = 0.0 call mo_diff(z_half_ag(:,:,kk:nlev), u_star, b_star, & k_m_mo(:,:,kk:nlev), k_t_mo(:,:,kk:nlev)) do k = kk, nlev where(z_half_ag(:,:,k) < sfclyr_h(:,:) .and. & z_half_ag(:,:,k) > 0.) k_t(:,:,k) = k_t_mo(:,:,k) k_m(:,:,k) = k_m_mo(:,:,k) endwhere enddo !----------------------------------------------------------------------- ! ! Diagnostics ! if ( id_fq_cv_int > 0 .or. id_fq_cv_top > 0 .or. id_n2 > 0 .or. & id_fq_turb > 0 .or. id_diss > 0 .or. id_s2 > 0 .or. & id_fq_cv_bot > 0 .or. id_fq_st > 0 .or. id_ri > 0 .or. & id_eddytau > 0 .or. id_tauinv > 0 .or. & id_leng > 0 .or. id_qaedt > 0 .or. & id_bprod > 0 .or. id_sprod > 0 .or. & id_trans > 0 .or. id_qcedt > 0 .or. & id_sigmas > 0 ) then mask3(:,:,1:(nlev+1)) = 1. if (present(kbot)) then where (z_half_ag < 0.) mask3(:,:,:) = 0. end where endif if ( id_fq_cv_int > 0 ) then where (turbtype(2,:,:,:) .eq. 1) tmpdat = 1. elsewhere tmpdat = 0. end where used = send_data (id_fq_cv_int,tmpdat,time, is, js, 1,& rmask=mask3 ) end if if ( id_fq_cv_top > 0 ) then where (turbtype(4,:,:,:) .eq. 1) tmpdat = 1. elsewhere tmpdat = 0. end where used = send_data (id_fq_cv_top,tmpdat,time, is, js, 1,& rmask=mask3 ) end if if ( id_fq_cv_bot > 0 ) then where (turbtype(3,:,:,:) .eq. 1) tmpdat = 1. elsewhere tmpdat = 0. end where used = send_data (id_fq_cv_bot,tmpdat,time, is, js, 1,& rmask=mask3 ) end if if ( id_fq_st > 0 ) then where (turbtype(1,:,:,:) .eq. 1) tmpdat = 1. elsewhere tmpdat = 0. end where used = send_data ( id_fq_st, tmpdat, time, is, js, 1, & rmask=mask3 ) end if if ( id_fq_turb > 0 ) then used = send_data ( id_fq_turb, isturb, time, is, js,1,& rmask=mask ) end if if ( id_n2 > 0 ) then used = send_data ( id_n2, n2, time, is, js, 1, & rmask=mask3 ) end if if ( id_s2 > 0 ) then used = send_data ( id_s2, s2, time, is, js, 1, & rmask=mask3 ) end if if ( id_ri > 0 ) then used = send_data ( id_ri, ri, time, is, js, 1, & rmask=mask3 ) end if if ( id_leng > 0 ) then used = send_data ( id_leng, leng, time, is, js, 1, & rmask=mask3 ) end if if ( id_bprod > 0 ) then used = send_data ( id_bprod, bprod, time, is, js, 1, & rmask=mask3 ) end if if ( id_sprod > 0 ) then used = send_data ( id_sprod, sprod, time, is, js, 1, & rmask=mask3 ) end if if ( id_trans > 0 ) then used = send_data ( id_trans, trans, time, is, js, 1, & rmask=mask3 ) end if if ( id_diss > 0 ) then used = send_data ( id_diss, diss, time, is, js, 1, & rmask=mask3 ) end if if ( id_radf > 0 ) then used = send_data ( id_radf, radf, time, is, js, 1, & rmask=mask3 ) end if if ( id_sh > 0 ) then used = send_data ( id_sh, sh, time, is, js, 1, & rmask=mask3 ) end if if ( id_sm > 0 ) then used = send_data ( id_sm, sm, time, is, js, 1, & rmask=mask3 ) end if if ( id_gh > 0 ) then used = send_data ( id_gh, gh, time, is, js, 1, & rmask=mask3 ) end if if ( id_evhc > 0 ) then used = send_data ( id_evhc, evhc, time, is, js, 1, & rmask=mask3 ) end if if ( id_tauinv > 0 ) then used = send_data ( id_tauinv, tauinv, time, is, js, 1,& rmask=mask3 ) end if if ( id_eddytau > 0 ) then used = send_data ( id_eddytau, eddytau, time, is, js, & 1, rmask=mask ) end if if ( id_sigmas > 0 ) then used = send_data ( id_sigmas, sigmas(is:ie,js:je,:), & time, is, js, 1, rmask=mask3 ) end if if (id_qaedt > 0) then tmpdat = 0.0 do k = 1, nlev tmpdat(:,:,k) = isturb(:,:,k)*qa_new(:,:,k) enddo used = send_data ( id_qaedt, tmpdat(:,:,1:nlev), & time, is, js, 1, rmask=mask ) end if if (id_qcedt > 0) then tmpdat = 0.0 do k = 1, nlev tmpdat(:,:,k) = isturb(:,:,k)*qc_new(:,:,k) enddo used = send_data ( id_qcedt, tmpdat(:,:,1:nlev), & time, is, js, 1, rmask=mask ) end if end if ! do diagnostics if !----------------------------------------------------------------------- ! ! close edt output file if data was written for this window !!RSH if (ipt .gt. 0 .and. jpt .gt. 0) call Close_File (dpu) !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine edt ! !======================================================================= !======================================================================= ! ! subroutine edt_end ! ! ! this subroutine writes out the restart field ! subroutine edt_end() !----------------------------------------------------------------------- ! ! variables ! ! -------- ! internal ! -------- ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! write out restart file ! call save_restart(edt_restart) !----------------------------------------------------------------------- ! ! close edt output file if data was written for this window if (do_print ) call Close_File (dpu) !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! subroutine end ! module_is_initialized = .false. end subroutine edt_end ! !======================================================================= !======================================================================= !======================================================================= ! ! ! Grenier-Bretherton turbulence subroutines follow. These include: ! ! ! sfdiag diagnoses the fraction of each interval ! between layer midpoints which is saturated; ! this fraction is used to calculate the ! proper buoyancy frequency, n2. ! ! trbintd calculates the subgrid vertical variability ! to sl, the liquid water static energy, and ! qt, the total water specific humidity. ! ! also calculates the richardson number ri ! from n2: ! ! n2 = ch*dslidz + cm*dqtdz ! ! where ! ! dslidz and dqtdz are the vertical slopes ! of sl and qt in the interior of the grid ! box and ch and cm are thermodynamic ! coefficients which are solely functions of ! temperature and pressure, but have very ! different values in saturated and unsaturated ! air. ! ! caleddy main subroutine which calculates diffusivity ! coefficients ! ! exacol determines whether the column has adjacent ! regions where Ri < 0 (unstable layers or ULs) ! and determine the indices kbase, ktop which ! delimit these unstable layers : ! ! ri(kbase) > 0 and ri(ktop) > 0, ! but ri(k) < 0 for ktop < k < kbase. ! ! zisocl solves for mean tke , W, Sh, and Sm for ! each convective layer. here, ! ! W = leng**2 * (-Sh*n2 + Sm*s2), where ! ! leng = mixing length, ! s2 = shear vector magnitude ! Sh = Galperin stab. function for buoyancy ! Sm = Galperin stab. function for momentum ! ! in addition, it merges adjacent convective ! layers when the intervening stable layers ! would not consume too much tke for the ! combined convective layer ! ! !======================================================================= !======================================================================= !======================================================================= subroutine sfdiag (qsl, esl, dqsldtl, hleff, qt, qtslope,sli, slislope,& p_full, z_full, p_half, z_half, sdevs, qxtop, qxbot,& qltop, sfuh, sflh, qc_new, qa_new, sfi) !----------------------------------------------------------------------- ! ! Purpose: ! Interface for computing cloud variables for use by turb. scheme ! ! Authors: B. Stevens and C. Bretherton (August 2000) ! !----------------------------------------------------------------------- ! ! variables ! ! ----- ! input ! ----- ! field 1-d arrays on model full levels, reals dimensioned ! (1:nlev), index running from top of atmosphere to bottom ! ! qsl saturation spec. humidity at the liquid-ice ! water temperature (kg water/kg air) ! esl saturation vapor pressure at the liquid-ice ! water temperature (Pa) ! dqsldtl temperature derivative of qsl (kg water/kg air/K) ! hleff effective latent heat of vaporization (J/kg) ! qt total water specific humidity (kg water/kg air) ! qtslope qt slope wrt pressure in thermo layer (kg/kg/Pa) ! sli ice-liq water static energy (J/kg) ! slislope sli slope wrt pressure in thermo layer (J/kg/Pa) ! p_full pressure (Pa) ! z_full height of full level above the surface (m) ! ! the following fields are on the model half levels, ! dimension(1:nlev+1) ! ! p_half pressure at half levels (Pa) ! z_half height of half model levels above the surface (m) ! sdevs standard deviation of water perturbation ! (kg water/kg air) ! ! ------ ! output ! ------ ! ! qxtop saturation excess at top of layer (kg wat/kg air) ! qxbot saturation excess at bottom of layer (kg wat/kg air) ! qltop cloud liquid at top of layer (kg wat/kg air) ! sfuh saturated fraction in upper half-layer ! sflh sflh saturated fraction in lower half-layer ! qc_new thermodynamically diagnosed condensate value ! qa_new thermodynamically diagnosed cloud fraction ! ! the following fields are on the model half levels, ! dimension(1:nlev+1) ! ! sfi interfacial saturated fraction ! real, intent(in) , dimension(:) :: qsl, esl, dqsldtl, hleff real, intent(in) , dimension(:) :: qt, qtslope, sli, slislope real, intent(in) , dimension(:) :: p_full, z_full real, intent(in) , dimension(:) :: p_half, z_half, sdevs real, intent(out), dimension(:) :: qxtop, qxbot, qltop, sfuh, sflh real, intent(out), dimension(:) :: qc_new, qa_new real, intent(out), dimension(:) :: sfi integer :: k ! vertical index integer :: nlev ! number of vertical levels real :: slitop, slibot ! sli at top/bot of layer real :: qttop , qtbot ! qt at top/bot of layer real :: qsltop, qslbot ! qsl at top/bot of layer real :: tlitop, tlibot ! liq wat temp at top/bot of layer real :: qxm ! sat excess at midpoint real :: qlm, qlbot ! liq wat at midpoint, and bottom real :: tlim ! tli at midpoint real :: dqsldp ! pressure derivative of qsl real :: sigmasf ! sdevs on full model levels real :: acoef ! 1./(1+L*dqsldT/cp_air) !----------------------------------------------------------------------- ! ! code ! nlev = size(p_full,1) sfi = 0.0 sfuh = 0.0 sflh = 0.0 qc_new = 0.0 qa_new = 0.0 qltop = 0.0 if (column_match) then write (dpu,'(a)') '-----------------------------------------'//& '-------------' write (dpu,'(a)') ' ENTERING SFDIAG ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k tlim dqsldp slitop tlitop qttop q'//& 'sltop qxtop qltop qxm qlm' write (dpu,'(a)') ' (K)(g/kg/100mb)(K) (K) (g/kg) ('//& 'g/kg) (g/kg) (g/kg) (g/kg) (g/kg)' write (dpu,'(a)') '-----------------------------------------'//& '----------------------------------' write (dpu,'(a)') ' ' 21 format(1X,i2,1X,f6.2,1X,f7.2,1X,8(f6.2,1X)) end if do k = 2,nlev !----------------------------------------------------------- ! calculate midpoint liquid-ice water temperature tlim = (sli(k) - grav*z_full(k))/cp_air !----------------------------------------------------------- ! calculate dqsl/dp dqsldp = -1.*qsl(k)*qsl(k)/d622/esl(k) !----------------------------------------------------------- ! Compute saturation excess at top and bottom of layer k !extrapolation to layer top slitop = sli(k) + slislope(k) * (p_half(k) - p_full(k)) qttop = qt (k) + qtslope (k) * (p_half(k) - p_full(k)) slitop = min ( max ( 0.25*sli(k), slitop), 4.*sli(k)) qttop = min ( max ( 0.25* qt(k), qttop ), 4.* qt(k)) tlitop = (slitop - grav*z_half(k))/cp_air qsltop = qsl(k) + dqsldp * (p_half(k) - p_full(k)) + & dqsldtl (k) * (tlitop - tlim ) qsltop = min ( max ( 0.25*qsl(k), qsltop ) , 4.*qsl(k) ) qxtop(k) = qttop - qsltop qltop(k) = max( 0., qxtop(k)/( 1.+ hleff(k)*dqsldtl(k)/cp_air)) !extrapolation to layer bottom slibot = sli(k) + slislope(k) * (p_half(k+1) - p_full(k)) qtbot = qt (k) + qtslope (k) * (p_half(k+1) - p_full(k)) slibot = min ( max ( 0.25*sli(k), slibot), 4.*sli(k)) qtbot = min ( max ( 0.25* qt(k), qtbot ), 4.* qt(k)) tlibot = (slibot - grav*z_half(k+1))/cp_air qslbot = qsl(k) + dqsldp * (p_half(k+1) - p_full(k))+ & dqsldtl (k) * (tlibot - tlim ) qslbot = min ( max ( 0.25*qsl(k), qslbot ) ,4.*qsl(k) ) qxbot(k) = qtbot - qslbot qlbot = max( 0., qxbot(k)/( 1.+hleff(k)*dqsldtl(k)/cp_air ) ) !----------------------------------------------------------- ! Compute saturation excess at midpoint of layer k qxm = qxtop(k) + (qxbot(k)-qxtop(k))* & (p_full(k)-p_half(k))/(p_half(k+1) - p_half(k)) qlm = max( 0., qxm / ( 1. + hleff(k)*dqsldtl(k)/cp_air ) ) if (column_match) then write(dpu,21) k,tlim,1000.*dqsldp*100.*100.,slitop/cp_air, & tlitop,1000.*qttop,1000.*qsltop,1000.*qxtop(k), & 1000.*qltop(k),1000.*qxm, 1000.*qlm write(dpu,21) k,tlim,1000.*dqsldp*100.*100.,slibot/cp_air, & tlibot,1000.*qtbot,1000.*qslbot,1000.*qxbot(k), & 1000.*qlbot,1000.*qxm, 1000.*qlm write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! ! TWO WAYS OF CALCULATING SATURATION FRACTION AND QA AND QC ! ! ! FIRST WAY: USE GAUSSIAN CLOUD MODEL ! ! if (do_gaussian_cloud) then !calculate sigmas on model full levels sigmasf = 0.5 * (sdevs(k) + sdevs(k+1)) acoef = 1. / ( 1. + hleff(k)*dqsldtl(k)/cp_air ) !sigmasf = max ( sigmasf, acoef*mesovar*qsl(k) ) sigmasf = acoef*mesovar*qsl(k) call gaussian_cloud(qxtop(k), qxm, qxbot(k), & acoef, sigmasf, qa_new(k), & qc_new(k),sfuh(k), sflh(k)) !------------------------------------------------------ ! Combine with sflh (still for layer k-1) to get ! interface layer saturation fraction ! ! N.B.: ! ! if sfuh(k)>sflh(k-1),sfi(k) = sflh(k-1) ! if sfuh(k) 0 ! or vice versa sfuh(k) = max(qxtop(k),qxm) / abs(qxtop(k) - qxm) end if !------------------------------------------------------ ! Combine with sflh (still for layer k-1) to get ! interfac layer sat frac ! ! N.B.: ! ! if sfuh(k)>sflh(k-1),sfi(k) = sflh(k-1) ! if sfuh(k) 0 or vice versa sflh(k) = max(qxbot(k),qxm) / abs(qxbot(k) - qxm) end if !------------------------------------------------------ !Compute grid volume mean condensate and cloud fraction qc_new(k) = 0.5 * ( sfuh(k) * 0.5 * (qltop(k) + qlm) )& + 0.5 * ( sflh(k) * 0.5 * (qlbot + qlm) ) qa_new(k) = 0.5 * ( sfuh(k) + sflh(k) ) end if end do !set surface saturated fraction equal to sflh(nlev) sfi(nlev+1) = sflh(nlev)! Sat frac in lowest half-layer. if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' k sfuh sflh sfi qc_n'//& 'ew qa_new' write (dpu,'(a)') ' (g/kg)' write (dpu,'(a)') '-----------------------------------------'//& '------------' write (dpu,'(a)') ' ' do k = nlev-n_print_levels,nlev write(dpu,22) k,sfuh(k),sflh(k),sfi(k),1000.*qc_new(k), & qa_new(k) enddo 22 format(1X,i2,1X,5(f9.4,1X)) write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine sfdiag ! !======================================================================= !======================================================================= ! ! subroutine trbintd ! subroutine trbintd (t, qv, qt, qc, sli, sliv, u, v, z_full, z_half, & p_full, p_half, sdevs, qsl, esl, hleff, dqsldtl, & slislope, qtslope, qxtop, qxbot, qltop, sfuh,sflh,& qc_new, qa_new, chu, chs, cmu, cms, n2, s2, ri) !----------------------------------------------------------------------- ! ! Purpose: ! time dependent initialization ! ! Method: ! Diagnosis of variables that do not depend on mixing assumptions or ! PBL depth. ! ! Authors: B. Stevens (extracted from pbldiff, August, 2000) ! C. Bretherton (Dec. 2000) ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! variables ! ! ----- ! input ! ----- ! field 1-d arrays on model full levels, reals dimensioned ! (1:nlev), index running from top of atmosphere to ! bottom ! ! t temperature (K) ! qv water vapor spec. humidity (kg vapor/kg air) ! qt total water specific humidity (kg water/kg air) ! qc condensed water spec. humidity (kg cond/kg air) ! sli ice-liq water static energy (J/kg) ! sliv ice-liq water virtual static energy (J/kg) ! u zonal wind (m/s) ! v meridional wind (m/s) ! z_full height of full level above the surface (m) ! p_full pressure (Pa) ! sdevs standard deviation of water perturbation s ! (kg water/kg air) ! qsl saturation spec. humidity at the liquid-ice ! water temperature (kg water/kg air) ! esl saturation vapor pressure at the liquid-ice ! water temperature (Pa) ! hleff effective latent heat (J/kg condensate) ! dqsldtl temperature derivative of qsl (kg water/kg air/K) ! ! the following fields are on the model half levels, ! dimension(1:nlev+1) ! p_half pressure at half levels (Pa) ! z_half height of half model levels above the surface (m) ! ! ------ ! output ! ------ ! ! slislope sli slope wrt pressure in thermo layer (J/kg/Pa) ! qtslope qt slope wrt pressure in thermo layer (kg/kg/Pa) ! qxtop saturation excess at the top of the layer ! (kg wat/kg air) ! qxbot saturation excess at the bottom of the layer ! (kg wat/kg air) ! qltop liquid water at top of thermo layer (kg condensate/ ! kg air) ! sfuh saturated fraction in upper half-layer ! sflh sflh saturated fraction in lower half-layer ! qc_new thermodynamically defined cloud condensate (kg/kg) ! qa_new thermodynamically defined cloud fraction ! ! the following fields are defined on model half levels, ! dimension(1:nlev+1) ! ! chu heat var. coef for dry states (1/m) ! chs heat var. coef for sat states (1/m) ! cmu moisture var. coef for dry states (kg/kg)*(m/s*s) ! cms moisture var. coef for sat states (kg/kg)*(m/s*s) ! n2 moist squared buoyancy freq (1/s*s) ! s2 squared deformation, or shear vector mag. (1/s*s) ! ri gradient Richardson number ! ! real, intent(in), dimension (:) :: t, qv, qt, qc, sli, sliv, u, v, & z_full, p_full, qsl, esl, hleff, & dqsldtl, sdevs, z_half, p_half real, intent(out), dimension (:) :: slislope, qtslope, qxtop, qxbot, & sfuh, sflh, qc_new, qa_new, qltop, & chu, chs, cmu, cms, n2, s2, ri ! internal variables integer :: kdim ! # of levels in the vertical integer :: k, km1, kp ! level indexes real :: rdz ! 1 / (delta z) between midpoints real :: dslidz ! delta sli / delta z at interface real :: dqtdz ! delta qt / delta z at interface real :: ch ! sfi weighted ch at the interface real :: cm ! sfi weighted cm at the interface real :: product ! temporary variable real :: dslidp_a ! sli slope across interface above real :: dqtdp_a ! qt slope across interface above real :: dslidp_b ! sli slope across interface below real :: dqtdp_b ! qt slope across interface below real, dimension(size(t,1)) :: bfact ! buoyancy factor in n2 calculation real, dimension(size(t,1)+1) :: sfi ! saturated fraction at interfaces !----------------------------------------------------------------------- ! ! code ! kdim = size(t,1) !----------------------------------------------------------------------- ! ! Thermodynamic coefficients for buoyancy flux - these ! are calculated at midpoints; they will be averaged to interfaces, ! where they will ultimately be used. At the surface, the coeff- ! icients are taken from the lowest midpoint. ! ! These formulas come from the following expression ! ! grav* tv' / tv = ch * sli' + cm * qt' ! ! chu and cmu are the values for unsaturated air, whereas ! chs and cms are the values for saturated air. bfact = grav/(t*(1.+zvir*qv - qc)) chu(1:kdim) = (1. + zvir*qt)*bfact/cp_air chs(1:kdim) = ( (1. + (1. + zvir)*dqsldtl*t) / & (1. + (hleff*dqsldtl/cp_air) ) ) * bfact/cp_air cmu(1:kdim) = zvir * bfact * t cms(1:kdim) = hleff * chs(1:kdim) - bfact * t chu(kdim+1) = chu(kdim) chs(kdim+1) = chs(kdim) cmu(kdim+1) = cmu(kdim) cms(kdim+1) = cms(kdim) if (column_match) then write (dpu,'(a)') '--------------------------------------------' write (dpu,'(a)') ' ENTERING TRBINTD ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k Tv bfact ' write (dpu,'(a)') ' (K) (m/K/s**2)' write (dpu,'(a)') '------------------------' write (dpu,'(a)') ' ' do k = kdim-n_print_levels,kdim write(dpu,17) k,t(k)*(1.+zvir*qv(k) - qc(k)),bfact(k) end do 17 format(1X,i2,1X,2(f9.4,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' k cp_air*chu cp_air*chs cmu cms' write (dpu,'(a)') ' (m/K/s**2)(m/K/s**2) (m/s**2) (m/s**2)' write (dpu,'(a)') '-------------------------------------------' write (dpu,'(a)') ' ' do k = kdim+1-n_print_levels,kdim+1 write(dpu,20) k,cp_air*chu(k),cp_air*chs(k),cmu(k),cms(k) enddo 20 format(1X,i2,1X,4(f9.4,1X)) write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! Compute slopes in conserved variab. sl, qt within thermo layer k. ! a indicates the 'above' gradient from layer k-1 to layer k and ! b indicates the 'below' gradient from layer k to layer k+1. ! We take the smaller (in absolute value) of these gradients ! as the slope within layer k. If they have opposite sign, gradient ! in layer k is taken to be zero. ! ! Slopes at endpoints determined by extrapolation slislope(kdim) = (sli (kdim) - sli (kdim-1))/ & (p_full(kdim) - p_full(kdim-1)) qtslope (kdim) = (qt (kdim) - qt (kdim-1))/ & (p_full(kdim) - p_full(kdim-1)) slislope(1) = (sli (2) - sli (2) )/ & (p_full(2) - p_full(1) ) qtslope (1) = (qt (2) - qt (2) )/ & (p_full(2) - p_full(1) ) dslidp_b = slislope(1) dqtdp_b = qtslope (1) do k = 2, kdim-1 kp = k + 1 dslidp_a = dslidp_b dqtdp_a = dqtdp_b dslidp_b = (sli(kp)-sli(k))/(p_full(kp)-p_full(k)) dqtdp_b = (qt (kp)-qt (k))/(p_full(kp)-p_full(k)) product = dslidp_a*dslidp_b if (product .le. 0.) then slislope(k) = 0. else if (product.gt.0. .and. dslidp_a.lt.0.) then slislope(k) = max(dslidp_a,dslidp_b) else if (product.gt.0. .and. dslidp_a.gt.0.) then slislope(k) = min(dslidp_a,dslidp_b) end if product = dqtdp_a*dqtdp_b if (product .le. 0.) then qtslope (k) = 0. else if (product.gt.0. .and. dqtdp_a.lt.0.) then qtslope (k) = max(dqtdp_a,dqtdp_b) else if (product.gt.0. .and. dqtdp_a.gt.0.) then qtslope (k) = min(dqtdp_a,dqtdp_b) end if end do ! loop over k if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k slislope qtslope ' write (dpu,'(a)') ' (K/100mb)(g/kg/100mb)' write (dpu,'(a)') '-------------------------' write (dpu,'(a)') ' ' do k = kdim-n_print_levels,kdim write(dpu,17) k,slislope(k)*100.*100./cp_air, & 1000.*qtslope(k)*100.*100. enddo write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! Compute saturation fraction in the interfacial layers for use in ! buoyancy flux computation. call sfdiag(qsl, esl, dqsldtl, hleff, qt, qtslope,sli, slislope,& p_full, z_full, p_half, z_half, sdevs, qxtop, qxbot,& qltop, sfuh, sflh, qc_new, qa_new, sfi) !----------------------------------------------------------------------- ! ! Compute shear squared (s2), squared buoyancy frequency (n2) and ! Ri. For the n2 calculation use gradients of sl and qt, weighted ! according to sfi, the fraction of the interfacial layer that is ! saturated. ! ! This loop has to be done in increasing levels to avoid over- ! writing chu, etc. arrays to interface values before we are done ! with their midpoint values. ! ! Note that n2,s2,ri are set to zero at interfaces k = 1 and ! k = kdim + 1, where they are not used. n2(kdim+1) = 0. s2(kdim+1) = 0. ri(kdim+1) = 0. if (column_match) then write (dpu,'(a)') '---------------------------------------------' write (dpu,'(a)') ' IN TRBINTD' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k rdz dslidz dqtdz chu chs'//& ' ch cmu cms cm' write (dpu,'(a)') ' (1/km) (K/km) (g/kg/km) (m/K/s'//& '**2) (m/s**2)' write (dpu,'(a)') '-----------------------------------------'//& '----------------------------------' write (dpu,'(a)') ' ' 22 format(1X,i2,1X,9(f7.3,1X)) end if do k = kdim, 2, -1 km1 = k - 1 rdz = 1. / (z_full(km1) - z_full(k)) dslidz = (sli(km1) - sli(k)) * rdz dqtdz = (qt (km1) - qt (k)) * rdz chu(k) = (chu(km1) + chu(k))*0.5 chs(k) = (chs(km1) + chs(k))*0.5 cmu(k) = (cmu(km1) + cmu(k))*0.5 cms(k) = (cms(km1) + cms(k))*0.5 ch = chu(k)*(1.-sfi(k)) + chs(k)*sfi(k) cm = cmu(k)*(1.-sfi(k)) + cms(k)*sfi(k) n2(k) = ch*dslidz + cm*dqtdz s2(k) = ((u(km1)-u(k))**2 + (v(km1)-v(k))**2)*(rdz**2) s2(k) = max(ntzero,s2(k)) ri(k) = n2(k) / s2(k) if (column_match) then write(dpu,22) k,1000.*rdz,1000.*dslidz/cp_air, & 1000.*1000.*dqtdz,cp_air*chu(k),cp_air*chs(k),cp_air*ch, & cmu(k),cms(k),cm end if end do n2(1) = 0. s2(1) = 0. ri(1) = 0. if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k sqrt(n2) sqrt(s2) ri' write (dpu,'(a)') ' (1/hr) (1/hr) ' write (dpu,'(a)') '---------------------------------' write (dpu,'(a)') ' ' do k = kdim-n_print_levels,kdim write(dpu,23) k,n2(k)*sqrt(abs(n2(k)))*3600./max(small, & abs(n2(k))),sqrt(s2(k))*3600.,ri(k) enddo 23 format(1X,i2,1X,3(f9.3,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') '--------------------------------------------' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine trbintd ! !======================================================================= !======================================================================= ! ! subroutine exacol ! subroutine exacol(ri, bflxs, ktop, kbase, ncvfin) !----------------------------------------------------------------------- ! ! object : determine whether the column has adjacent regions where ! Ri < 0 (unstable layers or ULs) and determine the indices ! kbase, ktop which delimit these unstable layers : ! ri(kbase) > 0 and ri(ktop) > 0, but ! ri(k) < 0 for ktop < k < kbase. ! ! author : H. Grenier 05/2000, ! C. Bretherton 08/2000 ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! real, intent(in), dimension(:) :: ri ! Moist gradient Ri. # real, intent(in) :: bflxs ! Surface buoyancy flux ! (m*m)/(s*s*s) integer, intent(out), dimension(:) :: kbase ! vertical index of UL base integer, intent(out), dimension(:) :: ktop ! vertical index of UL top integer, intent(out) :: ncvfin ! number of ULs ! ! internal variables ! integer :: k, kdim, ncv real :: riex(size(ri,1)) ! Column Ri profile extended to surface ! by taking riex > rimaxentr for bflxs < 0 ! riex < rimaxentr for bflxs > 0. !----------------------------------------------------------------------- ! ! Initialize variables ncvfin = 0 kdim = size(ri,1) - 1 ktop(:) = 0 kbase(:) = 0 riex(1:kdim) = ri(1:kdim) riex(kdim+1) = rimaxentr-bflxs ! Allows consistent treatment ! of surface with other intrfcs. ncv = 0 !----------------------------------------------------------------------- ! ! Work upward from surf interface k = kdim+1 do while ( k.gt.2 ) if (riex(k) .lt. rimaxentr) then !------------------------------------------------------ ! ! A new convective layer has been found. ! Define kbase as interface below first unstable one ! then decrement k until top unstable level is found. ! Set ktop to the first interface above unstable layer. ncv = ncv + 1 kbase(ncv) = min(k+1,kdim+1) do while (riex(k) .lt. rimaxentr .and. k.gt.2) k = k-1 end do ktop(ncv) = k else !------------------------------------------------------ ! ! Keep on looking for a CL. k = k-1 end if end do !-------------------------- ! ! Set total number of CLs ncvfin = ncv if (column_match) then write (dpu,'(a)') '--------------------------------------------' write (dpu,'(a)') ' IN EXACOL' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k riex' write (dpu,'(a)') '-------------' write (dpu,'(a)') ' ' do k = kdim+1-n_print_levels,kdim+1 write(dpu,28) k,riex(k) enddo 28 format(1X,i2,1X,f9.3) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a,i3)') ' ncvfin = ',ncvfin write (dpu,'(a)') ' ' write (dpu,'(a)') 'ncv ktop kbase' write (dpu,'(a)') '--------------------' write (dpu,'(a)') ' ' do ncv = 1, ncvfin write(dpu,29) ncv, ktop(ncv), kbase(ncv) enddo 29 format(1X,i2,3X,i2,5X,i2) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine exacol ! !======================================================================= !======================================================================= ! ! subroutine zisocl ! subroutine zisocl(u_star, bflxs, tkes, zm, ql, zi, n2, s2, ri, ncvfin, & kbase, ktop, belongcv, ebrk, wbrk, ghcl, shcl, smcl, & lbrk) !----------------------------------------------------------------------- ! ! object : find , , , , ktop(ncv), accounting for the ! presence of stably stratified layers inside the convective ! layer(CL, to be defined) but with r2 > ratinv. ! ! Re-arrange the indexing of arrays kbase/ktop if some CLs are ! found to be coupled such that ncv defines the index of each ! CL increasing with height. ! ! author : H. Grenier 05/08/2000 ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! real, intent(in) :: u_star ! friction velocity (m/s) real, intent(in) :: bflxs ! Surface buoyancy flux ! (m*m*m)/(s*s) real, intent(in) :: tkes ! TKE at the surface (m2/s2) real, intent(in), dimension(:) :: zm ! Layer midpoint height (m) real, intent(in), dimension(:) :: ql ! condensate spec. hum. (kg/kg) real, intent(in), dimension(:) :: zi ! Interface height (m) real, intent(in), dimension(:) :: n2 ! Moist squared buoy freq (s-2) real, intent(in), dimension(:) :: s2 ! Shear deformation (s-2) real, intent(in), dimension(:) :: ri ! Gradient Richardson number integer, intent(inout) :: ncvfin ! Total number of CLs integer, intent(inout), dimension(:) :: kbase ! Vert index of CL base integer, intent(inout), dimension(:) :: ktop ! Vert index of CL top logical, intent(out), dimension(:) :: belongcv ! = T if flux level in CL real, intent(out), dimension(:) :: ebrk ! vert ave. of TKE in CL real, intent(out), dimension(:) :: wbrk ! " of W^2 " real, intent(out), dimension(:) :: ghcl ! " of Gh " real, intent(out), dimension(:) :: shcl ! " of Sh " real, intent(out), dimension(:) :: smcl ! " of Sm " real, intent(out), dimension(:) :: lbrk ! CL depth not within entr ! layers ! internal variables logical :: extend ! True if CL is extended in zisocl logical :: bottom ! True if CL base at surface(kb = kdim+1) integer :: ncv ! Index enumerating convective layers in col integer :: incv integer :: k integer :: kdim ! number of full vertical levels integer :: kb ! Local index for kbase integer :: kt ! Local index for ktop integer :: ncvinit ! Value of ncv at routine entrance integer :: cntu ! counts upward no. of merged CLs integer :: cntd ! counts downward " " integer :: kbinc ! Index for incorporating underlying CL integer :: ktinc ! Index for incorporating overlying CL real :: ebar real :: wint real :: dwinc real :: dzinc real :: dwsurf real :: gh real :: sh real :: sm real :: l2n2 ! Vert. integral of l^2N^2 over CL real :: l2s2 ! Vert. integral of l^2S^2 over CL real :: dl2n2 ! Vert. int. of l^2N^2 over incorp. layer real :: dl2s2 ! Vert. int. of l^2S^2 over incorp. layer real :: lint ! CL depth excluding entrainment layers real :: lbulk ! Depth of the convective layer real :: lz ! Turbulent length scale real :: ricl ! Ri Number for the whole convective layer real :: zbot ! Height of CL base real :: l2rat ! Square of ratio of actual to initial CL depth real :: tmpr !----------------------------------------------------------------------- ! ! Initialize variables kdim = size(ql,1) ncv = 1 if (column_match) then write (dpu,'(a)') '-------------------------------------------' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ENTERING ZISOCL ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k lint lz l2n2 l2s'//& '2' write (dpu,'(a)') ' (m) (m) (m3/s2) (m3/'//& 's2)' write (dpu,'(a)') '-----------------------------------------'//& '------' end if !----------------------------------------------------------------------- ! ! Loop over convective layers to see if they need to be extended do while ( ncv .le. ncvfin ) ncvinit = ncv cntu = 0 cntd = 0 kb = kbase(ncv) kt = ktop(ncv) lbulk = zi(kt)-zi(kb) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,i3)') ' ncv = ', ncv write (dpu,'(a,i3)') ' kb = ', kb write (dpu,'(a,i3)') ' kt = ', kt write (dpu,'(a,f14.7,a)') ' zi(kt) = ', zi(kt), ' meters' write (dpu,'(a,f14.7,a)') ' zi(kb) = ', zi(kb), ' meters' write (dpu,'(a,f14.7,a)') ' lbulk = ', lbulk, ' meters' write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! ! Add contribution (if any) from surface interfacial layer ! to turbulent production and lengthscales. If there is ! positive surface buoyancy flux, the CL extends to the ! surface and there is a surface interfacial layer contri- ! bution to W and the CL interior depth. If there is neg- ! ative buoyancy flux, the surface interfacial layer is ! treated as energetically isolated from the rest of the CL ! and does not contribute to the layer-interior W or depth. ! This case also requires a redefinition of lbulk. bottom = kb .eq. kdim+1 if (bottom .and. (bflxs .ge. 0.)) then lint = zm(kdim) dwsurf = (tkes/b1)*zm(kdim) else if (bottom .and. (bflxs .lt. 0.)) then lint = 0. dwsurf = 0. lbulk = zi(kt)-zm(kdim) else lint = 0. dwsurf = 0. end if l2n2 = 0. l2s2 = 0. if (column_match .and. bottom) then write(dpu,30) kdim+1,lint,lz,l2n2,l2s2 30 format(1X,i2,1X,2(f8.4,1X),2(f12.9,1X)) end if !----------------------------------------------------------- ! ! Turbulence contribution from conv layer (CL) interior ! kt < k < kb, which at this point contains only unstable ! interfaces. Based on the CL interior stratification, ! initial guesses at the stability functions are made. If ! there is no CL interior interface, neutral stability is ! assumed for now. if (kt .lt. kb-1) then do k = kb-1, kt+1, -1 lz = lengthscale(zi(k),lbulk) l2n2 = l2n2 + lz*lz*n2(k)*(zm(k-1)-zm(k)) l2s2 = l2s2 + lz*lz*s2(k)*(zm(k-1)-zm(k)) lint = lint + (zm(k-1)-zm(k)) if (column_match) write(dpu,30) k,lint,lz,l2n2, & l2s2 enddo !------------------------------------------------------ ! ! Solve for bulk Sh, Sm, and wint over the CL interior ricl = min(l2n2/l2s2,ricrit) ! actually we should have ! ricl < 0 call galperin(ricl,gh,sh,sm) wint = -sh*l2n2 + sm*l2s2 + dwsurf ebar = b1*wint/lint else !----------------------------------------------------------- ! ! There is no CL interior interface. The only way that ! should happen at this point is if there is upward surface ! buoy flux but no unstable interior interfaces. In that ! case, the surface interface turbulent production terms are ! used as its CL 'interior'. if (bottom) then wint = dwsurf ebar = tkes !use neutral stability fns for layer extension call galperin(0.,gh,sh,sm) else call error_mesg ('edt_mod', & 'no convective layers found although ncv <= ncvfin', & FATAL) endif endif if(column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7)') ' ricrit = ',ricrit write (dpu,'(a,f14.7)') ' ricl = ', ricl write (dpu,'(a,f14.7)') ' gh = ', gh write (dpu,'(a,f14.7)') ' sh = ', sh write (dpu,'(a,f14.7)') ' sm = ', sm write (dpu,'(a,f14.7,a)') ' lbulk = ', lbulk, ' meters' write (dpu,'(a,f14.7,a)') ' wint = ', wint, ' m3/s2' write (dpu,'(a,f14.7,a)') ' dwsurf = ', dwsurf, ' m3/s2' write (dpu,'(a,f14.7,a)') ' ebar = ', ebar, ' m2/s2' write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! ! Try to extend the top of the convective layer. Compute ! possible contributions to TKE production and lengthscale ! were the CL top interfacial layer found by exacol incor- ! porated into the CL interior. extend = .false. ! will become true if CL top is extended dzinc = zm(kt-1)-zm(kt) lz = lengthscale(zi(kt),lbulk) dl2n2 = lz*lz*n2(kt)*dzinc dl2s2 = lz*lz*s2(kt)*dzinc dwinc = -sh*dl2n2 + sm*dl2s2 if (column_match) then write (dpu,'(a)') '------------------------------------'//& '-------' write (dpu,'(a)') ' ' write (dpu,'(a)') ' trying to extend a layer upwards ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k dzinc lz dl2n2 '//& ' dl2s2 dwinc' write (dpu,'(a)') ' (m) (m) (m3/s2) '//& ' (m3/s2) (m3/s2)' write (dpu,'(a)') '------------------------------------'//& '----------------------' write (dpu,'(a)') ' ' write(dpu,31) kt,dzinc,lz,dl2n2,dl2s2,dwinc 31 format(1X,i2,1X,2(f8.4,1X),3(f12.9,1X)) end if !----------------------------------------------------------- ! ! Test for incorporation of top layer kt into CL interior. ! If true, extend CL by incorporating top layers until test ! fails. l2n2 = -max(min(-l2n2,tkemax*lint/(b1*sh)),tkemin*lint/ (b1*sh)) tmpr = -rinc*dzinc*l2n2/(lint+(1-rinc)*dzinc) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') '-dl2n2 = ', -1.*dl2n2, ' m3/s2' write (dpu,'(a,f14.7,a)') ' l2n2 = ', l2n2, ' m3/s2' write (dpu,'(a,f14.7)') ' rinc = ', rinc write (dpu,'(a,f14.7,a)') ' tmpr = ', tmpr, ' m3/s2' !write (dpu,'(a,f14.7)') ' will layer be extended (-dl2'//& ! 'n2 .gt. tmpr)? ', real( -dl2n2 .gt. tmpr ) write (dpu,'(a)') ' ' end if do while (-dl2n2 .gt. tmpr) !------------------------------------------------------ ! Add contributions from layer kt to interior length- ! scale/TKE prod lint = lint + dzinc wint = wint + dwinc l2n2 = l2n2 + dl2n2 l2n2 = -max(min(-l2n2,tkemax*lint/(b1*sh)),tkemin*lint/(b1*sh)) l2s2 = l2s2 + dl2s2 kt = kt-1 extend = .true. if (kt .eq. 1) then call error_mesg ('edt_mod', & 'trying to extend convective layer at model top',& FATAL) end if !------------------------------------------------------ ! Check for existence of an overlying CL which might ! be merged. If such exists (ktinc > 1), check for ! merging by testing for incorporation of its top ! interior interface into current CL. If no such layer ! exists, ktop(ncv+cntu+1) will equal its default value ! of zero, so ktinc will be 1 and the test kt=ktinc ! will fail. ktinc = ktop(ncv+cntu+1)+1 if (kt .eq. ktinc) then ncvfin = ncvfin - 1 cntu = cntu + 1 end if !------------------------------------------------------ ! Compute possible lengthscale and TKE production ! contributions were layer kt incorporated into CL ! interior. Then go back to top of loop to test for ! incorporation. dzinc = zm(kt-1)-zm(kt) lz = lengthscale(zi(kt),lbulk) dl2n2 = lz*lz*n2(kt)*dzinc dl2s2 = lz*lz*s2(kt)*dzinc dwinc = -sh*dl2n2 + sm*dl2s2 if (column_match) write(dpu,31) kt,dzinc,lz,dl2n2, dl2s2,dwinc !------------------------------------------------------ ! Recalculate tmpr tmpr = -rinc*dzinc*l2n2/(lint+(1-rinc)*dzinc) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') '-dl2n2 = ', -1.*dl2n2, & ' m3/s2' write (dpu,'(a,f14.7,a)') ' l2n2 = ', l2n2, & ' m3/s2' write (dpu,'(a,f14.7)') ' rinc = ', rinc write (dpu,'(a,f14.7,a)') ' tmpr = ', tmpr, & ' m3/s2' !write (dpu,'(a,f14.7,a)') ' will layer be extende'//& ! 'd (-dl2n2 .gt. tmpr)? ', real(-dl2n2 .gt. tmpr) write (dpu,'(a)') ' ' end if end do ! Done with top extension of CL !----------------------------------------------------------- ! Shift indices appropriately if layers have been merged if (cntu .gt. 0) then do incv = 1, ncvfin - ncv kbase(ncv+incv) = kbase(ncv+cntu+incv) ktop(ncv+incv) = ktop(ncv+cntu+incv) end do end if !----------------------------------------------------------- ! ! Extend the CL base if possible. if (column_match .and. .not. bottom) then write (dpu,'(a)') '------------------------------------'//& '-------' write (dpu,'(a)') ' ' write (dpu,'(a)') ' trying to extend a layer downwards' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') 'k dzinc lz dl2n2 '//& ' dl2s2 dwinc' write (dpu,'(a)') ' (m) (m) (m3/s2) '//& ' (m3/s2) (m3/s2)' write (dpu,'(a)') '------------------------------------'//& '---------------------------' write (dpu,'(a)') ' ' end if if (.not. bottom) then !------------------------------------------------------ ! Compute possible contributions to TKE production and ! lengthscale, were the CL base interfacial layer found ! by exacol incorporated into the CL interior. dzinc = zm(kb-1)-zm(kb) lz = lengthscale(zi(kb),lbulk) dl2n2 = lz*lz*n2(kb)*dzinc dl2s2 = lz*lz*s2(kb)*dzinc dwinc = -sh*dl2n2 + sm*dl2s2 if (column_match) write(dpu,31) kb,dzinc,lz,dl2n2,dl2s2,dwinc !------------------------------------------------------ ! Test for incorporation of base layer kb into CL ! interior. If true, extend CL by incorporating base ! layers until test fails. tmpr = -rinc*dzinc*l2n2/(lint+(1-rinc)*dzinc) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') '-dl2n2 = ', -1.*dl2n2, & ' m3/s2' write (dpu,'(a,f14.7,a)') ' l2n2 = ', l2n2, & ' m3/s2' write (dpu,'(a,f14.7)') ' rinc = ', rinc write (dpu,'(a,f14.7,a)') ' tmpr = ', tmpr, & ' m3/s2' !write (dpu,'(a,f14.7,a)') ' will layer be extende'//& ! 'd (-dl2n2 .gt. tmpr)? ', real(-dl2n2 .gt. tmpr) write (dpu,'(a)') ' ' end if do while((-dl2n2.gt. tmpr) .and. (.not. bottom) ) !------------------------------------------------- ! Add contributions from layer kb to interior ! lengthscale/TKE prod lint = lint + dzinc wint = wint + dwinc l2n2 = l2n2 + dl2n2 l2n2 = -max(min(-l2n2,tkemax*lint/(b1*sh)),tkemin& *lint/(b1*sh)) l2s2 = l2s2 + dl2s2 !------------------------------------------------- ! Extend base of CL downward a layer kb = kb+1 extend = .true. !------------------------------------------------- ! Check for existence of an underlying CL which ! might be merged. If such exists (kbinc > 1), ! check for merging by testing for incorporation ! of its top interior interface into current CL. ! Note that this top 'interior' interface could be ! the surface. kbinc = 0 if (ncv .gt. 1) kbinc = ktop(ncv-1)+1 if (kb .eq. kbinc) then !-------------------------------------------- ! We are incorporating interior of CL ncv-1, ! so merge this CL into the current CL. ncv = ncv - 1 ncvfin = ncvfin - 1 cntd = cntd + 1 end if !------------------------------------------------- ! If CL would now reach the surface, check sign of ! surface buoyancy flux. If positive, add contri- ! butions of surface interfacial layer to TKE ! production and lengthscale. If negative, we ! regard the surface layer as stable and do not ! add surface interfacial layer contributions to ! the CL. In either case the surface interface is ! classified as part of the CL for bookkeeping ! purposes (to ensure no base entrainment calcula- ! tion is done). If we are merging with a surface- ! driven CL with no interior unstable interfaces, ! the above code will already have handled the ! merging book-keeping. bottom = kb .eq. kdim+1 if (bottom) then if (bflxs .gt. 0.) then dwsurf = (tkes/b1)*zm(kdim) lint = lint + zm(kdim) end if else !-------------------------------------------- ! Compute possible lengthscale and TKE prod- ! uction contributions were layer kb incor- ! porated into CL interior,then go back to ! top of loop to test for incorporation dzinc = zm(kb-1) - zm(kb) lz = lengthscale(zi(kb),lbulk) dl2n2 = lz*lz*n2(kb)*dzinc dl2s2 = lz*lz*s2(kb)*dzinc dwinc = -sh*dl2n2 + sm*dl2s2 end if if (column_match) write(dpu,31) kb,dzinc,lz,dl2n2,dl2s2,dwinc !------------------------------------------------- ! Recalculate tmpr tmpr = -rinc*dzinc*l2n2/(lint+(1-rinc)*dzinc) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') '-dl2n2 = ',-1.*dl2n2,' m3/s2' write (dpu,'(a,f14.7,a)') ' l2n2 = ', l2n2,' m3/s2' write (dpu,'(a,f14.7)') ' rinc = ', rinc write (dpu,'(a,f14.7,a)') ' tmpr = ', tmpr,' m3/s2' !write (dpu,'(a,f14.7,a)') ' will layer be ex'//& ! 'tended (-dl2n2 .gt. tmpr)? ', real(-dl2n2 & ! .gt. tmpr) write (dpu,'(a)') ' ' end if end do ! for downward extension if (bottom .and. ncv .ne. 1) then call error_mesg ('edt_mod', & 'bottom convective layer not indexed 1',& FATAL) end if end if ! Done with bottom extension of CL !----------------------------------------------------------- ! Shift indices if some layers with N2 < 0 have been found if (cntd .gt. 0) then do incv = 1, ncvfin - ncv kbase(ncv+incv) = kbase(ncvinit+incv) ktop(ncv+incv) = ktop(ncvinit+incv) end do end if !----------------------------------------------------------- ! Sanity check for positive wint. if (wint .lt. 0.) then call error_mesg ('edt_mod', & 'interior avg TKE < 0', FATAL) end if !----------------------------------------------------------- ! Recompute base and top indices, Ri_cl, Sh, Sm, and ! after layer extension if necessary. Ideally, we would ! recompute l2n2 and l2s2 to account for the incorrect lbulk ! used in the computation of lz, but we take the simpler ! approach of simply multiplying the lz's by the ratio of ! the actual PBL depth to lbulk. if (extend) then ktop (ncv) = kt kbase(ncv) = kb zbot = zi(kb) if (bottom .and. (bflxs.lt.0)) zbot = zm(kdim) l2rat = ((zi(kt) - zbot)/lbulk)**2 l2n2 = l2n2*l2rat l2s2 = l2s2*l2rat ricl = min(l2n2/l2s2,ricrit) call galperin(ricl,gh,sh,sm) !------------------------------------------------------ ! It is conceivable that even though the original wint ! was positive, it will be negative after correction. ! In this case, correct wint to be a small positive ! number wint = max(dwsurf + (-sh*l2n2 + sm*l2s2),0.01*wint) end if ! for extend if lbrk(ncv) = lint wbrk(ncv) = wint/lint ebrk(ncv) = b1*wbrk(ncv) ebrk(ncv) = max(min(ebrk(ncv),tkemax),tkemin) ghcl(ncv) = gh shcl(ncv) = sh smcl(ncv) = sm if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,i4)') ' FINAL RESULTS FOR CONVECTIVE LAYER', ncv write (dpu,'(a)') ' ' write (dpu,'(a,i4,i4)') ' ktop(ncv), kbase(ncv) = ', & ktop(ncv), kbase(ncv) write (dpu,'(a,f14.7,a)') ' lbrk(ncv) = ', lbrk(ncv),' m' write (dpu,'(a,f14.7,a)') ' wbrk(ncv) = ', wbrk(ncv),' m2/s2' write (dpu,'(a,f14.7,a)') ' sqrt(ebrk(ncv)) = ', & sqrt(ebrk(ncv)), ' m/s' write (dpu,'(a,f14.7)') ' ghcl(ncv) = ', ghcl(ncv) write (dpu,'(a,f14.7)') ' shcl(ncv) = ', shcl(ncv) write (dpu,'(a,f14.7)') ' smcl(ncv) = ', smcl(ncv) end if !----------------------------------------------------------- ! Increment counter for next CL ncv = ncv + 1 end do ! Loop over convective layers !---------------------------------------------------------------- ! ! set belongcv belongcv(:) = .false. do ncv = 1, ncvfin do k = ktop(ncv), kbase(ncv) belongcv(k) = .true. enddo enddo if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' END OF ZISOCL ' write (dpu,'(a)') ' ' write (dpu,'(a)') '-------------------------------------------' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine zisocl ! !======================================================================= !======================================================================= ! ! subroutine caleddy ! ! ! this is the main program provided by Chris Bretherton and Herve ! Grenier ! subroutine caleddy( u_star, kqfs, khfs, sl, & qt, ql, qa, slv, & u, v, pm, zm, & sfuh, sflh, slslope, qtslope, & qrl, hleff, density, qsl, & dqsldtl, qltop, pi, zi, & chu, chs, cmu, cms, & n2, s2, ri, pblh, & turbtype, kvh, kvm, tke, & leng, bprod, sprod, trans, & diss, isturb, adj_time_inv, sdevs, & radfvec, shvec, smvec, ghvec, & evhcvec) !----------------------------------------------------------------------- ! ! Driver routine to compute eddy diffusion coefficients for ! momentum, moisture, trace constituents and static energy. Uses ! first order closure for stable turbulent layers. For convective ! layers, an entrainment closure is used, coupled to a diagnosis ! of layer-average TKE from the instantaneous thermodynamic and ! velocity profiles. Convective layers are diagnosed by extending ! layers of moist static instability into adjacent weakly stably ! stratified interfaces, stopping if the stability is too strong. ! This allows a realistic depiction of dry convective boundary ! layers with a downgradient approach. ! ! Authors: Herve Grenier, 06/2000, Chris Bretherton 09/2000 ! !----------------------------------------------------------------------- ! ! Variable declarations ! ! Inputs ! real, intent(in) :: u_star ! surface friction velocity real, intent(in) :: kqfs ! kinematic surf constituent ! flux (kg/kg)*(m/s) real, intent(in) :: khfs ! kinematic surface heat flux ! (K*m/s) ! ! the following fields are defined on model full ! levels with dimension (1:nlev) ! real, intent(in), dimension(:) :: sl ! liq water static energy ! cp*T+g*z-L*ql (J/kg) real, intent(in), dimension(:) :: qt ! total water spec. hum. ! (kg water/kg air) real, intent(in), dimension(:) :: ql ! liq water spec. hum. ! (kg condensate/kg air) real, intent(in), dimension(:) :: qa ! cloud fraction (fraction) real, intent(in), dimension(:) :: slv ! liq water virtual static ! energy (J/kg) ! sl*(1 + .608*qt) real, intent(in), dimension(:) :: u ! u wind input (m/s) real, intent(in), dimension(:) :: v ! v wind input (m/s) real, intent(in), dimension(:) :: pm ! midpoint pressures (Pa) real, intent(in), dimension(:) :: zm ! layer midpoint height ! above sfc (m) real, intent(in), dimension(:) :: sfuh ! sat frac in upper half-lyr real, intent(in), dimension(:) :: sflh ! sat frac in lower half-lyr real, intent(in), dimension(:) :: slslope ! sl slope with respect to ! pressure in thermo lyr ! (J/kg/Pa) real, intent(in), dimension(:) :: qtslope ! qt slope with respect to ! pressure in thermo lyr ! (kg water/kg air/Pa) real, intent(in), dimension(:) :: qrl ! LW heating rate (K/s) real, intent(in), dimension(:) :: hleff ! effective latent heat of ! condensation (J/kg) real, intent(in), dimension(:) :: density ! air density (kg/m3) real, intent(in), dimension(:) :: qsl ! saturation specific hum. ! (kg water/kg air) real, intent(in), dimension(:) :: dqsldtl ! temperature derivative of ! qsl (kg water/kg air/K) real, intent(in), dimension(:) :: qltop ! cloud liquid at the top ! of the thermo layer ! ! the following fields are defined on model half levels, ! dimension(1:nlev+1) ! real, intent(in), dimension(:) :: pi ! interface pressures (Pa) real, intent(in), dimension(:) :: zi ! interface height above ! sfc (m) real, intent(in), dimension(:) :: chu ! Unsat sl (heat) coef (1/m) real, intent(in), dimension(:) :: chs ! Sat sl (heat) coef (1/m) real, intent(in), dimension(:) :: cmu ! Unsat qt (moisture) coef ! (kg/kg)*(m/s*s) real, intent(in), dimension(:) :: cms ! Sat qt (moisture) coef ! (kg/kg)*(m/s*s) real, intent(in), dimension(:) :: n2 ! Moist squared buoy freq ! (1/sec**2) real, intent(in), dimension(:) :: s2 ! Squared deformation (s-2) ! (1/sec**2) real, intent(in), dimension(:) :: ri ! gradient Richardson number ! ! Outputs ! real, intent(out) :: pblh ! planetary boundary layer height (m) ! ! the following field is defined on model full levels, ! dimension(1:nlev) ! real, intent(out), dimension(:) :: isturb ! is full layer part of a ! turbulent layer? ! ! the following fields are defined on model half levels, ! dimension(1:nlev+1) ! integer, intent(out), dimension(:,:) :: turbtype! Interface turb. type ! 1 = stable turb ! 2 = CL interior ! 3 = bottom entr intfc ! 4 = upper entr intfc real, intent(out), dimension(:) :: kvh ! diffusivity for heat ! and tracers (m*m/s) real, intent(out), dimension(:) :: kvm ! diffusivity for mom. ! (m*m/s) real, intent(out), dimension(:) :: tke ! turb. kin. energy ! (m*m)/(s*s) real, intent(out), dimension(:) :: leng ! turbulent length scale ! (m) real, intent(out), dimension(:) :: bprod ! Buoyancy production ! (m*m)/(s*s*s) real, intent(out), dimension(:) :: sprod ! shear production ! (m*m)/(s*s*s) real, intent(out), dimension(:) :: trans ! TKE transport ! (m*m)/(s*s*s) real, intent(out), dimension(:) :: diss ! TKE dissipation ! (m*m)/(s*s*s) real, intent(out), dimension(:) :: adj_time_inv ! inverse adjustment ! time for turbulent ! layer (1/sec) real, intent(out), dimension(:) :: sdevs ! std. dev. of water ! perturbation ! (kg water/kg air) ! defined on model half ! levels real, intent(out), dimension(:) :: radfvec ! Buoyancy production ! from lw radiation ! (m*m)/(s*s*s) real, intent(out), dimension(:) :: shvec ! Galperin heat stab. ! fn. (none) real, intent(out), dimension(:) :: smvec ! Galperin mom. stab. ! fn. (none) real, intent(out), dimension(:) :: ghvec ! Galperin stability ! ratio (none) real, intent(out), dimension(:) :: evhcvec ! Evaporative cooling ! entrainment factor ! (none) ! ! Internal variables ! logical :: in_CL ! True if interfaces k,k+1 ! both in same CL logical :: any_stable ! Are there any stable ! turbulent layers? logical :: cloudtop ! Is the interface at ! cloudtop? logical, dimension(size(sl,1)+1) :: belongcv ! True for interfaces ! interior to convective ! layer (CL) logical, dimension(size(sl,1)+1) :: belongst ! True for interfaces ! interior to a stable ! turbulent layer integer :: k ! vertical index integer :: ks ! vertical index integer :: kk ! vertical index integer :: kdim ! number of full vertical levels integer :: ncvfin ! Total number of CL in column integer :: ncvf ! Total number of CL in column prior to ! addition of one layer rad-driven CLs integer :: ncv ! index of current CL integer :: ncvnew ! index of added one layer rad-driven CL integer :: ncvsurf ! if nonzero, index of CL including ! surface integer :: kb, kt ! kbase and ktop for current CL integer, dimension(size(sl,1)+1) :: kbase ! vert. index for base ! interface of CL integer, dimension(size(sl,1)+1) :: ktop ! vert. index for top ! interface of CL real :: bflxs ! Surface buoyancy flux (m2/s3) real :: tkes ! Surface TKE real :: jtzm ! Interface layer thickness atop conv layer (CL) ncv real :: jtsl ! Jump in s_l atop " real :: jtqt ! Jump in q_t atop " real :: jtbu ! Jump in buoyancy atop " real :: jtu ! Jump in zonal wind atop " real :: jtv ! Jump in meridional wind atop " real :: jt2slv ! 2-layer Jump in s_lv atop " real :: radf ! buoy flx jump at cloudtop from lw rad flx div real :: jbzm ! Interface layer thickness at base of CL ncv real :: jbsl ! Jump in s_l at base " real :: jbqt ! Jump in qt at base " real :: jbbu ! Jump in buoyancy at base " real :: jbu ! Jump in zonal wind at base " real :: jbv ! Jump in merid. wind at base " real :: ch ! buoy flux coefs for sl, qt in a half-layer real :: cm ! real :: n2h ! Moist squared buoy freq for a half-layer (s-2) real :: ckh ! Galperin stability function for heat real :: ckm ! Galperin stability function for momentum real :: gh ! Normalised buoyancy production (m*m*m)/(s*s) real :: lbulk ! Depth of turbulent layer (m) real :: trma ! intermediate variables real :: vus real :: vub real :: trmp real :: trmq real :: angle real :: qq real :: rootp real :: evhc ! (1+E) with E = evap. cool. efficiency [nd] real :: vys ! n2h/n2 at upper inversion [nd] real :: vyb ! Same at lower inversion [nd] real :: kentr ! effective entrainment diffusivity we*dz (m*m/s) real :: lwp ! liquid water path in layer kt (kg cond/m2) real :: opt_depth ! optical depth of layer kt real :: radinvfrac ! frac of lw cooling in layer kt put at inv. real :: qtsltmp ! dqt/dz (kg/kg/m) real :: slsltmp ! dsl/dz (J/kg/m) real :: qsltmp ! qsl (kg water/kg air) real :: dqsldtltmp ! dqsldtl (kg water/kg air/K) real :: hlefftmp ! effective latent heat (J/kg water) real :: qstartmp ! q* at upper/lower inversion (kg water/kg air) real :: bstartmp ! b* at upper/lower inversion (m/s2) real :: temperature ! actual temperature (K) real :: tmpfrac ! fraction of cloud at top of convective layer ! exposed to clear air above using maximum ! overlap assumption real, dimension(size(sl,1)+1) :: ebrk,wbrk,lbrk,ghcl,shcl,smcl real, dimension(size(sl,1)+1) :: wcap ! W (m2/s2) real, dimension(size(sl,1)+1) :: rcap ! e/ !----------------------------------------------------------------------- ! ! Initialize to zero outputs needed at all interfaces, but ! calculated only at turbulent interfaces. Only kvh and kvm are ! outputs, the other arrays are zeroed for plotting or diagnostic ! purposes. kdim = size(sl,1) kvh(:) = 0.0 kvm(:) = 0.0 wcap(:) = 0.0 leng(:) = 0.0 rcap(:) = 0.0 bprod(:) = 0.0 sprod(:) = 0.0 trans(:) = 0.0 diss(:) = 0.0 tke(:) = 0.0 turbtype(:,:) = 0 adj_time_inv(:) = missing_value isturb(:) = 0.0 sdevs(:) = missing_value radfvec(:) = 0.0 shvec(:) = missing_value smvec(:) = missing_value ghvec(:) = missing_value evhcvec(:) = missing_value !----------------------------------------------------------------------- ! ! Optional printout if (column_match) then write (dpu,'(a)') '--------------------------------------------' write (dpu,'(a)') ' ENTERING CALEDDY ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') 'checking inputs: ' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' u_star = ',u_star, ' m/s' write (dpu,'(a,f14.7,a)') ' latent heat flux = ',density(kdim)& *hleff(kdim)*kqfs,' W/m2' write (dpu,'(a,f14.7,a)') ' sensible heat flux = ',density(kdim)& *cp_air*khfs,' W/m2' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k sl/cp_air qt qc q'//& 'a slv/cp_air' write (dpu,'(a)') ' (K) (g/kg) (g/kg) '//& ' (K) ' write (dpu,'(a)') '-----------------------------------------'//& '------------' write (dpu,'(a)') ' ' do kk = kdim-n_print_levels,kdim write(dpu,224) kk,sl (kk)/cp_air,1000.*qt(kk),1000.*ql(kk), & qa(kk),slv(kk)/cp_air end do 224 format(1X,i2,1X,5(f9.4,1X)) 24 format(1X,i2,1X,4(f9.4,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' k u v z_full p_full' write (dpu,'(a)') ' (m/s) (m/s) (m) (mb) ' write (dpu,'(a)') '-------------------------------------------' write (dpu,'(a)') ' ' do kk = kdim-n_print_levels,kdim write(dpu,24) kk,u(kk),v(kk),zm(kk),pm(kk)/100. enddo write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k sfuh sflh slslope qtslope' write (dpu,'(a)') ' (K/100mb)(g/kg/100mb)' write (dpu,'(a)') '---------------------------------------------' write (dpu,'(a)') ' ' do kk = kdim-n_print_levels,kdim write(dpu,24) kk,sfuh(kk),sflh(kk),slslope(kk)*100.*100./cp_air& ,1000.*qtslope(kk)*100.*100. enddo write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') 'k rho qrl hleff' write (dpu,'(a)') ' (kg/m3) (K/day) (MJ/kg)' write (dpu,'(a)') '---------------------------------' write (dpu,'(a)') ' ' do kk = kdim-n_print_levels,kdim write(dpu,124) kk,density(kk),qrl(kk)*86400,hleff(kk)/1.e+06 124 format(1X,i2,1X,3(f9.4,1X)) enddo write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k p_half z_half chu chs cmu '//& ' cms sqrt(n2) sqrt(s2) ri' write (dpu,'(a)') ' (mb) (m) (m/K/s**2) (m/s*'//& '*2) (1/hr) (1/hr) ' write (dpu,'(a)') '-----------------------------------------'//& '--------------------------------------' write (dpu,'(a)') ' ' do kk = kdim+1-n_print_levels,kdim+1 write(dpu,26) kk,pi(kk)/100.,zi(kk),cp_air*chu(kk),cp_air*chs(kk), & cmu(kk),cms(kk),n2(kk)*sqrt(abs(n2(kk)))*3600./& max(small,abs(n2(kk))),sqrt(s2(kk))*3600.,ri(kk) enddo 26 format(1X,i2,1X,f7.2,1X,f7.1,1X,2(f7.4,1X),2(f6.3,1X),3(f9.3,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') '-----------------------------------------'//& '-------------------------------------' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! calculate 'surface' (actually lowest half-layer) buoyancy flux. ! ! Note that sensible heat flux in W/m2 = density_air * cp * khfs ! latent heat flux in W/m2 = density_air * L * kqfs ! ! units of bflxs below is m2/s3 ! ! Buoyancy flux in W/m2 = rho * bflxs / ch ! ch = chu(kdim+1)*(1-sflh(kdim)) + chs(kdim+1)*sflh(kdim) cm = cmu(kdim+1)*(1-sflh(kdim)) + cms(kdim+1)*sflh(kdim) bflxs = ch*cp_air*khfs + cm*kqfs bprod(kdim+1) = bflxs !----------------------------------------------------------------------- ! ! Diagnostic surface TKE used in calculating layer-average TKE. ! ! Chris B. note: ! ! Originally tkes was computed from the commented out line below. ! This commented-out line can make tkes small or even negative. ! The replacement has at least equal justifiability, trying to ! represent tke at the surface rather than at zm(kdim). Near the ! surface is where averaging W to get is not really just- ! ifiable. ! ! tkes = (b1*(ustar**3+vonkarm*zm(kdim)*bflxs))**2./3. ! ! [Steve Klein note: ! ! The replacement line appears to create an inconsistency in the ! sense that the buoyancy flux in the lower half of level one does ! not effectively contribute to dwsurf, the contribution to W, ! in the lower half of level one. Hence diagnosing shear production ! and dissipation at the surface is a bit of a ruse.] tkes = max(min(b123*u_star**2,tkemax),tkemin) tke(kdim+1) = tkes diss(kdim+1) = tkes**(3/2)/b1/zm(kdim) sprod(kdim+1) = tkes**(3/2)/b1/zm(kdim) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' ch (surface) = ',cp_air*ch,' m/K/s**2' write (dpu,'(a,f14.7,a)') ' cm (surface) = ',cm,' m/s**2' write (dpu,'(a,f14.7,a)') ' bflxs = ',bflxs,' m2/s3' write (dpu,'(a,f14.7,a)') ' buoyancy flux = ',density(kdim)* & bflxs/ch,' W/m2' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' u_star = ', u_star,' m/s' write (dpu,'(a,f14.7,a)') ' surface tke sqrt = ',sqrt(tkes), & ' m/s' write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! Examine each column and determine whether it is convective. call exacol(ri, bflxs, ktop, kbase, ncvfin) !----------------------------------------------------------------------- ! ! CONVECTIVE LAYER (CL) computations ! ! If some convective layers have been found, determine their bulk ! properties (, , , and indices ktop/bkase for upper/ ! lower inversion ). ncvsurf = 0 if (ncvfin .gt. 0) then call zisocl(u_star, bflxs, tkes, zm, ql, zi, n2, s2, ri, & ncvfin, kbase, ktop, belongcv, ebrk, wbrk, ghcl, & shcl, smcl, lbrk) !------------------------------------------------------------- ! CLs found by zisocl are in order of height, so if any CL ! contains the surface, it will be CL1. if (kbase(1) .eq. kdim+1) ncvsurf = 1 else belongcv(:) = .false. end if if (column_match) write (dpu,'(a,i3)') ' ncvsurf = ', ncvsurf !----------------------------------------------------------------------- ! ! Find single-level radiatively-driven cloud-topped convective ! layers (SRCLs). SRCLs extend through a single thermo layer k, ! with entrainment at interfaces k and k+1 (unless k+1 is the ! surface, in which case surface shear generation contributes to ! the layer-averaged energy). The conditions for an SRCL are: ! 1. cloud at level k ! 2. no cloud at level k+1 (else assuming that some fraction ! of the longwave flux div in layer k is concentrated at ! the top interface is invalid. ! 3. Longwave radiative cooling (shortwave heating is assumed ! uniformly distributed through layer k, so not relevant to ! buoyancy production of TKE) ! 4. Internal stratification n2h of half-layer from level k to ! interface k is unstable using similar method as in sfdiag, ! but applied to internal slopes of sl, qt in layer k. ! 5. Interfaces k, k+1 not both in the same existing convective ! layer. ! 6. k >= 2 ! 7. Ri at interface k > ricrit, otherwise stable turb mixing ! will broadly distribute the cloud top in the vertical, ! preventing localized radiative radiative destabilization ! at the interface height. ncv = 1 ncvf = ncvfin do k = kdim, 2, -1 cloudtop = .false. if (ql(k) .gt. qcminfrac*qsl(k) .and. use_qcmin .and. & ql(k-1).lt. qcminfrac*qsl(k-1)) cloudtop = .true. if (.not.use_qcmin .and. qa(k).gt.qa(k-1)) cloudtop = .true. if (qrl(k).lt. 0. .and. ri(k).gt.ricrit .and. cloudtop) then ch = (1 -sfuh(k))*chu(k) + sfuh(k)*chs(k) cm = (1 -sfuh(k))*cmu(k) + sfuh(k)*cms(k) n2h = ch*slslope(k) + cm*qtslope(k) if (n2h.le.0.) then !------------------------------------------------- ! Test if k and k+1 are part of the same preexist- ! ing CL. If not, find appropriate index for new ! SRCL. Note that this calculation makes use of ! ncv set from prior passes through the k do loop in_CL = .false. do while (ncv .le. ncvf) if (ktop(ncv) .le. k) then !--------------------------------------- ! If kbase > k, k and k+1 are part of ! same prior CL if (kbase(ncv) .gt. k) in_CL = .true. !--------------------------------------- ! exit from do-loop once CL top at/above ! intfc k. exit else !--------------------------------------- ! Go up one CL ncv = ncv + 1 end if end do ! ncv !------------------------------------------------- ! Add a new SRCL if (.not.in_CL) then ncvfin = ncvfin+1 ncvnew = ncvfin ktop(ncvnew) = k kbase(ncvnew) = k+1 belongcv(k) = .true. belongcv(k+1) = .true. if (k.lt.kdim) then ebrk(ncvnew) = 0. lbrk(ncvnew) = 0. shcl(ncvnew) = 0. smcl(ncvnew) = 0. else !--------------------------------------- ! surface radiatively driven fog if (bflxs.gt.0.) then !---------------------------------- ! unstable surface layer ! incorporate surface TKE into ebrk(ncvnew) = tkes lbrk(ncvnew) = zm(k) adj_time_inv(kdim+1) = sqrt(tkes)/zm(k) else !---------------------------------- ! stable surface layer ! don't incorporate surface TKE ebrk(ncvnew) = 0. lbrk(ncvnew) = 0. end if shcl(ncvnew) = 0. smcl(ncvnew) = 0. ncvsurf = ncvnew end if ! k < kdim end if ! new SRCL end if ! n2h < 0 end if ! qrl < 0, ri(k) > ricrit and cloudtop end do ! k do loop, end of SRCL section !----------------------------------------------------------------------- ! ! For each CL, compute length scale, r^2, e, Kh and Km if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' IN CALEDDY CALCULATION OF L, R2, e, Kh, a'//& 'nd Km' write (dpu,'(a)') ' ' end if do ncv = 1, ncvfin kt = ktop(ncv) kb = kbase(ncv) lbulk = zi(kt)-zi(kb) if (column_match) then write (dpu,'(a,i3)') ' convective layer #',ncv write (dpu,'(a,i4,i4)') ' kt, kb = ', kt,kb write (dpu,'(a,f14.7,a)') ' lbulk = ', lbulk, ' m' write (dpu,'(a)') ' ' end if do k = min(kb,kdim), kt, -1 leng(k) = lengthscale(zi(k),lbulk) wcap(k) = (leng(k)**2)*(-shcl(ncv)*n2(k)+ & smcl(ncv)*s2(k)) end do ! k do loop if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' k leng wcap ' write (dpu,'(a)') ' (m) (m2/s2)' write (dpu,'(a)') '----------------------------' write (dpu,'(a)') ' ' do k = min(kb,kdim),kt,-1 write(dpu,33) k,leng(k),wcap(k) enddo 33 format(1X,i2,1X,f9.4,1X,f14.9) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! Calculate jumps at the lower inversion if (kb .lt. kdim+1) then jbzm = zm(kb-1) - zm(kb) jbsl = sl(kb-1) - sl(kb) jbqt = qt(kb-1) - qt(kb) jbbu = n2(kb) * jbzm jbbu = max(jbbu,jbumin) jbu = u(kb-1) - u(kb) jbv = v(kb-1) - v(kb) ch = (1 -sflh(kb-1))*chu(kb) + sflh(kb-1)*chs(kb) cm = (1 -sflh(kb-1))*cmu(kb) + sflh(kb-1)*cms(kb) n2h = (ch*jbsl + cm*jbqt)/jbzm vyb = n2h*jbzm/jbbu vub = min(1.,(jbu**2+jbv**2)/(jbbu*jbzm) ) else !------------------------------------------------------ ! Zero bottom entrainment contribution for CL extending ! down to sfc vyb = 0. vub = 0. end if if (column_match .and. kb .lt. kdim+1) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' jumps at lower inversion ' write (dpu,'(a)') ' ------------------------ ' write (dpu,'(a)') ' ' write (dpu,'(a,i3)') ' inversion half level = ',kb write (dpu,'(a,f14.7,a)') ' jbzm = ',jbzm,' m' write (dpu,'(a,f14.7,a)') ' jbsl = ',jbsl/cp_air,' K' write (dpu,'(a,f14.7,a)') ' jbqt = ',jbqt*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' jbbu = ',jbbu,' m/s2' write (dpu,'(a,f14.7,a)') ' jbumin = ',jbumin,' m/s2' write (dpu,'(a,f14.7,a)') ' jbu = ',jbu,' m/s' write (dpu,'(a,f14.7,a)') ' jbv = ',jbv,' m/s' write (dpu,'(a,f14.7,a)') ' ch = ',ch*cp_air, ' m/K/s2' write (dpu,'(a,f14.7,a)') ' cm = ',cm, ' m/s2' write (dpu,'(a,f14.7,a)') ' sqrt(n2h) = ',n2h* & sqrt(abs(n2h))*3600./abs(n2h),' 1/sec' write (dpu,'(a,f14.7,a)') ' vyb = ', vyb write (dpu,'(a,f14.7,a)') ' vub = ', vub write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! Calculate jumps at the upper inversion ! Note the check to force jtbu to be greater than or equal ! to jbumin. jtzm = zm(kt-1) - zm(kt) jtsl = sl(kt-1) - sl(kt) jtqt = qt(kt-1) - qt(kt) jtbu = n2(kt) * jtzm jtbu = max(jtbu,jbumin) jtu = u(kt-1) - u(kt) jtv = v(kt-1) - v(kt) ch = (1 -sfuh(kt))*chu(kt) + sfuh(kt)*chs(kt) cm = (1 -sfuh(kt))*cmu(kt) + sfuh(kt)*cms(kt) n2h = (ch*jtsl + cm*jtqt)/jtzm ! Ratio of buoy flux to w'(b_l)' vys = n2h*jtzm/jtbu ! Inverse of shear prodution divided by buoyancy production vus = min(1.,(jtu**2+jtv**2)/(jtbu*jtzm)) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' jumps at upper inversion ' write (dpu,'(a)') ' ------------------------ ' write (dpu,'(a)') ' ' write (dpu,'(a,i3)') ' inversion half level = ',kt write (dpu,'(a,f14.7,a)') ' jtzm = ',jtzm,' m' write (dpu,'(a,f14.7,a)') ' jtsl = ',jtsl/cp_air,' K' write (dpu,'(a,f14.7,a)') ' jtqt = ',jtqt*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' jtbu = ',jtbu,' m/s2' write (dpu,'(a,f14.7,a)') ' jbumin = ',jbumin,' m/s2' write (dpu,'(a,f14.7,a)') ' jtu = ',jtu,' m/s' write (dpu,'(a,f14.7,a)') ' jtv = ',jtv,' m/s' write (dpu,'(a,f14.7,a)') ' ch = ',ch*cp_air, ' m/K/s2' write (dpu,'(a,f14.7,a)') ' cm = ',cm, ' m/s2' write (dpu,'(a,f14.7,a)') ' sqrt(n2h) = ',n2h* & sqrt(abs(n2h))*3600./abs(n2h),' 1/sec' write (dpu,'(a,f14.7,a)') ' vys = ', vys write (dpu,'(a,f14.7,a)') ' vus = ', vus write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! ! Calculate evaporative entrainment enhancement factor evhc. ! We take the full inversion strength to be jt2slv, where ! jt2slv = slv(kt-2) - slv(kt), and kt - 1 is in the ! ambiguous layer. However, for a cloud-topped CL overlain ! by another convective layer, it is possible that ! slv(kt-2) < slv(kt). To avoid negative or excessive evhc, ! we lower-bound jt2slv and upper-bound evhc. evhc = 1. cloudtop = .false. if (ql(kt) .gt. qcminfrac*qsl(kt) .and. use_qcmin .and. & ql(kt-1).lt. qcminfrac*qsl(kt-1)) cloudtop= .true. if (.not.use_qcmin .and. qa(kt).gt.qa(kt-1))cloudtop= .true. if (cloudtop) then if (use_qcmin) then tmpfrac = 1. else tmpfrac = 1. - (qa(kt-1)/qa(kt)) end if jt2slv = slv(max(kt-2,1)) - slv(kt) jt2slv = max(jt2slv, jbumin*slv(kt-1)/grav) if (use_extrapolated_ql) then evhc = 1.+tmpfrac*a2l*a3l*hleff(kt)*qltop(kt)/ & jt2slv else evhc = 1.+tmpfrac*a2l*a3l*hleff(kt)*ql(kt) / & jt2slv end if evhc = min(evhc,evhcmax) evhcvec(kt) = evhc end if if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' computing entrainment enhancem'//& 'ent factor ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' slv(max(kt-2,1))/cp_air = ',& slv(max(kt-2,1))/cp_air,' K' write (dpu,'(a,f14.7,a)') ' slv(kt)/cp_air = ', slv(kt)/& cp_air,' K' write (dpu,'(a,f14.7,a)') ' jt2slv(kt)/cp_air = ',jt2slv& /cp_air,' K' write (dpu,'(a,f14.7,a)') ' ql(kt-1) = ',ql(kt-1)* & 1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' ql(kt) = ',ql(kt)*1000., & ' g/kg' write (dpu,'(a,f14.7,a)') ' qltop(kt) = ',qltop(kt)* & 1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' qa(kt-1) = ',qa(kt-1) write (dpu,'(a,f14.7,a)') ' qa(kt) = ', qa(kt) write (dpu,'(a,f14.7,a)') ' evhc = ',evhc end if !----------------------------------------------------------- ! ! Radiative forcing at the upper inversion if at a cloud top if (cloudtop) then !------------------------------------------------------ ! estimate longwave opt depth in layer kt lwp = ql(kt) * (pi(kt+1) - pi(kt)) / grav opt_depth = 156*lwp !------------------------------------------------------ ! Approximation to LW cooling frac at inversion ! The following formula is a polynomial approximation ! to exact solution which is ! radinvfrac = 1 - 2/opt_depth + 2/(exp(opt_depth)-1)) radinvfrac = opt_depth*(4.+opt_depth) / & (6.*(4.+opt_depth) + opt_depth**2) !------------------------------------------------------ ! units of radf = (m*m)/(s*s*s) radf = max(-radinvfrac*qrl(kt)*(zi(kt)-zi(kt+1)),0.)* & cp_air * chs(kt) else lwp = 0. opt_depth = 0. radinvfrac = 0. radf = 0. end if !----------------------------------------------------------- ! ! Solve cubic equation ! r^3 + trmp*r + trmq = 0, r = sqrt ! to estimate for multilayer convection. Note, that if ! the CL goes to the surface, vyb and vub are zero, and ebrk ! and lbrk have already incorporated the surface interfacial ! layer, so the formulas below still apply. For a SRCL, ! there are no interior interfaces so ebrk = lbrk = 0. ! The cases are: ! (1) no cloudtop cooling (radf=0) -- trmq = 0, ! r = sqrt(-trmp) ! (2) radf > 0 but no interior CL interface -- trmp = 0, ! trmq < 0 ! (3) radf > 0 and interior CL interface(s) -- trmp < 0, ! trmq < 0 trma = 1. - (b1*a1l/lbulk)*(evhc*(-vys+vus)*(zi(kt)-zm(kt))& + (-vyb+vub)*(zm(kb-1)-zi(kb)) ) trma = max(trma,0.5) ! Prevents runaway entrainment instab. trmp = -ebrk(ncv) *(lbrk(ncv)/lbulk)/trma trmq = -b1*radf*(leng(kt)/lbulk)*(zi(kt)-zm(kt))/trma qq = (trmp/3.)**3+(trmq/2.)**2 if (trmq .lt. 0.) then if (qq .gt. 0.) then rootp = (-trmq/2.+sqrt(qq))**(1./3.) if (trmp .lt. 0.) then !-------------------------------------------- ! case 3 (in case 2, added term is zero) rootp = rootp + (-trmq/2.-sqrt(qq))**(1./3.) end if else !------------------------------------------------- ! also part of case 3 angle = acos(-trmq/2./sqrt(-(trmp/3)**3)) rootp = 2.*sqrt(-trmp/3.)*cos(angle/3.) end if else !------------------------------------------------------ ! case 1: radf = 0, so trmq = 0 rootp = sqrt(-trmp) endif if (column_match) then write (dpu,'(a,f14.7)') ' trma = ',trma write (dpu,'(a,f14.7,a)') ' trmp = ',trmp,' m2/s2' write (dpu,'(a,f14.7,a)') ' trmq = ',trmq,' m3/s3' write (dpu,'(a,f14.7,a)') ' qq = ', qq,' m6/s6' write (dpu,'(a,f14.7)') ' b1 = ', b1 write (dpu,'(a,f14.7)') ' a1l = ', a1l write (dpu,'(a,f14.7)') ' vys = ', vys write (dpu,'(a,f14.7)') ' vus = ', vus write (dpu,'(a,f14.7)') ' vyb = ', vyb write (dpu,'(a,f14.7)') ' vub = ', vub write (dpu,'(a,f14.7,a)') ' cwp = ', 1000.*lwp, & ' g/m2' write (dpu,'(a,f14.7)') ' opt_depth = ', opt_depth write (dpu,'(a,f14.7)') ' radinvfrac = ',radinvfrac write (dpu,'(a,f14.7,a)') ' radf = ', radf,' m2/s3' end if !----------------------------------------------------------- ! limit CL-avg TKE used for entrainment ebrk(ncv) = rootp**2 ebrk(ncv) = max(min(ebrk(ncv),tkemax),tkemin) wbrk(ncv) = ebrk(ncv)/b1 if (column_match) then write (dpu,'(a,f14.7,a)') ' rootp**2 = ', rootp**2, ' m2/s2' write (dpu,'(a,f14.7,a)') ' ebrk = ', ebrk(ncv), ' m2/s2' write (dpu,'(a,f14.7,a)') ' wbrk = ', wbrk(ncv), ' m2/s2' end if !----------------------------------------------------------- ! ebrk should be greater than zero so if it is not a FATAL ! call is implemented if (ebrk(ncv) .eq. 0.) then call error_mesg ('edt_mod', & 'convective layer average tke = 0', FATAL) end if !----------------------------------------------------------- ! Compute adjustment time for layer equal to lbulk / ! Should this be divided by "c" = b1/mu, which for ! default value = 5.8/70 = 0.083 do k = kb, kt, -1 adj_time_inv(k) = sqrt(ebrk(ncv))/lbulk enddo !----------------------------------------------------------- ! We approximate TKE = at entrainment interfaces con- ! sistent with entrainment closure. rcap(kt) = 1. rcap(kb) = 1. !----------------------------------------------------------- ! Calculate ratio rcap = e/ in convective layer interior. ! Bound it by limits rmin = 0.1 to rmax = 2.0 to take care ! of some pathological cases. if ((kb-kt).gt.1) then do k = kb-1, kt+1, -1 rcap(k) = (mu*leng(k)/lbulk + wcap(k)/wbrk(ncv)) / & (mu*leng(k)/lbulk + 1. ) rcap(k) = min(max(rcap(k),rmin), rmax) end do end if !----------------------------------------------------------- ! Compute TKE throughout CL, and bound by tkemin & tkemax. ! ! Question by Steve Klein: ! ! Does tke(kb) properly account if kb is the top interface ! of another convective layer? ! do k = kb, kt, -1 tke(k) = max(min(ebrk(ncv)*rcap(k),tkemax),tkemin) end do !----------------------------------------------------------- ! Compute CL interior diffusivities, buoyancy and shear ! production if ((kb-kt).gt.1) then do k = kb-1, kt+1, -1 kvh(k) = leng(k)*sqrt(tke(k))*shcl(ncv) kvm(k) = leng(k)*sqrt(tke(k))*smcl(ncv) bprod(k) = - kvh(k)*n2(k) sprod(k) = kvm(k)*s2(k) trans(k) = mu*(ebrk(ncv)-tke(k))* & adj_time_inv(k)/b1 diss(k) = sqrt(tke(k)*tke(k)*tke(k))/b1/ & leng(k) shvec(k) = shcl(ncv) smvec(k) = smcl(ncv) ghvec(k) = ghcl(ncv) turbtype(2,k) = 1 isturb(k) = 1. isturb(k-1) = 1. !------------------------------------------------- ! compute sdevs ! ! The approximation used is that qt and qs(Tl) are ! correlated with strength kappa. The sign of the ! correlation is positive if the vertical ! gradients to qt and sli are of the same sign, ! and negative otherwise. if (k .eq. kdim+1) call error_mesg ('edt_mod', & 'trying to compute sdevs at the surface', & FATAL) if (k .eq. 1) call error_mesg ('edt_mod', & 'trying to compute sdevs at the model top', & FATAL) qtsltmp = (qt(k-1) - qt(k))/(zm(k-1) - zm(k)) slsltmp = (sl(k-1) - sl(k))/(zm(k-1) - zm(k)) qsltmp = 0.5 * ( qsl (k-1) + qsl (k) ) dqsldtltmp = 0.5 * ( dqsldtl(k-1) + dqsldtl(k) ) hlefftmp = 0.5 * ( hleff (k-1) + hleff (k) ) sdevs(k) = ( (mesovar*qsltmp)**2.0) + & ( ( (kvh(k)/sqrt(tke(k))) * & (qtsltmp-(kappa*slsltmp*dqsldtltmp/cp_air)))**2.0) sdevs(k) = sqrt(sdevs(k)) / & (1.+hlefftmp*dqsldtltmp/cp_air) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a,i4)') ' sigmas for level ',k write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ', 1000.* & sdevs(k), ' g/kg' write (dpu,'(a,f14.7)') ' acoef = ', 1. / & ( 1. + hlefftmp*dqsldtltmp/cp_air ) write (dpu,'(a,f14.7,a)') ' sigmas/a = ', & 1000.*sdevs(k)*(1. +hlefftmp*dqsldtltmp/cp_air),& ' g/kg' write (dpu,'(a,f14.7)' ) ' mesovar = ', mesovar write (dpu,'(a,f14.7,a)') ' mesovar*qsl = ', & mesovar*qsltmp*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' turb.fluct = ',1000.& *(kvh(k)/sqrt(tke(k)))*(qtsltmp- & (kappa*slsltmp*dqsldtltmp/cp_air)) ,' g/kg' write (dpu,'(a,f14.7,a)') ' (kvh/sqrt(tke))* '//& ' qtsltmp = ',1000.*(kvh(k)/sqrt(tke(k)))* & qtsltmp ,' g/kg' write (dpu,'(a,f14.7,a)') ' (kvh/sqrt(tke))* '//& ' (kappa*slsltmp*dqsldtltmp/cp_air) = ',1000.* & (kvh(k)/sqrt(tke(k)))*(kappa*slsltmp* & dqsldtltmp/cp_air) ,' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if end do !------------------------------------------------------ ! ! set sdevs at kt and kb equal to kt+1 and kb-1 values ! respectively ! sdevs(kt) = sdevs(kt+1) sdevs(kb) = sdevs(kb-1) end if !----------------------------------------------------------- ! Compute diffusivity we*dz and some diagnostics at the ! upper inversion. Limit entrainment rate below the free ! entrainment limit a1l * sqrt(e) kentr = jtzm * a1l * sqrt(ebrk(ncv)) * & min(evhc * ebrk(ncv)/(jtbu*leng(kt)),1.) kvh(kt) = kentr kvm(kt) = kentr bprod(kt) = -kentr*n2h+radf sprod(kt) = kentr*s2(kt) trans(kt) = mu*(ebrk(ncv)-tke(kt))*adj_time_inv(kt)/b1 diss(kt) = sqrt(tke(kt)*tke(kt)*tke(kt))/b1/leng(kt) radfvec(kt) = radf turbtype(4,kt) = 1 isturb(kt) = 1. !----------------------------------------------------------- ! set isturb to 1 in the ambiguous layer isturb(kt-1) = 1. if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' at upper inversion: ' write (dpu,'(a,f14.7,a)') ' kentr = ', kentr, ' m2/s' write (dpu,'(a,f14.7)') ' mu = ', mu write (dpu,'(a,f14.7,a)') ' 1/adj_time_inv = ', & 1./adj_time_inv(kt), ' sec' write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! compute sdevs at CL top ! ! FOR SOME REASON THIS DERIVATION DOES NOT WORK ! ! note that a water perturbation q* and buoyancy ! perturbation b* are computed from the below equations: ! ! q* = (w'qt') / sqrt(tke) ! = - kentr * (dqt/dz) / sqrt(tke) ! ! b* = (w'b') / sqrt(tke) ! = bprod / sqrt(tke) ! ! qtsltmp = (qt(kt-1) - qt(kt))/(zm(kt-1) - zm(kt)) qstartmp = - kentr * qtsltmp / sqrt(tke(kt)) bstartmp = bprod(kt) / sqrt(tke(kt)) temperature = pm(kt)/density(kt)/rdgas !sdevs(kt) = ((mesovar*qsl(kt))**2.0) + & ! ((qstartmp-dqsldtl(kt)*bstartmp*temperature/grav)**2.0) !sdevs(kt) = sqrt(sdevs(kt))/(1.+hleff(kt)*dqsldtl(kt)/cp_air) if ((kb-kt).le.1) sdevs(kt) = mesovar*qsl(kt)/ & (1.+hleff(kt)*dqsldtl(kt)/cp_air) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a,i4)') ' sigmas at kt level # ',kt write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ', 1000.*sdevs(kt), & ' g/kg' write (dpu,'(a,f14.7)') ' acoef = ', 1. /(1.+hleff(kt)* & dqsldtl(kt)/cp_air ) write (dpu,'(a,f14.7,a)') ' sigmas/a = ', 1000.*sdevs(kt)*& ( 1. + hleff(kt)*dqsldtl(kt)/cp_air ), ' g/kg' write (dpu,'(a,f14.7)' ) ' mesovar = ', mesovar write (dpu,'(a,f14.7,a)') ' mesovar*qsl = ',mesovar* & qsl(kt)*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' turb.fluct = ',1000.* & (qstartmp-dqsldtl(kt)*bstartmp*temperature/grav) , & ' g/kg' write (dpu,'(a,f14.7,a)') ' temperature = ',temperature, ' K' write (dpu,'(a,f14.7,a)') ' qstartmp = ',1000.*qstartmp , ' g/kg' write (dpu,'(a,f14.7,a)') ' bstartmp = ',bstartmp ,' m/s2' write (dpu,'(a,f14.7,a)') ' jtbu = ',jtbu ,' m/s2' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------- ! Compute Kh, Km and some diagnostics at the lower inversion if (kb .lt. kdim+1) then kentr = jbzm * a1l * sqrt(ebrk(ncv)) * & min(ebrk(ncv)/(jbbu*leng(kb)),1.) !------------------------------------------------------ ! The addition of the previous value of kvh and khm ! handles the case of 2 CLs entraining into each other kvh(kb) = kentr + kvh(kb) kvm(kb) = kentr + kvm(kb) bprod(kb) = bprod(kb) - kvh(kb)*n2(kb) sprod(kb) = sprod(kb) + kvm(kb)*s2(kb) trans(kb) = trans(kb) + mu*(ebrk(ncv)-tke(kb))* & adj_time_inv(kb)/b1 diss(kb) = diss(kb)+sqrt(tke(kb)*tke(kb)*tke(kb)) & /b1/leng(kb) turbtype(3,kb) = 1 isturb(kb-1) = 1. !------------------------------------------------------ ! set isturb to 1 in the ambiguous layer isturb(kb) = 1. if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' at lower inversion: ' write (dpu,'(a,f14.7,a)') ' kentr = ', kentr, ' m2/s' write (dpu,'(a,f14.7)') ' mu = ', mu write (dpu,'(a,f14.7,a)') ' adj_time_inv = ', & 1./adj_time_inv(kb), ' sec' write (dpu,'(a)') ' ' end if !------------------------------------------------------ ! compute sdevs at CL bottom ! ! note that same formulas are used here as for CL top qtsltmp = (qt(kb-1) - qt(kb))/(zm(kb-1) - zm(kb)) qstartmp = - kentr * qtsltmp / sqrt(tke(kb)) bstartmp = bprod(kb) / sqrt(tke(kb)) temperature = pm(kb-1)/density(kb-1)/rdgas !sdevs(kb) = ((mesovar*qsl(kb-1))**2.0) + ((qstartmp-& ! dqsldtl(kb-1)*bstartmp*temperature/grav)& ! **2.0) !sdevs(kb) = sqrt(sdevs(kb))/(1.+hleff(kb-1)* & ! dqsldtl(kb-1)/cp_air) if ((kb-kt).le.1) sdevs(kb) = mesovar*qsl(kb-1)/ & (1.+hleff(kb-1)*dqsldtl(kb-1)/cp_air) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a,i5)') ' sigmas at lower inversion l'//& 'evel # ',kb write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ', 1000.* & sdevs(kb), ' g/kg' write (dpu,'(a,f14.7) ') ' acoef = ', 1. /(1.+ & hleff(kb-1)*dqsldtl(kb-1)/cp_air ) write (dpu,'(a,f14.7,a)') ' sigmas/a = ', 1000.* & sdevs(kb)*( 1. + hleff(kb-1)*dqsldtl(kb-1)/cp_air ), & ' g/kg' write (dpu,'(a,f14.7)' ) ' mesovar = ', mesovar write (dpu,'(a,f14.7,a)') ' mesovar*qsl = ',mesovar* & qsl(kb-1)*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' turb.fluct = ',1000.* & (qstartmp-dqsldtl(kb-1)*bstartmp*temperature/ & grav) , ' g/kg' write (dpu,'(a,f14.7,a)') ' temperature = ', & temperature,' K' write (dpu,'(a,f14.7,a)') ' qstartmp = ',1000.* & qstartmp ,' g/kg' write (dpu,'(a,f14.7,a)') ' bstartmp = ',bstartmp , ' m/s2' write (dpu,'(a,f14.7,a)') ' jbbu = ',jbbu ,' m/s2' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if end if !----------------------------------------------------------- ! put minimum threshold on TKE to prevent possible division ! by zero. if (kb .lt. kdim+1) then wcap(kb) = bprod(kb)*leng(kb)/sqrt(max(tke(kb),tkemin)) else wcap(kb) = tkes / b1 end if wcap(kt) = bprod(kt)*leng(kt) / sqrt(max(tke(kt),tkemin)) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') 'k leng rcap wcap tk'//& 'e' write (dpu,'(a)') ' (m) (m2/s2) (m2'//& '/s2)' write (dpu,'(a)') '------------------------------------'//& '----' write (dpu,'(a)') ' ' do k = kb,kt,-1 write(dpu,34) k,leng(k),rcap(k),wcap(k),tke(k) enddo 34 format(1X,i2,1X,4(f8.4,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k kvh kvm bprod '//& ' sprod trans diss' write (dpu,'(a)') ' (m2/s) (m2/s) (m2/s3) '//& ' (m2/s3) (m2/s3) (m2/s3)' write (dpu,'(a)') '------------------------------------'//& '-------------------------------------' write (dpu,'(a)') ' ' do k = kb,kt,-1 write(dpu,35) k,kvh(k),kvm(k),bprod(k),sprod(k),trans(k), & diss(k) enddo 35 format(1X,i2,1X,2(f8.4,1X),4(f12.9,1X)) write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !---------------------------------------------------------------- ! End big loop over ncv end do !----------------------------------------------------------------------- ! ! If the lowest CL reaches the surface, define the PBL ! depth as the CL top. ! if (ncvsurf .gt. 0) then pblh = zi(ktop(ncvsurf)) if(bflxs.ge.0.) then turbtype(2,kdim+1) = 1 else turbtype(3,kdim+1) = 1 end if else pblh = 0. end if !----------------------------------------------------------------------- ! ! STABLE TURBULENT LAYERS ! !---------------------------------------------------------------- ! find turbulent lengthscales in all stable turbulent layers belongst(1) = .false. ! k = 1 assumed nonturbulent any_stable = .false. do k = 2, kdim belongst(k) = (ri(k) .lt. ricrit) .and. (.not. belongcv(k)) if (belongst(k)) any_stable = .true. if (belongst(k) .and. (.not.belongst(k-1)) ) then kt = k ! Top of stable turb layer else if (.not. belongst(k) .and. belongst(k-1)) then kb = k-1 ! Base of stable turb layer lbulk = zm(kt-1) - zm(kb) do ks = kt, kb leng(ks)=lengthscale(zi(ks),lbulk) adj_time_inv(ks) = 1. / lbulk end do end if end do ! k !---------------------------------------------------------------- ! Now look whether stable turb layer extends to ground. Note that ! interface kdim+1 is assumed to always be stable-turbulent if it ! is not convective. Note that if it is convective, kdim will ! also be convective, so the above loop will have finished ! finding all elevated turbulent layers. belongst(kdim+1) = .not. belongcv(kdim+1) if (belongst(kdim+1)) then turbtype(1,kdim+1) = 1 if (belongst(kdim)) then !------------------------------------------------------ ! surface stable layer includes interior stable turb- ! ulent interface kt already defined above. Note that ! zm(kb) = 0. lbulk = zm(kt-1) else !------------------------------------------------------ ! surface stable BL with no interior turbulence kt = kdim+1 end if lbulk = zm(kt-1) pblh = lbulk ! PBL Height <-> lowest stable turb. layer do ks = kt,kdim+1 leng(ks) = lengthscale(zi(ks),lbulk) adj_time_inv(ks) = 1./ lbulk end do adj_time_inv(kdim+1) = adj_time_inv(kdim+1)*sqrt(tkes) end if ! for belongst if !---------------------------------------------------------------- ! Calculate tke, kvh, kvm if (column_match .and. any_stable) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' STABLE LAYERS ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' k kvh kvm leng t'//& 'ke bprod sprod diss' write (dpu,'(a)') ' (m2/s) (m2/s) (m) (m2'//& '/s2) (m2/s3) (m2/s3) (m2/s3) ' write (dpu,'(a)') '------------------------------------'//& '-----------------------------------------' write (dpu,'(a)') ' ' end if do k = 2, kdim if (belongst(k)) then turbtype(1,k) = 1 isturb(k) = 1. isturb(k-1) = 1. call galperin(ri(k),gh,ckh,ckm) ! note that original cb code did not limit the maximum ! gh to 0.0233 so that we should check that if the ! limiter was invoked !!RSH drop this fail for now (maybe ok during early spinup ?? ) ! if (gh .eq. 0.0233) call error_mesg ('edt_mod',& ! 'galperin stability fn = 0.0233',& ! FATAL) ! ri(k) should be less than ricrit if (ri(k) .gt. ricrit) call error_mesg ('edt_mod',& 'ri(k) > ricrit, but belongst(k) = T',FATAL) tke(k) = b1*(leng(k)**2)*(-ckh*n2(k)+ckm*s2(k)) tke(k) = max(min(tke(k),tkemax),tkemin) kvh(k) = leng(k) * sqrt(tke(k)) * ckh kvm(k) = leng(k) * sqrt(tke(k)) * ckm bprod(k) = - kvh(k)*n2(k) sprod(k) = kvm(k)*s2(k) diss (k) = sqrt(tke(k)*tke(k)*tke(k))/b1/leng(k) adj_time_inv(k) = adj_time_inv(k)*sqrt(tke(k)) shvec(k) = ckh smvec(k) = ckm ghvec(k) = gh if (column_match) then write(dpu,37) k,kvh(k),kvm(k),leng(k),tke(k),bprod(k),& sprod(k),diss(k) 37 format(1X,i2,1X,3(f8.4,1X),4(f12.9,1X)) end if !------------------------------------------------------ ! compute sdevs ! (note same formulas as used in CL calculations) ! if (k .eq. kdim+1) call error_mesg ('edt_mod', & 'trying to compute sdevs at the surface', FATAL) if (k .eq. 1) call error_mesg ('edt_mod', & 'trying to compute sdevs at the model top', FATAL) qtsltmp = (qt(k-1) - qt(k))/(zm(k-1) - zm(k)) slsltmp = (sl(k-1) - sl(k))/(zm(k-1) - zm(k)) qsltmp = 0.5 * ( qsl (k-1) + qsl (k) ) dqsldtltmp = 0.5 * ( dqsldtl(k-1) + dqsldtl(k) ) hlefftmp = 0.5 * ( hleff (k-1) + hleff (k) ) sdevs(k) = ( (mesovar*qsltmp)**2.0) + & ( ( (kvh(k)/sqrt(tke(k))) * & (qtsltmp-(kappa*slsltmp*dqsldtltmp/cp_air)) )& **2.0) sdevs(k) = sqrt(sdevs(k)) / (1.+hlefftmp*dqsldtltmp/cp_air) if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' write (dpu,'(a)') ' sigmas .... ' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ',1000.*sdevs(k), ' g/kg' write (dpu,'(a,f14.7) ') ' acoef = ', 1. / ( 1. + & hlefftmp*dqsldtltmp/cp_air ) write (dpu,'(a,f14.7,a)') ' sigmas/a = ', 1000.* & sdevs(k)* ( 1. + hlefftmp*dqsldtltmp/cp_air ), ' g/kg' write (dpu,'(a,f14.7)' ) ' mesovar = ', mesovar write (dpu,'(a,f14.7,a)') ' mesovar*qsl = ',mesovar* & qsltmp*1000.,' g/kg' write (dpu,'(a,f14.7,a)') ' turb.fluct = ',1000.* & (kvh(k)/sqrt(tke(k)))*(qtsltmp-(kappa*slsltmp* & dqsldtltmp/cp_air)) ,' g/kg' write (dpu,'(a,f14.7,a)') ' (kvh/sqrt(tke))*qtsltm'//& 'p = ',1000.*(kvh(k)/sqrt(tke(k)))*qtsltmp ,' g/kg' write (dpu,'(a,f14.7,a)') ' (kvh/sqrt(tke))*(kappa'//& '*slsltmp*dqsldtltmp/cp_air) = ',1000.*(kvh(k)/sqrt & (tke(k)))*(kappa*slsltmp*dqsldtltmp/cp_air),' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if end if ! belongs to stable layer end do ! k loop !----------------------------------------------------------------------- ! ! set tke at surface equal to diagnostic variable tkes tke(kdim+1) = tkes if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' pblh = ', pblh, ' m' write (dpu,'(a,f14.7,a)') ' tkes = ', tkes, ' m2/s2' write (dpu,'(a)') ' ' do k = kdim-n_print_levels,kdim write (dpu,'(a,i4,a,f14.7)') 'k = ',k,'; isturb =',isturb(k) enddo write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine caleddy ! !======================================================================= !======================================================================= ! ! subroutine galperin ! subroutine galperin(ricl,gh,sh,sm) ! !----------------------------------------------------------------------- ! ! Given a Richardson number ricl, return the stability functions ! sh and sm calculated according Galperin (1982). ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! real, intent(in) :: ricl real, intent(out) :: gh, sh, sm ! internal variables real :: ri, trma, trmb, trmc, det !----------------------------------------------------------------------- ! ! code ! ri = min(ricl,0.163) trma = alph3*alph4*ri+2.*b1*(alph2-alph4*alph5*ri) trmb = ri*(alph3+alph4)+2.*b1*(-alph5*ri+alph1) trmc = ri det = max(trmb*trmb-4.*trma*trmc,0.) gh = (-trmb + sqrt(det))/2./trma gh = max(gh,-0.28) gh = min(gh,0.0233) sh = alph5 / (1.+alph3*gh) sm = (alph1 + alph2*gh)/(1.+alph3*gh)/(1.+alph4*gh) !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine galperin ! !======================================================================= !======================================================================= ! ! function lengthscale ! function lengthscale(height,depth) ! !----------------------------------------------------------------------- ! ! Calculate the turbulent length scale given the depth of the ! turbulent layer. Near the surface, the lengthscale asymptotes ! to vonkarm*height. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! real :: lengthscale real, intent(in) :: height,depth !----------------------------------------------------------------------- ! ! code ! lengthscale = vonkarm*height/(1.+(vonkarm*height/(tunl*depth))) !----------------------------------------------------------------------- ! ! function end ! end function lengthscale ! !======================================================================= !======================================================================= ! ! subroutine gaussian_cloud ! ! ! this subroutine computes a new cloud fraction and cloud conden- ! sate using the Gaussian cloud modell of Mellor (1977) and ! Sommeria and Deardorff (1977). ! ! In this derivation a background mesoscale variability to total ! water variability equivalent to a fraction, mesovar, of qsl is ! assumed. ! subroutine gaussian_cloud (qxtop, qxmid, qxbot, acoef, sigmasf, qalyr, & qclyr, sfuh, sflh) !----------------------------------------------------------------------- ! ! variables ! ! ----- ! input ! ----- ! ! qxtop saturation excess at the top of the layer ! (kg wat/kg air) ! qxmid saturation excess at the midpoint of the layer ! (kg wat/kg air) ! qxbot saturation excess at the bottom of the layer ! (kg wat/kg air) ! sigmasf standardard devations of water perturbation s ! (kg water / kg air) ! acoef thermo coefficient = 1./(1.+ L*dqsdT/cp) ! ! ------ ! output ! ------ ! ! qalyr layer mean cloud fraction (fraction) ! qclyr layer mean cloud condensate (kg condensate/kg air) ! sfuh saturated fraction of the upper half of the layer ! sflh saturated fraction of the lower half of the layer ! ! -------- ! internal ! -------- ! ! !----------------------------------------------------------------------- real, intent(in) :: qxtop, qxmid, qxbot, acoef, sigmasf real, intent(out):: qalyr, qclyr, sfuh, sflh real :: q1top,q1mid,q1bot,cftop,cfmid,cfbot,qctop,qcmid,qcbot real :: qcuh,qclh !----------------------------------------------------------------------- ! ! initialize variables qalyr = 0.0 qclyr = 0.0 sfuh = 0.0 sflh = 0.0 q1top = 0.0 q1mid = 0.0 q1bot = 0.0 cftop = 0.0 cfmid = 0.0 cfbot = 0.0 qctop = 0.0 qcmid = 0.0 qcbot = 0.0 qcuh = 0.0 qclh = 0.0 !----------------------------------------------------------------------- ! ! compute cloud fraction and cloud condensate at layer top, ! midpoint and bottom q1top = acoef * qxtop / sigmasf cftop = max ( 0.5 * ( 2. - erfcc( q1top / sqrt(2.) ) ) , 0.0 ) qctop = cftop * q1top + ( exp(-0.5*q1top*q1top) / sqrt(2*fpi) ) qctop = sigmasf * max ( qctop, 0. ) q1mid = acoef * qxmid / sigmasf cfmid = max ( 0.5 * ( 2. - erfcc( q1mid / sqrt(2.) ) ) , 0.0 ) qcmid = cfmid * q1mid + ( exp(-0.5*q1mid*q1mid) / sqrt(2*fpi) ) qcmid = sigmasf * max ( qcmid, 0. ) q1bot = acoef * qxbot / sigmasf cfbot = max ( 0.5 * ( 2. - erfcc( q1bot / sqrt(2.) ) ) , 0.0 ) qcbot = cfbot * q1bot + ( exp(-0.5*q1bot*q1bot) / sqrt(2*fpi) ) qcbot = sigmasf * max ( qcbot, 0. ) sfuh = 0.5 * (cftop + cfmid) sflh = 0.5 * (cfmid + cfbot) qcuh = 0.5 * (qctop + qcmid) qclh = 0.5 * (qcmid + qcbot) qalyr = 0.5 * ( sfuh + sflh ) qclyr = 0.5 * ( qcuh + qclh ) if (qalyr .lt. qcminfrac) then qalyr = 0. qclyr = 0. end if if (column_match) then write (dpu,'(a)') ' ' write (dpu,'(a)') ' ========================================'//& '======= ' write (dpu,'(a)') ' GAUSSIAN CLOUD MODEL ' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' sigmas = ', 1000.*sigmasf, ' g/kg' write (dpu,'(a,f14.7)' ) ' acoef = ', acoef write (dpu,'(a,f14.7,a)') ' qxtop = ', qxtop*1000.,' g/kg' write (dpu,'(a,f14.7)') ' q1top = ', q1top write (dpu,'(a,f14.7)') ' cftop = ', cftop write (dpu,'(a,f14.7,a)') ' qctop = ', 1000.*qctop,' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' qxmid = ', qxmid*1000.,' g/kg' write (dpu,'(a,f14.7)') ' q1mid = ', q1mid write (dpu,'(a,f14.7)') ' cfmid = ', cfmid write (dpu,'(a,f14.7,a)') ' qcmid = ', 1000.*qcmid,' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7,a)') ' qxbot = ', qxbot*1000.,' g/kg' write (dpu,'(a,f14.7)') ' q1bot = ', q1bot write (dpu,'(a,f14.7)') ' cfbot = ', cfbot write (dpu,'(a,f14.7,a)') ' qcbot = ', 1000.*qcbot,' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7)') ' sfuh = ', sfuh write (dpu,'(a,f14.7)') ' sflh = ', sflh write (dpu,'(a,f14.7,a)') ' qcuh = ', 1000.*qcuh, ' g/kg' write (dpu,'(a,f14.7,a)') ' qclh = ', 1000.*qclh, ' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a,f14.7)') ' qalyr = ', qalyr write (dpu,'(a,f14.7,a)') ' qclyr = ', 1000.*qclyr, ' g/kg' write (dpu,'(a)') ' ' write (dpu,'(a)') ' ' end if !----------------------------------------------------------------------- ! ! subroutine end ! end subroutine gaussian_cloud ! !======================================================================= !======================================================================= ! ! function erfcc ! function erfcc(x) ! !----------------------------------------------------------------------- ! ! This numerical recipes routine calculates the complementary ! error function. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! real :: erfcc real, intent(in) :: x real :: t,z !----------------------------------------------------------------------- ! ! code ! z=abs(x) t=1./(1.+0.5*z) erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+t* & (.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+t* & (1.48851587+t*(-.82215223+t*.17087277))))))))) if (x.lt.0.) erfcc=2.-erfcc !----------------------------------------------------------------------- ! ! function end ! end function erfcc ! !======================================================================= end module edt_mod