! !IBM* SOURCEFORM(FREE(F90)) ! module gw_drag ! !--------------------------------------------------------------------------------- ! Purpose: ! ! Module to compute the forcing due to parameterized gravity waves. Both an ! orographic and an internal source spectrum are considered. ! ! Author: Byron Boville ! ! B.Foster, Oct, 2003: Adapted for timegcm1 ! !--------------------------------------------------------------------------------- ! ! gw_params module (gw_params.F) is shared by timegcm mgw module and this module. use params_module,only: nlev,dlon,nlonp4 use gw_share,only: & ! timegcm gw_params.F pcnst,pnats, & ! number of constituents pgwv, & ! number of waves allowed cpair, & ! specific heat of dry air pr_num, & ! prandtl number as in waccm/vertical_diffusion.F90 km_fac, & ! molecular viscosity constant orog,sghdat ! Orography data (see orogdat.F) use addfld_module,only: addfld implicit none save private ! Make default type private to the module ! ! PUBLIC: interfaces ! public gw_inti ! Initialization public gw_intr ! interface to actual parameterization ! ! PRIVATE: Rest of the data and interfaces are private to this module ! integer,parameter :: r8 = selected_real_kind(12) ! 8 byte real !+waccm+ ! integer, parameter :: pgwv = 4 ! number of waves allowed !-waccm- integer :: kbotbg, kbotoro ! interface of gwd source integer :: ktopbg, ktoporo ! top interface of gwd region real(r8) :: alpha(0:nlev) ! newtonian cooling coefficients real(r8) :: c(-pgwv:pgwv) ! list of wave phase speeds real(r8) :: dback ! background diffusivity real(r8) :: effkwv ! effective wavenumber (fcrit2*kwv) real(r8) :: effgw_oro = 1.0_r8 ! tendency efficiency real(r8) :: effgw_spec=.125_r8 ! tendency efficiency real(r8) :: fracldv ! fraction of stress deposited in low level region real(r8) :: g ! acceleration of gravity real(r8) :: kwv ! effective horizontal wave number real(r8) :: lzmax ! maximum vertical wavelength at source real(r8) :: mxasym ! max asymmetry between tau(c) and tau(-c) real(r8) :: mxrange ! max range of tau for all c real(r8) :: n2min ! min value of bouyancy frequency real(r8) :: fcrit2 ! critical froude number real(r8) :: oroko2 ! 1/2 * horizontal wavenumber real(r8) :: oroe ! efficiency factor for orography real(r8) :: orohmin ! min surface displacment height for orographic waves real(r8) :: orovmin ! min wind speed for orographic waves real(r8) :: r ! gas constant for dry air real(r8) :: rog ! r / g real(r8) :: taubgnd ! background source strength (/tauscal) real(r8) :: taumin ! minimum (nonzero) stress real(r8) :: tauscal ! scale factor for background stress source real(r8) :: tndmax ! maximum wind tendency real(r8) :: umcfac ! factor to limit tendency to prevent reversing u-c real(r8) :: ubmc2mn ! min (u-c)**2 real(r8) :: zldvcon ! constant for determining zldv from tau0 !+waccm+ integer, parameter :: nalph=66 ! no. of levels of pre-calculated newtonian cooling integer, parameter :: lchnk=0 ! nominal defination of lchnk real(r8) :: alpha0(nalph) ! newtonian cooling (1/day) real(r8) :: palph(nalph) ! pressure levs of pre-calculated newtonian cooling (hpa) data alpha0 / 1.896007 , 1.196965 , 0.7251356 , 0.6397463 , & 0.5777858 , 0.5712274 , 0.6836302 , 0.6678557 , & 0.5683219 , 0.4754283 , 0.3960519 , 0.332022 , & 0.2497581 , 0.168667 , 0.1323903 , 0.1257139 , & 0.1069889 , 0.09873954 , 0.09215571 , 0.09398635 , & 0.1061087 , 0.1294598 , 0.1544743 , 0.1648226 , & 0.1687332 , 0.1691513 , 0.1664987 , 0.159048 , & 0.149292 , 0.1351563 , 0.1174998 , 0.09913579 , & 0.08300615 , 0.0707 , 0.0615588 , 0.0542623 , & 0.0478562 , 0.04132157 , 0.03454087 , 0.02296682 , & 0.006723819, 0.02164464 , 0.05756261 , 0.003844868 , & 0.02929285 , 0.006627098, 0.04558291 , 0.02042176 , & 0.00000000 , 0.005880283, 0.00689498 , 0.01343466 , & 0.00000000 , 0.03415992 , 0.02855049 , 0.01688839 , & 0.0272628 , 0.02772121 , 0.02135626 , 0.04863235 , & 0.04568304 , 0.00000000 , 0.009604108, 0.00000000 , & 0.00000000 , 0.00000000 / data palph / 5.11075e-06 , 9.8269e-06 , 1.620185e-05 , 2.671225e-05 , & 4.4041e-05 , 7.261275e-05 , 0.000119719 , 0.00019738 , & 0.0003254225, 0.0005365325 , 0.0008846025 , 0.001458458 , & 0.002404575 , 0.00397825 , 0.006556825 , 0.01081382 , & 0.017898 , 0.02955775 , 0.04873075 , 0.07991075 , & 0.1282732 , 0.19812 , 0.292025 , 0.4101675 , & 0.55347 , 0.73048 , 0.9559475 , 1.244795 , & 1.61285 , 2.079325 , 2.667425 , 3.404875 , & 4.324575 , 5.4654 , 6.87285 , 8.599725 , & 10.70705 , 13.26475 , 16.35175 , 20.05675 , & 24.479 , 29.728 , 35.92325 , 43.19375 , & 51.6775 , 61.5205 , 72.8745 , 85.65715 , & 100.5147 , 118.2503 , 139.1154 , 163.6621 , & 192.5399 , 226.5132 , 266.4812 , 313.5013 , & 368.818 , 433.8952 , 510.4553 , 600.5242 , & 696.7963 , 787.7021 , 867.1607 , 929.6489 , & 970.5548 , 992.5561 / !-waccm- contains !=============================================================================== subroutine gw_inti (cpairx, cpwv, gx, rx, hypi, pver) !----------------------------------------------------------------------- ! Time independent initialization for multiple gravity wave parameterization. !----------------------------------------------------------------------- ! use history, only: addfld, add_default, phys_decomp !------------------------------Arguments-------------------------------- integer,intent(in) :: pver real(r8), intent(in) :: cpairx ! specific heat of dry air (constant p) real(r8), intent(in) :: cpwv ! specific heat of water vapor (constant p) real(r8), intent(in) :: gx ! acceleration of gravity real(r8), intent(in) :: rx ! gas constant for dry air real(r8), intent(in) :: hypi(pver+1) ! reference interface pressures !---------------------------Local storage------------------------------- integer :: k !+waccm+ integer :: m character(len=80) fnum,fname ! dummies !-waccm- !----------------------------------------------------------------------- ! Copy model constants ! cpair = cpairx g = gx r = rx ! Set MGWD constants kwv = 6.28e-5 ! 100 km wave length dback = 0.05 ! background diffusivity fcrit2 = 0.5 ! critical froude number squared oroe = 0.125 ! orographic efficiency factor tauscal= 0.001 ! scale factor for background stress !+waccm+ ! ! Tunable for time-gcm: ! taubgnd= 5.6 ! background stress amplitude ! taubgnd= 0.5 ! background stress amplitude taubgnd= 0.4 ! background stress amplitude !-waccm- fracldv= 0.7 ! fraction of tau0 diverged in low level region zldvcon= 10. ! constant for determining zldv lzmax = 7.E3 ! maximum vertical wavelength at source (m) ! Set phase speeds do k = -pgwv, pgwv ! c(k) = 10. * k ! 0, +/- 10, +/- 20, ... m/s c(k) = 15. * k ! 0, +/- 15, +/- 30, ... m/s end do ! if (masterproc) then ! write(6,*) ' ' ! write(6,*) 'GW_INTI: pgwv = ', pgwv ! write(6,*) 'GW_INTI: c(l) = ', c ! write(6,*) ' ' ! end if write(6,"('gw_inti: pgwv = ',i3)") pgwv ! !+waccm+ ! pre-calculated newtonian damping: ! * convert to 1/s ! * ensure it is not smaller than 1e-6 ! * convert palph from hpa to pa do k=1,nalph alpha0(k) = alpha0(k) / 86400. if (alpha0(k) .lt. 1.e-6) alpha0(k) = 1.e-6 palph(k) = palph(k) *1.e2 enddo ! interpolate to current vertical grid and obtain alpha call lininterp (alpha0 ,palph, 1 , nalph , alpha , hypi , pver+1) ! if (masterproc) then ! subscript out of bounds error here if -C is set ! write (6,*) 'gw_inti: newtonian damping (1/day):' ! write (6,fmt='(a4,a12,a10)') ' k ',' hypi ',' alpha ' ! do k=0,pver ! write (6,fmt='(i4,1e12.5,1f10.2)') k,hypi(k+1),alpha(k)*86400. ! enddo ! endif !-waccm- ! set radiative damping times !!$ do k = 0, pver !!$ alpha(k) = 1.e-6 ! about 10 days. !!$ end do ! Min and max values to keep things reasonable ! mxasym = 0.1 ! max factor of 10 from |tau(c)| to |tau(-c)| ! mxrange= 0.001 ! factor of 100 from max to min |tau(c)| mxasym = 1.e-12 ! max factor of 10 from |tau(c)| to |tau(-c)| mxrange= 1.e-12 ! factor of 100 from max to min |tau(c)| n2min = 1.e-8 ! min value of Brunt-Vaisalla freq squared orohmin= 10. ! min surface displacement for orographic wave drag orovmin= 2. ! min wind speed for orographic wave drag taumin = 1.e-13 ! min stress considered > 0 tndmax = 500. / 86400. ! max permitted tendency (500 m/s/day) umcfac = 0.5 ! max permitted reduction in u-c ubmc2mn= 0.01 ! min value of (u-c)^2 ! Determine other derived constants oroko2 = 0.5 * kwv effkwv = oroe * fcrit2 * kwv rog = r/g ! Determine the bounds of the background and orographic stress regions ! (cam vertical order, i.e., k=0 at top, k=pver at bottom) ktopbg = 0 kbotbg = pver ktoporo = 0 kbotoro = pver ! do k = 0, pver ! if (hypi(k+1) .lt. 10000.) kbotbg = k ! spectrum source at 100 mb ! if (hypi(k+1) .lt. 3000.) ktoporo = k ! end do !+waccm+ !!$ ktoporo = 0 !-waccm- ! if (masterproc) then ! write (6,*) 'KTOPBG =',ktopbg ! write (6,*) 'KBOTBG =',kbotbg ! write (6,*) 'KTOPORO =',ktoporo ! write (6,*) 'KBOTORO =',kbotoro ! end if return end subroutine gw_inti !=============================================================================== subroutine gw_intr ( & u,v,t,q, & s,zm,sgh,pblh,dt,landfrac,pmid,pint,& lnpint,pdel,rpdel,rlat, & utend,vtend,ttend,difk,qtend, & pcols,pver,lon0,lon1,lat) integer,intent(in) :: pcols,pver,lon0,lon1,lat real(r8),intent(in) :: rlat ! current latitude in radians real(r8),dimension(pcols,pver),intent(in) :: & u,v,t,s,zm,pmid,pdel,rpdel real(r8),dimension(pcols,0:pver),intent(in) :: pint,lnpint real(r8),dimension(pcols,pver,pcnst+pnats),intent(in) :: q real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography real(r8), intent(in) :: pblh(pcols) ! planetary boundary layer height real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: landfrac(pcols) ! Land fraction real(r8),dimension(pcols,pver),intent(out) :: & utend,vtend,ttend ! tendencies of zonal wind, meridional wind, and ! temperature (units: m/s/s, m/s/s, J/kg/s) real(r8),dimension(pcols,0:pver),intent(out) :: & difk ! eddy viscosity (interface levels, unit: m^2/s). real(r8),dimension(pcols,pver,pcnst+pnats),intent(out) :: qtend !---------------------------Local storage------------------------------- integer :: lchnk ! chunk identifier integer :: ncol ! number of atmospheric columns integer :: i,k ! loop indexes integer :: kldv(pcols) ! top interface of low level stress divergence region integer :: kldvmn ! min value of kldv integer :: ksrc(pcols) ! index of top interface of source region integer :: ksrcmn ! min value of ksrc real(r8) :: ttgw(pcols,pver) ! temperature tendency (J/kg/s) real(r8) :: utgw(pcols,pver) ! zonal wind tendency real(r8) :: vtgw(pcols,pver) ! meridional wind tendency real(r8) :: difkgw(pcols,0:pver) ! eddy diffusion due to wave breaking real(r8) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency real(r8) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency real(r8) :: rdpldv(pcols) ! 1/dp across low level divergence region real(r8) :: rhoi(pcols,0:pver) ! interface density real(r8) :: tau(pcols,-pgwv:pgwv,0:pver) ! wave Reynolds stress real(r8) :: tau0x(pcols) ! c=0 sfc. stress (zonal) real(r8) :: tau0y(pcols) ! c=0 sfc. stress (meridional) real(r8) :: ti(pcols,0:pver) ! interface temperature real(r8) :: ubi(pcols,0:pver) ! projection of wind at interfaces real(r8) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8) :: xv(pcols) ! unit vectors of source wind (x) real(r8) :: yv(pcols) ! unit vectors of source wind (y) !+waccm+ real(r8) egwdffi(pcols,0:pver) ! effective gw diffusivity at interfaces real(r8) gwut(pcols,pver,-pgwv:pgwv) ! storage for gw tendency character(len=8) :: fnum,fname ! dummy characters integer ic,l,m ! dummy integers real(r8) uzi(pcols,0:pver) ! zonal wind vertical shear at interfaces real(r8) uzm(pcols,pver) ! zonal wind vertical shear at midpoints real(r8) nzm(pcols,pver) ! vertical grad. of brunt-vaisalla freq. at midpoints real(r8) nzi(pcols,0:pver) ! vertical grad. of brunt-vaisalla freq. at interfaces real(r8) dqgw(pcols,pver) ! constituents tendencies real(r8) topflx(pcols) ! flux at top interface real(r8) :: ttgwdf(pcols,pver) ! dse tendency from diffusion term real(r8) :: ttgwke(pcols,pver) ! dse tendency from kinetic energy conversion !-waccm- real(r8) :: kvh(pcols,0:pver) ! molecular thermal diffusivity (m^2/s) real(r8) :: mkvisc ! molecular kinematic viscosity c*tint**(2/3)/rho real(r8) :: tint ! interface temperature real(r8) :: kvhmin,kvhmax !----------------------------------------------------------------------------- ncol = pcols ! ! call addfld('T_CAM',' ',' ',t,'lon',lon0,lon1,'lev',1,pver,lat) ! call addfld('U_CAM',' ',' ',u,'lon',lon0,lon1,'lev',1,pver,lat) ! call addfld('V_CAM',' ',' ',v,'lon',lon0,lon1,'lev',1,pver,lat) ! call addfld('S_CAM',' ',' ',s,'lon',lon0,lon1,'lev',1,pver,lat) ! call gw_prof(lchnk, ncol, & u ,v ,t ,pmid ,pint ,rhoi ,ni ,ti ,nm ,& !+waccm+ uzi ,nzi ,uzm ,nzm ,pcols ,pver) !-waccm- if (pgwv >0) then ! Determine the wave source for a background spectrum at ~100 mb call gw_bgnd (lchnk , ncol , & u , v , t , pmid , pint , & pdel , rpdel , lnpint , kldv , kldvmn , & ksrc , ksrcmn , rdpldv , tau , ubi , & ubm , xv , yv , PGWV , kbotbg , & rlat ,pcols , pver , lon0,lon1,lat) ! ! btf 10/8/03: Calculate kvh roughly as in vertical_diffusion.F90: ! (rhoi was returned by gw_prof above) ! Comment in vertical_diffusion.F90: ! Molecular diffusion turned on above ~60 km (50 Pa) if model top is above ~90 km ! kvh(:,:) = 0._r8 do k = 2, pver kvhmin = 1.e36 ; kvhmax = -1.e36 do i = 1, ncol tint = 0.5*(t(i,k) + t(i,k-1)) mkvisc = km_fac * tint**(2./3.) / rhoi(i,k) kvh(i,k) = kvh(i,k) + mkvisc*pr_num if (kvh(i,k) < kvhmin) kvhmin = kvh(i,k) if (kvh(i,k) > kvhmax) kvhmax = kvh(i,k) enddo ! write(6,"('gw_intr: k=',i3,' kvh(:,k)=',/,(6e12.4))") k,kvh(:,k) ! write(6,"('gw_intr: k=',i3,' kvh min,max=',2e12.4))") k,kvhmin,kvhmax enddo ! call addfld('RHOI_CAM',' ',' ',rhoi(:,1:pver),'lon',lon0,lon1, & ! 'lev',1,pver,lat) ! call addfld('KVH_CAM',' ',' ',kvh(:,1:pver),'lon',lon0,lon1, & ! 'lev',1,pver,lat) ! call addfld('ZM_CAM',' ',' ',zm(:,1:pver),'lon',lon0,lon1, & ! 'lev',1,pver,lat) ! ! When calculated at all levels, kvh range is about 25.0 to 1.6e7 call gw_drag_prof (lchnk , ncol , & PGWV , kbotbg , ktopbg , u , v , & t , pint , pdel , rpdel , lnpint , & rhoi , ni , ti , nm , dt , & kldv , kldvmn , ksrc , ksrcmn , rdpldv , & tau , ubi , ubm , xv , yv , & effgw_spec , utgw , vtgw , tau0x , tau0y , & difkgw, & !+waccm+ kvh , gwut , uzi , nzi, pcols, pver) !-waccm- ! Tau is output by gw_bgnd and gw_drag_prof -- save to timegcm secondary history: ! real(r8) :: tau(pcols,-pgwv:pgwv,0:pver) ! Reynolds stress (density*m^2/s^2) ! call addfld('TAU0_CAM',' ',' ',tau(:,0,1:pver),'lon',lon0,lon1,& ! 'lev',1,pver,lat) ! call gw_ediff (lchnk ,ncol & ! ,pgwv ,ktopbg ,kbotbg ,gwut ,ubm & ! ,uzm ,nm , nzm , t ,pint & ! ,ksrc ,egwdffi , pcols, pver) do m = 1, pcnst+pnats ! call gw_const_tend (lchnk ,ncol & ! ,pgwv ,kbotbg ,ktopbg & ! ,q(1,1,m) ,pint ,pmid & ! ,t ,dt ,rpdel & ! ,dqgw ,egwdffi ,topflx, pcols,pver ) ! ! btf 10/8/03: commented out the following zero set of dqgw: ! dqgw(:,:)=0.0 do k = 1, pver do i = 1 , ncol qtend(i,k,m) = dqgw(i,k) enddo enddo ! write(fnum,fmt='(i4)') 100+m ! fname='QTGWSP'//fnum(3:4) ! call outfld (fname,dqgw(1:pcols,1:pver),pcols,lchnk) ! write(fnum,fmt='(i4)') 100+m ! fname='QTPXGW'//fnum(3:4) ! call outfld (fname,topflx(1:pcols),pcols,lchnk) enddo ! another constituent call gw_tmpr_tend (lchnk ,ncol & ,pgwv ,kbotbg ,ktopbg ,nm ,pint ,pmid & ,t ,difkgw ,ubm ,dt ,rpdel & ,zm ,ttgw ,pcols ,pver) ! call gw_temp_tend (lchnk ,ncol & ! ,pgwv ,kbotbg ,ktopbg & ! ,s ,pint ,pmid ,t & ! ,gwut ,ubm ,dt ,rpdel & ! ,egwdffi ,ttgwdf ,ttgwke ,topflx, pcols,pver) ! call outfld ('TTPXGW',topflx(1:pcols),pcols,lchnk) !-waccm- ! add the momentum tendencies to the output tendency arrays do k = 1, pver do i = 1, ncol utend(i,k) = utgw(i,k) & !+time-gcm+ + utend(i,k) !-time-gcm- vtend(i,k) = vtgw(i,k) & !+time-gcm+ + vtend(i,k) !-time-gcm- ttend(i,k) = ttgw(i,k) & !+time-gcm+ + ttend(i,k) !-time-gcm- end do enddo do k = 0,pver do i = 1,ncol difk(i,k) = difkgw(i,k) & !+time-gcm+ + difk(i,k) !-time-gcm- end do end do ! call addfld('TTENDCAM',' ',' ',ttend,'lon',lon0,lon1,'lev',1,pver,lat) ! call addfld('UTENDCAM',' ',' ',utend,'lon',lon0,lon1,'lev',1,pver,lat) ! call addfld('VTENDCAM',' ',' ',vtend,'lon',lon0,lon1,'lev',1,pver,lat) ! Write output fields to history file ! call outfld ('UTGWSPEC', utgw , pcols, lchnk) ! call outfld ('VTGWSPEC', vtgw , pcols, lchnk) ! call outfld ('TTGWSDF', ttgwdf / cpair, pcols, lchnk) ! call outfld ('TTGWSKE', ttgwke / cpair, pcols, lchnk) !+waccm+ ! ic=0 ! do l=-pgwv,pgwv ! ic=ic+1 ! write(fnum,fmt='(i2)') ic ! fname='UTGWSPC'//fnum(2:2) ! call outfld (fname,gwut(1:pcols,1:pver,l),pcols,lchnk) ! enddo ! fname='EKGWSPEC' ! call outfld (fname,egwdffi(1:pcols,0:pver),pcols,lchnk) !-waccm- else ! zero net tendencies if no spectrum computed utend = 0. vtend = 0. ttend = 0. !+waccm+ qtend = 0. !-waccm- end if !----------------------------------------------------------------------------- ! Orographic stationary gravity wave !----------------------------------------------------------------------------- ! ! btf 10/8/03: Below orographic calls commented out for timegcm: ! ! Determine the orographic wave source ! call gw_oro (lchnk, ncol, & ! state%u , state%v , state%t , sgh , state%pmid , & ! state%pint , state%pdel , state%zm , nm , pblh , & ! kldv , kldvmn , ksrc , ksrcmn , rdpldv , & ! tau , ubi , ubm , xv , yv ) ! Solve for the drag profile ! call gw_drag_prof (lchnk, ncol, & ! 0 , kbotoro , ktoporo , state%u , state%v , & ! state%t , state%pint , state%pdel , state%rpdel, state%lnpint,& ! rhoi , ni , ti , nm , dt , & ! kldv , kldvmn , ksrc , ksrcmn , rdpldv , & ! tau , ubi , ubm , xv , yv , & ! effgw_oro , utgw , vtgw , tau0x , tau0y , & !+waccm+ ! kvh , gwut , uzi , nzi) !-waccm- ! Add the orographic tendencies to the spectrum tendencies ! Compute the temperature tendency from energy conservation (includes spectrum). ! do k = 1, pver ! do i = 1, ncol ! utend(i,k) = utend(i,k) + utgw(i,k) * landfrac(i) ! vtend(i,k) = vtend(i,k) + vtgw(i,k) * landfrac(i) ! ttend(i,k) = -(utend(i,k) * (u(i,k) + utend(i,k)*0.5*dt) & ! +vtend(i,k) * (v(i,k) + vtend(i,k)*0.5*dt)) ! ttgw(i,k) = ttend(i,k) / cpair ! end do ! end do ! Set flags for nonzero tendencies, q not yet affected by gwd ! ptend%name = "vertical diffusion" !+waccm+ ! ptend%lq(:) = .true. !-waccm- ! ptend%ls = .true. ! ptend%lu = .true. ! ptend%lv = .true. ! Write output fields to history file ! call outfld ('UTGWORO', utgw, pcols, lchnk) ! call outfld ('VTGWORO', vtgw, pcols, lchnk) ! call outfld ('TTGWORO', ttgw, pcols, lchnk) ! call outfld ('TAUGWX', tau0x, pcols, lchnk) ! call outfld ('TAUGWY', tau0y, pcols, lchnk) ! call outfld ('SGH ', sgh, pcols, lchnk) end subroutine gw_intr !=============================================================================== subroutine gw_prof (lchnk, ncol, u, v, t, pm, pi, rhoi, ni, ti, nm , & !+waccm+ uzi, nzi , uzm, nzm, pcols,pver) !-waccm- !----------------------------------------------------------------------- ! Compute profiles of background state quantities for the multiple ! gravity wave drag parameterization. ! ! The parameterization is assumed to operate only where water vapor ! concentrations are negligible in determining the density. !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures real(r8), intent(out) :: rhoi(pcols,0:pver) ! interface density real(r8), intent(out) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency real(r8), intent(out) :: ti(pcols,0:pver) ! interface temperature real(r8), intent(out) :: nm(pcols,pver) ! midpoint brunt-vaisalla frequency !+waccm+ real(r8), intent(out) :: uzi(pcols,0:pver) ! zonal wind vertical shear at interfaces real(r8), intent(out) :: uzm(pcols,pver) ! zonal wind vertical shear at midpoints real(r8), intent(out) :: nzi(pcols,0:pver) ! verical grad. of brunt-vaisalla freq. at interfaces real(r8), intent(out) :: nzm(pcols,pver) ! verical grad. of brunt-vaisalla freq. at midpoints !-waccm- !---------------------------Local storage------------------------------- integer :: i,k ! loop indexes real(r8) :: dtdp real(r8) :: n2 ! Brunt-Vaisalla frequency squared !----------------------------------------------------------------------------- ! Determine the interface densities and Brunt-Vaisala frequencies. !----------------------------------------------------------------------------- ! The top interface values are calculated assuming an isothermal atmosphere ! above the top level. k = 0 do i = 1, ncol ti(i,k) = t(i,k+1) rhoi(i,k) = pi(i,k) / (r*ti(i,k)) ni(i,k) = sqrt (g*g / (cpair*ti(i,k))) end do ! Interior points use centered differences do k = 1, pver-1 do i = 1, ncol ti(i,k) = 0.5 * (t(i,k) + t(i,k+1)) rhoi(i,k) = pi(i,k) / (r*ti(i,k)) dtdp = (t(i,k+1)-t(i,k)) / (pm(i,k+1)-pm(i,k)) n2 = g*g/ti(i,k) * (1./cpair - rhoi(i,k)*dtdp) ni(i,k) = sqrt (max (n2min, n2)) end do end do ! Bottom interface uses bottom level temperature, density; next interface ! B-V frequency. k = pver do i = 1, ncol ti(i,k) = t(i,k) rhoi(i,k) = pi(i,k) / (r*ti(i,k)) ni(i,k) = ni(i,k-1) ! write(6,"('gw_prof: i=',i3,' rhoi(i,pver) bot bndry =',e12.4)") & ! i,rhoi(i,pver) ! write(6,"('gw_prof: i=',i3,' ti(i,pver) bot bndry =',6e12.4)") & ! i,ti(i,pver) ! write(6,"('gw_prof: i=',i3,' pi(i,pver) bot bndry =',6e12.4)") & ! i,pi(i,pver) end do !----------------------------------------------------------------------------- ! Determine the midpoint Brunt-Vaisala frequencies. !----------------------------------------------------------------------------- do k=1,pver do i=1,ncol nm(i,k) = 0.5 * (ni(i,k-1) + ni(i,k)) end do end do !+waccm+ ! calculate vertical wind shear at interfaces ! assume zero shear at top and bottom ! ! btf 10/8/03: Changed "/g" to "*g" in uzi and nzi, as per Hanli and Fab: ! ! Hanli: ! >I guess this is to write dz in terms of dp. But if using dp/dz=-rho*g, it ! >seems uzi should be - uzi(i,k) *rhoi(i,k)*g. ! Fab: ! yes, you're right. Who knows why I put the "/" instead of "*" ! uzi(:,0) = 0.0 uzi(:,pver) = 0.0 do k = 1,pver-1 do i=1, ncol uzi(i,k) = (u(i,k+1) - u(i,k))/(pm(i,k+1)-pm(i,k)) uzi(i,k) = - uzi(i,k) *rhoi(i,k)*g enddo enddo ! calculate vertical gradient of n at interfaces ! assume zero gradient at top and bottom nzi(:,0) = 0.0 nzi(:,pver) = 0.0 do k = 1,pver-1 do i=1, ncol nzi(i,k) = (nm(i,k+1) - nm(i,k))/(pm(i,k+1)-pm(i,k)) nzi(i,k) = - nzi(i,k) *rhoi(i,k)*g enddo enddo ! calculate vertical wind shear at midpoints: ! average of interface values uzi uzm(:,:) = 0.0 do k = 1,pver do i=1, ncol uzm(i,k) = (uzi(i,k-1) + uzi(i,k))/2. enddo enddo ! calculate vertical gradient of n at midpoints ! assume zero gradient at top and bottom nzm(:,:) = 0.0 do k = 1,pver do i=1, ncol nzm(i,k) = (nzi(i,k-1) - nzi(i,k))/ 2. enddo enddo !-waccm- return end subroutine gw_prof !=============================================================================== subroutine gw_oro (lchnk, ncol, & u, v, t, sgh, pm, pi, dpm, zm, nm, pblh, & kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, pcols,pver) !----------------------------------------------------------------------- ! Orographic source for multiple gravity wave drag parameterization. ! ! The stress is returned for a single wave with c=0, over orography. ! For points where the orographic variance is small (including ocean), ! the returned stress is zero. !------------------------------Arguments-------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures real(r8), intent(in) :: sgh(pcols) ! standard deviation of orography real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1)) real(r8), intent(in) :: zm(pcols,pver) ! midpoint heights real(r8), intent(in) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency real(r8), intent(in) :: pblh(pcols) ! planetary boundary layer height integer, intent(out) :: kldv(pcols) ! top interface of low level stress div region integer, intent(out) :: kldvmn ! min value of kldv integer, intent(out) :: ksrc(pcols) ! index of top interface of source region integer, intent(out) :: ksrcmn ! min value of ksrc real(r8), intent(out) :: rdpldv(pcols) ! 1/dp across low level divergence region real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress real(r8), intent(out) :: ubi(pcols,0:pver) ! projection of wind at interfaces real(r8), intent(out) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(out) :: xv(pcols) ! unit vectors of source wind (x) real(r8), intent(out) :: yv(pcols) ! unit vectors of source wind (y) !---------------------------Local storage------------------------------- integer :: i,k ! loop indexes real(r8) :: pil ! don't we have pi somewhere? real(r8) :: lzsrc ! vertical wavelength at source real(r8) :: hdsp(pcols) ! surface streamline displacment height (2*sgh) real(r8) :: sghmax ! max orographic sdv to use real(r8) :: tauoro(pcols) ! c=0 stress from orography real(r8) :: zldv(pcols) ! top of the low level stress divergence region real(r8) :: nsrc(pcols) ! b-f frequency averaged over source region real(r8) :: psrc(pcols) ! interface pressure at top of source region real(r8) :: rsrc(pcols) ! density averaged over source region real(r8) :: usrc(pcols) ! u wind averaged over source region real(r8) :: vsrc(pcols) ! v wind averaged over source region !--------------------------------------------------------------------------- ! Average the basic state variables for the wave source over the depth of ! the orographic standard deviation. Here we assume that the apropiate ! values of wind, stability, etc. for determining the wave source are ! averages over the depth of the atmosphere pentrated by the typical mountain. ! Reduces to the bottom midpoint values when sgh=0, such as over ocean. !--------------------------------------------------------------------------- k = pver do i = 1, ncol ksrc(i) = k-1 psrc(i) = pi(i,k-1) rsrc(i) = pm(i,k)/(r*t(i,k)) * dpm(i,k) usrc(i) = u(i,k) * dpm(i,k) vsrc(i) = v(i,k) * dpm(i,k) nsrc(i) = nm(i,k)* dpm(i,k) hdsp(i) = 2.0 * sgh(i) end do do k = pver-1, pver/2, -1 do i = 1, ncol if (hdsp(i) > sqrt(zm(i,k)*zm(i,k+1))) then ksrc(i) = k-1 psrc(i) = pi(i,k-1) rsrc(i) = rsrc(i) + pm(i,k) / (r*t(i,k))* dpm(i,k) usrc(i) = usrc(i) + u(i,k) * dpm(i,k) vsrc(i) = vsrc(i) + v(i,k) * dpm(i,k) nsrc(i) = nsrc(i) + nm(i,k)* dpm(i,k) end if end do end do do i = 1, ncol rsrc(i) = rsrc(i) / (pi(i,pver) - psrc(i)) usrc(i) = usrc(i) / (pi(i,pver) - psrc(i)) vsrc(i) = vsrc(i) / (pi(i,pver) - psrc(i)) nsrc(i) = nsrc(i) / (pi(i,pver) - psrc(i)) ubi(i,pver) = sqrt (usrc(i)**2 + vsrc(i)**2) xv(i) = usrc(i) / ubi(i,pver) yv(i) = vsrc(i) / ubi(i,pver) end do ! Project the local wind at midpoints onto the source wind. do k = 1, pver do i = 1, ncol ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i) end do end do ! Compute the interface wind projection by averaging the midpoint winds. ! Use the top level wind at the top interface. do i = 1, ncol ubi(i,0) = ubm(i,1) end do do k = 1, pver-1 do i = 1, ncol ubi(i,k) = 0.5 * (ubm(i,k) + ubm(i,k+1)) end do end do !--------------------------------------------------------------------------- ! Determine the depth of the low level stress divergence region, as ! the max of the source region depth and 1/2 the vertical wavelength at the ! source. !--------------------------------------------------------------------------- pil = acos(-1.) do i = 1, ncol lzsrc = max(2. * pil * usrc(i) / nsrc(i), lzmax) zldv(i) = max(hdsp(i), 0.5 * lzsrc) end do ! find the index of the top of the low level divergence region kldv(:) = pver-1 do k = pver-1, pver/2, -1 do i = 1, ncol if (zldv(i) .gt. sqrt(zm(i,k)*zm(i,k+1))) then kldv(i) = k-1 end if end do end do ! Determine the orographic c=0 source term following McFarlane (1987). ! Set the source top interface index to pver, if the orographic term is zero. do i = 1, ncol if ((ubi(i,pver) .gt. orovmin) .and. (hdsp(i) .gt. orohmin)) then sghmax = fcrit2 * (ubi(i,pver) / nsrc(i))**2 tauoro(i) = oroko2 * min(hdsp(i)**2, sghmax) * rsrc(i) * nsrc(i) * ubi(i,pver) else tauoro(i) = 0. ksrc(i) = pver kldv(i) = pver end if end do ! Set the phase speeds and wave numbers in the direction of the source wind. ! Set the source stress magnitude (positive only, note that the sign of the ! stress is the same as (c-u). do i = 1, ncol tau(i,0,pver) = tauoro(i) end do ! Determine the min value of kldv and ksrc for limiting later loops ! and the pressure at the top interface of the low level stress divergence ! region. ksrcmn = pver kldvmn = pver do i = 1, ncol ksrcmn = min(ksrcmn, ksrc(i)) kldvmn = min(kldvmn, kldv(i)) if (kldv(i) .ne. pver) then rdpldv(i) = 1. / (pi(i,kldv(i)) - pi(i,pver)) end if end do if (fracldv .le. 0.) kldvmn = pver return end subroutine gw_oro !=============================================================================== subroutine gw_bgnd (lchnk, ncol, & u, v, t, pm, pi, dpm, rdpm, piln, & kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, & ngwv, kbot, rlat, pcols, pver, lon0,lon1,lat) !----------------------------------------------------------------------- ! Driver for multiple gravity wave drag parameterization. ! ! The parameterization is assumed to operate only where water vapor ! concentrations are negligible in determining the density. !----------------------------------------------------------------------- ! use phys_grid, only: get_rlat_all_p !------------------------------Arguments-------------------------------- use gw_share,only: faca,faca1,facnh,facsh,facoron,facoros, & wspnh,wspsh,uth use init_module,only: glon ! integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: lat ! current latitude index integer, intent(in) :: lon0,lon1 ! subdomain lons integer, intent(in) :: kbot ! index of bottom (source) interface integer, intent(in) :: ngwv ! number of gravity waves to use real(r8), intent(in) :: u(pcols,pver+1) ! midpoint zonal wind real(r8), intent(in) :: v(pcols,pver+1) ! midpoint meridional wind real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressures real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1)) real(r8), intent(in) :: rdpm(pcols,pver) ! 1. / (pi(k)-pi(k-1)) real(r8), intent(in) :: piln(pcols,0:pver) ! ln(interface pressures) real(r8), intent(in) :: rlat ! current latitude in radians integer, intent(out) :: kldv(pcols) ! top interface of low level stress divergence region integer, intent(out) :: kldvmn ! min value of kldv integer, intent(out) :: ksrc(pcols) ! index of top interface of source region integer, intent(out) :: ksrcmn ! min value of ksrc real(r8), intent(in) :: rdpldv(pcols) ! 1/dp across low level divergence region real(r8), intent(out) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress real(r8), intent(out) :: ubi(pcols,0:pver) ! projection of wind at interfaces real(r8), intent(out) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(out) :: xv(pcols) ! unit vectors of source wind (x) real(r8), intent(out) :: yv(pcols) ! unit vectors of source wind (y) !---------------------------Local storage------------------------------- integer :: i,k,l ! loop indexes integer :: idir(pcols) ! direction indicator integer :: l1 ! direction indicator real(r8) :: tauback ! max background stress real(r8) :: tauback1 ! max background stress with local time dependence real(r8) :: taueq ! max background stress at the equator real(r8) :: usrc(pcols) ! u wind averaged over source region real(r8) :: vsrc(pcols) ! v wind averaged over source region real(r8) :: al0 ! Used in lat dependence of GW spec. real(r8) :: al1 ! Used in lat dependence of GW spec. real(r8) :: al2 ! Used in lat dependence of GW spec. real(r8) :: dlat0 ! Used in lat dependence of GW spec. real(r8) :: flat_gw ! The actual lat dependence of GW spec. real(r8) :: flat_gw1 ! this one cut the high latitude amplitude real(r8) :: pi_g ! 3.14........ real(r8) :: taueqw ! taueq latitudinal width real(r8) :: wspd ! phase speed spectral peak real(r8) :: whw ! phase speed spectral width real(r8) :: wspde ! phase speed spectral peak for the equatorial distribution real(r8) :: whwe ! phase speed spectral width for the equatorial distribution real(r8) :: slt ! solar local time (hours) real(r8) :: flt ! gravity wave distribution with local time real(r8) :: tautmp ! temperory variable real(r8) :: xran ! local random number real,external :: fslt ! ! 3/5/07 btf: ! If random==true, then use random number generator for lt dependence ! in gw_bgnd ! ! logical :: random = .true. logical :: random = .false. real :: fxran(pcols) real :: fxran_diag(nlev,lon0:lon1) ! diag only real :: myglon(nlonp4) integer :: mynlon=nlonp4-4 myglon(3:nlonp4-2) = glon(1:mynlon) myglon(1:2) = myglon(mynlon+1:mynlon+2) myglon(nlonp4-1:nlonp4) = myglon(3:4) !--------------------------------------------------------------------------- ! Determine the source layer wind and unit vectors, then project winds. !--------------------------------------------------------------------------- ! Just use the source level interface values for the source ! wind speed and direction (unit vector). idir = 1 k = kbot do i = 1, ncol ksrc(i) = k kldv(i) = k usrc(i) = 0.5*(u(i,k-1)+u(i,k)) vsrc(i) = 0.5*(v(i,k-1)+v(i,k)) ubi(i,kbot) = sqrt (usrc(i)**2 + vsrc(i)**2) xv(i) = usrc(i) / ubi(i,k) yv(i) = vsrc(i) / ubi(i,k) if(xv(i).le.-0.1) idir(i) = -1 end do ! Project the local wind at midpoints onto the source wind. do k = 1, kbot do i = 1, ncol ubm(i,k) = u(i,k) * xv(i) + v(i,k) * yv(i) end do end do ! Compute the interface wind projection by averaging the midpoint winds. ! Use the top level wind at the top interface. do i = 1, ncol ubi(i,0) = ubm(i,1) end do do k = 1, kbot-1 do i = 1, ncol ubi(i,k) = 0.5 * (ubm(i,k) + ubm(i,k+1)) end do end do !----------------------------------------------------------------------- ! Gravity wave sources !----------------------------------------------------------------------- pi_g=4.*atan(1.) al0=20.*pi_g/180. dlat0=10.*pi_g/180. al1=80.*pi_g/180. al2 = 45.*pi_g/180. ! Determine the background stress at c=0 tauback = taubgnd * tauscal * faca ! taueqw = 10.*3.14159/180. taueqw = 20.*3.14159/180. taueq = taubgnd * tauscal*faca1*exp(-(rlat/taueqw)**2) wspd= (0.5*(1.+tanh(rlat/dlat0))*wspnh & + 0.5*(1.+tanh(-rlat/dlat0))*wspsh) ! *abs(sin(2.*rlat)) ! wspd = wspnh ! if(abs(rlat) .gt. al2) wspd = wspnh*abs(sin(2.*rlat)) whw = 50. whwe = 50. wspde = 5. ! wspde = 50. ! Include dependence on latitude: ! The lat function was obtained by RR Garcia as ! currently used in his 2D model ! [Added by F. Sassi on May 30, 1997] ! ! call get_rlat_all_p(lchnk, ncol, rlat) ! rlat is passed in from timegcm flat_gw= 0.5*(1.+tanh((rlat-al0)/dlat0))*facnh & + 0.5*(1.+tanh(-(rlat+al0)/dlat0))*facsh ! *abs(sin(2.*rlat)) ! flat_gw=max(flat_gw,0.05) ! flat_gw = 1. tauback1=tauback*flat_gw ! ! 3/5/07 btf: Set logical random=.true. (module data above) to use random ! number generator for local time dependence (see flt below) ! if (random) then ! use random number generator for lt dependence do i=1,ncol call random_number(xran) fxran(i) = xran*2.0 enddo do k=1,nlev fxran_diag(k,:) = fxran(1:ncol) enddo ! ! Note tgcmproc_f90 has problems plotting fxran, but tgcmproc_idl does plot it. call addfld('FXRAN',' ',' ',fxran_diag,'lev',1,nlev,'lon',lon0,lon1,lat) endif ! ! Set the phase speeds and wave numbers in the direction of the source wind. ! Set the source stress magnitude (positive only, note that the sign of the ! stress is the same as (c-u). ! flt = 1. ! no local time dependence do l = 1, ngwv do i = 1, ncol l1 = l*idir(i) ! slt = solar local time (function fslt is in util.F) slt = fslt(slt,uth,myglon(lon0+i-1),1) if (random) flt = fxran(i) ! use random for local time dependence tauback1 = tauback*flt taueq = taueq*flt tau(i,l1,kbot) = (exp(-((c(l1)-idir(i)*wspd)/whw)**2) & *tauback1 & +exp(-((c(l1)-idir(i)*wspde)/whwe)**2) & *taueq)*abs(xv(i)) & +tauback1*exp(-(c(l1)/whw)**2) & *abs(yv(i))*0.01 tau(i,-l1,kbot) = (exp(-((c(l1)+idir(i)*wspd)/whw)**2) & *tauback1 & +exp(-((c(l1)+idir(i)*wspde)/whwe)**2) & *taueq)*abs(xv(i)) & +tauback1*exp(-(c(l1)/whw)**2) & *abs(yv(i))*0.01 end do end do do l = 0,0 do i = 1, ncol l1 = 0 flt = 1. ! no local time dependence ! tauback1 = tauback*flt tauback1 = tauback tau(i,l1,kbot) = (exp(-(c(l1)/whw)**2) & *tauback1)*abs(xv(i)) & +tauback1*exp(-(c(l1)/whw)**2) & *abs(yv(i))*0.01 enddo enddo if(rlat.ge.0) tau(:,0,kbot) = tau(:,0,kbot)*facoron if(rlat.lt.0) tau(:,0,kbot) = tau(:,0,kbot)*facoros ! ! At standard 5 deg horizontal resolution, use orography data ! sghdat (see orogdat.F), as in tgcm24: ! ! if (dlon == 5.0) then ! tautmp = 1.e-4 ! tautmp = 1.e-6 ! ! if(rlat.gt.0.6981.and.rlat.lt.1.57) tautmp = 1.4e-3 & ! *exp(-(rlat-1.134)**2/0.2684**2)*facoron ! if(rlat.gt.-1.57.and.rlat.lt.-0.6981) tautmp = 1.4e-3 & ! *exp(-(rlat+1.134)**2/0.2684**2)*facoros ! ! do i = 1, ncol ! ! slt = solar local time (function fslt is in util.F) ! slt = fslt(slt,uth,myglon(lon0+i-1),1) ! flt = 1.0 ! if(rlat.le.-1.2) flt = 0.1 ! tau(i,0,kbotoro) = tautmp+2.0e-4*sghdat(i,lat)/1240.*flt ! tau(i,0,kbotoro) = tautmp+1.0e-3*sghdat(i,lat)/1240.*flt ! tau(i,0,kbotoro) = tautmp+2.0e-3*sghdat(i,lat)/1240.*flt ! tau(i,0,kbotoro) = tautmp+4.0e-3*sghdat(i,lat)/1240.*flt ! end do ! write(6,"('gw_bgnd: tau(:,0,kbotoro)=',/,(6e12.4))") tau(:,0,kbotoro) ! ! Non-standard horizontal resolution (probably 2.5) -> cannot use sghdat. ! else ! ! Orographic waves leaked above 30mb. ! ! write(6,"('gw_bgnd: NOT using orography data sghdat..')") tautmp = 1.e-6 if(rlat.gt.0.6981.and.rlat.lt.1.57) & tautmp = 1.0e-3*exp(-(rlat-1.134)**2/0.2684**2)*facoron if(rlat.gt.-1.57.and.rlat.lt.-0.6981) & tautmp = 1.0e-3*exp(-(rlat+1.134)**2/0.2684**2)*facoros do i = 1,ncol tau(i,0,kbot) = max(tautmp,1.e-6)+tau(i,0,kbot) enddo ! Determine the min value of kldv and ksrc for limiting later loops ! and the pressure at the top interface of the low level stress divergence ! region. ksrcmn = pver kldvmn = pver return end subroutine gw_bgnd !=============================================================================== subroutine gw_drag_prof (lchnk ,ncol , & ngwv ,kbot ,ktop ,u ,v , & t ,pi ,dpm ,rdpm ,piln , & rhoi ,ni ,ti ,nm ,dt , & kldv ,kldvmn ,ksrc ,ksrcmn,rdpldv, & tau ,ubi ,ubm ,xv ,yv , & effgw ,ut ,vt ,tau0x ,tau0y , & difk, & !+waccm+ kvh ,gwut ,uzi ,nzi, pcols,pver) !-waccm- !----------------------------------------------------------------------- ! Solve for the drag profile from the multiple gravity wave drag ! parameterization. ! 1. scan up from the wave source to determine the stress profile ! 2. scan down the stress profile to determine the tendencies ! => apply bounds to the tendency ! a. from wkb solution ! b. from computational stability constraints ! => adjust stress on interface below to reflect actual bounded tendency !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: kbot ! index of bottom (source) interface integer, intent(in) :: ktop ! index of top interface of gwd region integer, intent(in) :: ngwv ! number of gravity waves to use integer, intent(in) :: kldv(pcols) ! top interface of low level stress divergence region integer, intent(in) :: kldvmn ! min value of kldv integer, intent(in) :: ksrc(pcols) ! index of top interface of source region integer, intent(in) :: ksrcmn ! min value of ksrc real(r8), intent(in) :: u(pcols,pver) ! midpoint zonal wind real(r8), intent(in) :: v(pcols,pver) ! midpoint meridional wind real(r8), intent(in) :: t(pcols,pver) ! midpoint temperatures real(r8), intent(in) :: pi(pcols,0:pver) ! interface pressures real(r8), intent(in) :: dpm(pcols,pver) ! midpoint delta p (pi(k)-pi(k-1)) real(r8), intent(in) :: rdpm(pcols,pver) ! 1. / (pi(k)-pi(k-1)) real(r8), intent(in) :: piln(pcols,0:pver) ! ln(interface pressures) real(r8), intent(in) :: rhoi(pcols,0:pver) ! interface density real(r8), intent(in) :: ni(pcols,0:pver) ! interface Brunt-Vaisalla frequency real(r8), intent(in) :: ti(pcols,0:pver) ! interface temperature real(r8), intent(in) :: nm(pcols,pver) ! midpoint Brunt-Vaisalla frequency real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: rdpldv(pcols) ! 1/dp across low level divergence region real(r8), intent(in) :: ubi(pcols,0:pver) ! projection of wind at interfaces real(r8), intent(in) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(in) :: xv(pcols) ! unit vectors of source wind (x) real(r8), intent(in) :: yv(pcols) ! unit vectors of source wind (y) real(r8), intent(in) :: effgw ! tendency efficiency real(r8), intent(inout) :: tau(pcols,-pgwv:pgwv,0:pver)! wave Reynolds stress real(r8), intent(out) :: ut(pcols,pver) ! zonal wind tendency real(r8), intent(out) :: vt(pcols,pver) ! meridional wind tendency real(r8), intent(out) :: difk(pcols,0:pver) ! eddy diffusion coefficient due to wave breaking ! real(r8), intent(out) :: difkm(pcols,0:pver) ! difk+molecular diffusion real(r8), intent(out) :: tau0x(pcols) ! c=0 sfc. stress (zonal) real(r8), intent(out) :: tau0y(pcols) ! c=0 sfc. stress (meridional) !+waccm+ real(r8), intent(in) :: kvh(pcols,0:pver) ! molecular thermal diffusivity real(r8), intent(in) :: uzi(pcols,0:pver) ! vertical shear at interfaces real(r8), intent(in) :: nzi(pcols,0:pver) ! vertical shear at interfaces real(r8), intent(out) :: gwut(pcols,pver,-pgwv:pgwv) ! storage for gw tendency !-waccm- !---------------------------local storage------------------------------- integer :: i,k,l ! loop indexes real(r8) :: d(pcols) ! "total" diffusivity real(r8) :: dsat(pcols,-pgwv:pgwv) ! saturation diffusivity real(r8) :: dscal ! fraction of dsat to use real(r8) :: mi ! imaginary part of vertical wavenumber real(r8) :: taudmp ! stress after damping real(r8) :: taumax(pcols) ! max(tau) for any l real(r8) :: tausat(pcols,-pgwv:pgwv) ! saturation stress real(r8) :: ubmc(pcols,-pgwv:pgwv) ! (ub-c) real(r8) :: ubmc2 ! (ub-c)**2 real(r8) :: ubt(pcols,pver) ! ubar tendency real(r8) :: ubtl ! ubar tendency from wave l real(r8) :: ubtlsat ! saturation tendency real(r8) :: wrk real(r8) :: difkm(pcols,0:pver) ! difk+molecular diffusion ! Initialize gravity wave drag tendencies to zero do k=1,pver do i=1,pcols ut(i,k) = 0. vt(i,k) = 0. end do end do !+waccm+ gwut(:,:,:)=0.0 !-waccm- !--------------------------------------------------------------------------- ! Compute the stress profiles and diffusivities !--------------------------------------------------------------------------- ! Loop from bottom to top to get stress profiles do k = kbot-1, ktop, -1 ! Determine the absolute value of the saturation stress and the diffusivity ! for each wave. ! Define critical levels where the sign of (u-c) changes between interfaces. d(:ncol) = dback do l = -ngwv, ngwv do i = 1, ncol ubmc(i,l) = ubi(i,k) - c(l) tausat(i,l) = abs (effkwv * rhoi(i,k) * ubmc(i,l)**3 / (2.*ni(i,k)) ) if (tausat(i,l) .le. taumin) tausat(i,l) = 0.0 if (ubmc(i,l) / (ubi(i,k+1) - c(l)) .le. 0.0) tausat(i,l) = 0.0 !+waccm+ ! ! btf 10/8/03: changed dsat statement, as per Hanli and Fab ! Hanli wrote: ! >I don't quite understand how alpha gets in there. For one thing, alpha is ! >newtonian cooling and its unit is something like 1/s, while the rest of ! >the terms in the parenthesis have unit 1/m. ! Fab wrote: ! comment that line block, uncomment the following, and you'll be fine. ! ! Commented this dsat: ! dsat(i,l) = effkwv * ubmc(i,l)**4 / ni(i,k)**3 & ! *(0.5 / (rog * ti(i,k)) - 1.5 * uzi(i,k) / ubmc(i,l) & ! +0.5 * nzi(i,k)/ni(i,k) - alpha(k) ) !-waccm- ! ! Uncommented this dsat: dsat(i,l) = (ubmc(i,l) / ni(i,k))**2 * & (effkwv * ubmc(i,l)**2 / (rog * ti(i,k) * ni(i,k))*.5 - alpha(k)) dscal = min (1.0_r8, tau(i,l,k+1) / (tausat(i,l)+taumin)) !+time-gcm+ d(i) = max( d(i), dscal * dsat(i,l)) + kvh(i,k) difkm(i,k) = d(i) difk(i,k) = max((difkm(i,k) - kvh(i,k)),0.0) ! Note this difk doesn't contain molecular ! diffusion and different from the waccm version !-time-gcm- end do end do ! Compute stress for each wave. The stress at this level is the min of ! the saturation stress and the stress at the level below reduced by damping. ! The sign of the stress must be the same as at the level below. do l = -ngwv, ngwv do i = 1, ncol ubmc2 = max(ubmc(i,l)**2, ubmc2mn) mi = ni(i,k) / (2. * kwv * ubmc2) * (alpha(k) + ni(i,k)**2/ubmc2 * d(i)) wrk = -2.*mi*rog*t(i,k+1)*(piln(i,k+1) - piln(i,k)) if( wrk >= -150._r8 ) then taudmp = tau(i,l,k+1) * exp( wrk ) else taudmp = 0._r8 end if if (taudmp .le. taumin) taudmp = 0. tau(i,l,k) = min (taudmp, tausat(i,l)) end do end do ! The orographic stress term does not change across the source region ! Note that k ge ksrcmn cannot occur without an orographic source term if (k .ge. ksrcmn) then do i = 1, ncol if (k .ge. ksrc(i)) then tau(i,0,k) = tau(i,0,pver) end if end do end if ! Require that the orographic term decrease linearly (with pressure) ! within the low level stress divergence region. This supersedes the above ! requirment of constant stress within the source region. ! Note that k ge kldvmn cannot occur without an orographic source term, since ! kldvmn=pver then and k<=pver-1 if (k .ge. kldvmn) then do i = 1, ncol if (k .ge. kldv(i)) then tau(i,0,k) = min (tau(i,0,k), tau(i,0,pver) * & (1. - fracldv * (pi(i,k)-pi(i,pver)) * rdpldv(i))) end if end do end if ! Apply lower bounds to the stress if ngwv > 0. if (ngwv .ge. 1) then ! Determine the max value of tau for any l do i = 1, ncol taumax(i) = tau(i,-ngwv,k) end do do l = -ngwv+1, ngwv do i = 1, ncol taumax(i) = max(taumax(i), tau(i,l,k)) end do end do do i = 1, ncol taumax(i) = mxrange * taumax(i) end do ! Set the min value of tau for each wave to the max of mxrange*taumax or ! mxasym*tau(-c) do l = 1, ngwv do i = 1, ncol tau(i, l,k) = max(tau(i, l,k), taumax(i)) tau(i, l,k) = max(tau(i, l,k), mxasym*tau(i,-l,k)) tau(i,-l,k) = max(tau(i,-l,k), taumax(i)) tau(i,-l,k) = max(tau(i,-l,k), mxasym*tau(i, l,k)) end do end do l = 0 do i = 1, ncol tau(i,l,k) = max(tau(i,l,k), mxasym * 0.5 * (tau(i,l-1,k) + tau(i,l+1,k))) end do end if end do difk(:,kbot) = difk(:,kbot-1) ! Put an upper bound on the stress at the top interface to force some stress ! divergence in the top layer. This prevents the easterlies from running ! away in the summer mesosphere, since most of the gravity wave activity ! will pass through a top interface at 75--80 km under strong easterly ! conditions. !++BAB fix to match ccm3.10 !!$ do l = -ngwv, ngwv !!$ do i = 1, ncol !!$ tau(i,l,0) = min(tau(i,l,0), 0.5*tau(i,l,1)) !!$ end do !!$ end do !--BAB fix to match ccm3.10 !--------------------------------------------------------------------------- ! Compute the tendencies from the stress divergence. !--------------------------------------------------------------------------- ! Loop over levels from top to bottom do k = ktop+1, kbot ! Accumulate the mean wind tendency over wavenumber. do i = 1, ncol ubt (i,k) = 0.0 end do do l = -ngwv, ngwv do i = 1, ncol ! Determine the wind tendency including excess stress carried down from above. ubtl = g * (tau(i,l,k)-tau(i,l,k-1)) * rdpm(i,k) ! Require that the tendency be no larger than the analytic solution for ! a saturated region [proportional to (u-c)^3]. ubtlsat = effkwv * abs((c(l)-ubm(i,k))**3) /(2.*rog*t(i,k)*nm(i,k)) ubtl = min(abs(ubtl), ubtlsat) ! Apply tendency limits to maintain numerical stability. ! 1. du/dt < |c-u|/dt so u-c cannot change sign (u^n+1 = u^n + du/dt * dt) ! 2. du/dt < tndmax so that ridicuously large tendency are not permitted ubtl = min(ubtl, umcfac * abs(c(l)-ubm(i,k)) / dt) ubtl = min(ubtl, tndmax) ! accumulate the mean wind tendency over wavenumber. !+waccm+ gwut(i,k,l) = sign(ubtl, c(l)-ubm(i,k)) * effgw !-waccm- ubt (i,k) = ubt (i,k) + sign(ubtl, c(l)-ubm(i,k)) ! Redetermine the effective stress on the interface below from the wind ! tendency. If the wind tendency was limited above, then the new stress ! will be small than the old stress and will cause stress divergence in ! the next layer down. This has the effect of smoothing large stress ! divergences downward while conserving total stress. tau(i,l,k) = tau(i,l,k-1) + ubtl * dpm(i,k) / g end do end do ! Project the mean wind tendency onto the components and scale by "efficiency". do i = 1, ncol ! ut(i,k) = ubt(i,k) * xv(i) * effgw ! vt(i,k) = ubt(i,k) * yv(i) * effgw ut(i,k) = ubt(i,k) * xv(i) vt(i,k) = ubt(i,k) * yv(i) end do ! End of level loop end do !----------------------------------------------------------------------- ! Project the c=0 stress (scaled) in the direction of the source wind for recording ! on the output file. !----------------------------------------------------------------------- do i = 1, ncol tau0x(i) = tau(i,0,kbot) * xv(i) * effgw tau0y(i) = tau(i,0,kbot) * yv(i) * effgw end do return end subroutine gw_drag_prof !======================================================================================= subroutine gw_const_tend (lchnk ,ncol & ,ngwv ,kbot ,ktop ,q ,pi ,pm & ,t ,dt ,rdpm ,dq ,egwdffi ,topflx & ,pcols, pver) ! ! calculates constituent tendencies due to gw breaking ! ! method: ! a constituent flux on interfaces is given by ! ! ! rho * (w'q') = rho * Deff qz ! ! ! where (all evaluated on interfaces): ! ! rho = density ! qz = constituent vertical gradient ! Deff = effective diffusivity ! ! ! an effective diffusivity is calculated by adding up the diffusivities from all waves ! (see gw_ediff) ! tendency is calculated by invoking lu decomposition and solving as for a regular ! diffusion equation. ! ! author : sassi - jan 2001 !--------------------------------------------------------------------------------- ! B.Foster: not available in timegcm: ! use vertical_diffusion, only: vd_lu_decomp, vd_lu_solve, cfluxtop implicit none !---------------------------input arguments--------------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: kbot ! index of bottom (source) midpoint integer, intent(in) :: ktop ! index of top midpoint of gwd region integer, intent(in) :: ngwv ! number of gravity waves to use real(r8), intent(in) :: q(pcols,pver) ! constituent real(r8), intent(in) :: t(pcols,pver) ! temperature at midpoints real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: rdpm(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: pi(pcols,pver+1) ! interface pressure real(r8), intent(in) :: pm(pcols,pver+1) ! midpoint pressure real(r8), intent(in) :: egwdffi(pcols,pver+1) ! effective gw diffusivity at interfaces !--------------------------output arguments------------------------------------- real(r8), intent(out) :: dq(pcols,pver) ! constituent tendencies at midpoints real(r8), intent(out) :: topflx(pcols) ! constituent flux at top interface !--------------------------local worspace--------------------------------------- integer i,k,l real(r8) tmpi2(pcols,pver+1) ! factor real(r8) ca(pcols,pver) ! -upper diagonal real(r8) cc(pcols,pver) ! -lower diagonal real(r8) term(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) real(r8) ze(pcols,pver) ! term in tri-diag. matrix system real(r8) qconst(pcols,pver) ! temporary storage for constituent real(r8) rho(pcols,pver+1) real(r8) zero(pcols) !------------------------------------------------------------------------------- zero = 0.0 dq(:,:) = 0.0 !$$$ if (ngwv.eq.0) return ! Compute factor for diffusion matrices at interfaces (except top), ! Compute rho at interfaces p/RT, Ti = (Tm_k + Tm_k-1)/2, interface k- do k = 2, pver do i = 1, ncol rho(i,k) = pi(i,k) * 2. / (r * (t(i,k) + t(i,k-1))) end do end do do i = 1, ncol rho(i,pver+1) = pi (i,pver+1) / (r * t(i,pver)) end do do i = 1, ncol rho(i,1) = pi (i,1) / (r * t(i,1)) end do ! Calculate dt * (g*rho)^2/dp at interior interfaces, interface k- tmpi2(:,:)=0.0 do k = ktop+2 , kbot+1 do i = 1, ncol tmpi2(i,k) = dt * (g * rho(i,k))**2 / (pm(i,k) - pm(i,k-1)) end do end do qconst = 0.0 do k = 1, pver do i = 1, ncol qconst(i,k) = q(i,k) enddo enddo ! Constant flux approximation at top interface topflx(:) = 0.0 !!$ call cfluxtop ( & !!$ lchnk ,ncol , & !!$ qconst ,egwdffi ,pm ,rdpm ,dt , & !!$ pi ,t ,topflx) ! Decompose and solve the diffusion matrix ! call vd_lu_decomp( & ! lchnk ,ncol , & ! zero ,egwdffi ,tmpi2 ,rdpm ,dt , & ! ca ,cc ,term ,ze ,ktop+1 , & ! kbot+1 ) ! call vd_lu_solve( & ! lchnk ,ncol , & ! qconst(1,1),ca ,ze ,term , & ! ktop+1 ,kbot+1 ) ! Evaluate tendency to be reported back do k = ktop+1, kbot do i = 1, ncol dq(i,k) = (qconst(i,k) - q(i,k)) / dt enddo enddo return end subroutine gw_const_tend !===================================================================================== subroutine gw_ediff (lchnk ,ncol & ,ngwv ,ktop ,kbot ,gwut ,ubm & ,uzm ,nm ,nzm ,t ,pi & ,ksrc ,egwdffi ,pcols ,pver) ! call gw_ediff (lchnk ,ncol & ! ,pgwv ,ktopbg ,kbotbg ,gwut ,ubm & ! ,uzm ,nm , nzm , t ,pint & ! ,ksrc ,egwdffi , pcols, pver) ! ! Calculate effective diffusivity associated with GW forcing. ! ! Author: F. Sassi, Jan 31, 2001 ! implicit none !-------------------------------Input Arguments------------------------------------------ integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: kbot ! index of bottom (source) midpoint integer, intent(in) :: ktop ! index of top midpoint of gwd region integer, intent(in) :: ngwv ! number of gravity waves to use integer, intent(in) :: ksrc(pcols) ! index of top interface of source region real(r8), intent(in) :: t(pcols,pver) ! midpoint temperature real(r8), intent(in) :: pi(pcols,pver+1) ! interfaces pressure real(r8), intent(in) :: gwut(pcols,pver,-pgwv:pgwv) ! GW zonal wind tendencies at midpoint real(r8), intent(in) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(in) :: uzm(pcols,pver) ! zonal wind shear at midpoints real(r8), intent(in) :: nm(pcols,pver) ! Brunt-Vaisalla frequency real(r8), intent(in) :: nzm(pcols,pver) ! verical grad. of brunt-vaisalla freq. at midpoints !-----------------------------Output Argument------------------------------------------- real(r8), intent(out) :: egwdffi(pcols,0:pver) ! effective gw diffusivity at interfaces !-----------------------------Local workspace------------------------------------------- real(r8),parameter :: delta_min=0. real(r8) delta ! dissipation rate real(r8) egwdffm(pcols,pver) ! effective gw diffusivity at midpoints real(r8),parameter :: egwdffm_min=150. integer i,k,l real(r8), parameter :: prndl=1. ! Prandtl no. real(r8) :: coef real(r8) :: nms2(pver) ! zonal average nm^2 real(r8) :: dummy(pver) ! dummy real(r8) :: rho(pcols,pver+1) ! density at interfaces real(r8) :: dscale ! density scale height !--------------------------------------------------------------------------------------- ! zero egwdffi(:,:) = 0. egwdffm(:,:) = 0. if (ngwv == 0) return do k=2, pver do i=1, ncol rho(i,k)= pi(i,k) / (r * 0.5 * (t(i,k-1) + t(i,k)) ) enddo enddo do i=1, ncol rho(i,1) = pi(i,1) / (r * t(i,1)) rho(i,pver+1) = pi(i,pver+1) / (r * t(i,pver)) enddo ! Calculate effective diffusivity at midpoints do k = ktop+1, kbot do i = 1, ncol do l = -ngwv, ngwv dscale = g * ( rho(i,k) - rho(i,k+1) ) / (pi(i,k) - pi(i,k+1)) dscale = 1. / dscale delta = gwut(i,k,l) / ( (c(l)-ubm(i,k)) * (1. - dscale * nzm(i,k) / & nm(i,k)) + 3. * dscale * uzm(i,k) ) delta=max(delta,delta_min) egwdffm(i,k) = egwdffm(i,k) + & prndl * 0.5 * (ubm(i,k)-c(l))**2 * delta / nm(i,k)**2 enddo ! another wave enddo ! another long enddo ! another level ! limit diffusivity within some reasonable range do i=1, ncol do k=ktop+1, kbot egwdffm(i,k)=min(egwdffm_min,egwdffm(i,k)) enddo enddo ! Interpolate effective diffusivity to interfaces do i = 1, ncol ! Assume zero at top and bottom interfaces do k=ktop+1, kbot-1 egwdffi(i,k) = (egwdffm(i,k)+egwdffm(i,k+1)) / 2. enddo enddo ! Do not calculate diffusivities below source level do k= ktop,kbot do i = 1, ncol if (k .ge. ksrc(i)) egwdffi(i,k) = 0.0 enddo enddo return end subroutine gw_ediff !======================================================================================== subroutine gw_tmpr_tend (lchnk ,ncol & ,ngwv ,kbot ,ktop ,nm ,pi ,pm & ,t ,difkgw ,ubm ,dt ,rdpm & ,zm ,ttend , pcols, pver) ! ! Calculates temperature tendency due to GW breaking based on Liu (2000). ! ! ! Author : Liu (adapted from TIME-GCM: mgw.F) !--------------------------------------------------------------------------------- ! B.Foster: not available in timegcm: ! use vertical_diffusion, only: vd_lu_decomp, vd_lu_solve, cfluxtop use gw_share,only: prndtlt implicit none !---------------------------Input Arguments--------------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: kbot ! index of bottom (source) midpoint integer, intent(in) :: ktop ! index of top midpoint of gwd region integer, intent(in) :: ngwv ! number of gravity waves to use real(r8), intent(in) :: difkgw(pcols,0:pver) ! eddy diffusion coefficient due to wave breaking real(r8), intent(in) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(in) :: t(pcols,pver) ! temperature at midpoints real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: rdpm(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: pi(pcols,pver+1) ! interface pressure real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressure real(r8), intent(in) :: zm(pcols,pver) ! geopotential heigh (meters) on midpoints real(r8), intent(in) :: nm(pcols,pver) ! midpoint Brunt-Vasaila frequency !--------------------------Output Arguments------------------------------------- real(r8), intent(out) :: ttend(pcols,pver) ! T tendencies at midpoints (J/kg/s) !--------------------------Local Worspace--------------------------------------- integer i,k,m,l real(r8) tmpi2(pcols,pver+1) ! pi(k)/(.5*(tm1(k)+tm1(k-1)) real(r8) hflx(pcols,pver) ! heat flux induced by wave breaking. real(r8) ttenda(pcols,pver) ! temperature tendency part a (K/s) real(r8) ttendb(pcols,pver) ! temperature tendency part b (K/s) !------------------------------------------------------------------------------- !+time-gcm+ do k=ktop+1,kbot hflx(:,k) = -.25*(1.-2.*.3)*t(:,k)*nm(:,k)*nm(:,k) & *(difkgw(:,k-1)+difkgw(:,k))/9.8/prndtlt enddo do k=ktop+2,kbot-1 ttenda(:,k) = -(hflx(:,k-1)-hflx(:,k+1))/ & (zm(:,k-1)-zm(:,k+1))+ & ( (1.-.286)*9.8/287./t(:,k) + & (t(:,k-1)-t(:,k+1))/(zm(:,k-1)-zm(:,k+1)) & /t(:,k) )*hflx(:,k) enddo ttenda(:,ktop+1)=ttenda(:,ktop+2) ttenda(:,kbot)=ttenda(:,kbot-1) do k = ktop+1,kbot ttendb(:,k) = nm(:,k)*nm(:,k)*.25* & (difkgw(:,k-1)+difkgw(:,k))/1004. enddo ttend = (ttenda+ttendb) * 1004. ! convert from K/s to J/kg/s !-time-gcm- return end subroutine gw_tmpr_tend subroutine gw_temp_tend (lchnk ,ncol & ,ngwv ,kbot ,ktop ,dse ,pi ,pm & ,t ,gwut ,ubm ,dt ,rdpm & ,egwdffi ,dttdf ,dttke ,topflx, pcols, pver) ! ! Calculates temperature tendency due to GW breaking ! ! Method: ! GW mixing is calculated using the effective diffusivity obtained from gw_const_tend: ! it is used to diffuse dry static energy. ! ! Author : Sassi - Jan 2001 !--------------------------------------------------------------------------------- ! B.Foster: not available in timegcm: ! use vertical_diffusion, only: vd_lu_decomp, vd_lu_solve, cfluxtop implicit none !---------------------------Input Arguments--------------------------------------- integer, intent(in) :: pcols,pver integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: kbot ! index of bottom (source) midpoint integer, intent(in) :: ktop ! index of top midpoint of gwd region integer, intent(in) :: ngwv ! number of gravity waves to use real(r8), intent(in) :: gwut(pcols,pver,-pgwv:pgwv) ! GW zonal wind tendencies at midpoint real(r8), intent(in) :: ubm(pcols,pver) ! projection of wind at midpoints real(r8), intent(in) :: t(pcols,pver) ! temperature at midpoints real(r8), intent(in) :: dt ! time step real(r8), intent(in) :: rdpm(pcols,pver) ! 1./pdel (thickness bet interfaces) real(r8), intent(in) :: pi(pcols,pver+1) ! interface pressure real(r8), intent(in) :: pm(pcols,pver) ! midpoint pressure real(r8), intent(in) :: egwdffi(pcols,pver+1) ! effective gw diffusivity at interfaces real(r8), intent(in) :: dse(pcols,pver) ! midpoint dry static energy !--------------------------Output Arguments------------------------------------- real(r8), intent(out) :: dttdf(pcols,pver) ! dse tendencies at midpoints: diffusion term real(r8), intent(out) :: dttke(pcols,pver) ! dse tendencies at midpoints: kin. energy conversion term real(r8), intent(out) :: topflx(pcols) ! dse flux at top interface !--------------------------Local Worspace--------------------------------------- integer i,k,m,l real(r8) delta ! dissipation rate real(r8) tmpi2(pcols,pver+1) ! pi(k)/(.5*(tm1(k)+tm1(k-1)) real(r8) ca(pcols,pver) ! -upper diagonal real(r8) cc(pcols,pver) ! -lower diagonal real(r8) term(pcols,pver) ! 1./(1. + ca(k) + cc(k) - cc(k)*ze(k-1)) real(r8) ze(pcols,pver) ! term in tri-diag. matrix system real(r8) temp(pcols,pver) ! temporary storage for constituents real(r8) rho(pcols,pver+1) real(r8) zero(pcols) !------------------------------------------------------------------------------- zero = 0.0 dttdf = 0.0 dttke = 0.0 if (ngwv .eq. 0) return ! ! Compute factor for diffusion matrices at interior interfaces ! Note that [ktop,kbot] are model interfaces (beginning at zero), ! whereas in vd_lu_decomp they are expected as midpoints ! Compute rho at interfaces p/RT, Ti = (Tm_k + Tm_k-1)/2, interface k- do k = 2, pver do i = 1, ncol rho(i,k) = pi(i,k) * 2. / (r * (t(i,k) + t(i,k-1))) end do end do do i = 1, ncol rho(i,pver+1) = pi (i,pver+1) / (r * t(i,pver)) end do do i = 1, ncol rho(i,1) = pi (i,1) / (r * t(i,1)) end do ! Calculate dt * (g*rho)^2/dp at interior interfaces, interface k- tmpi2(:,:)=0.0 do k = ktop+2 , kbot+1 do i = 1, ncol tmpi2(i,k) = dt * (g * rho(i,k))**2 / (pm(i,k) - pm(i,k-1)) end do end do !!$ write (99,*) 'TMPI2:' !!$ do k=1,pver+1 !!$ write (99,*) k,tmpi2(1,k),egwdffi(1,k) !!$ enddo do k = 1, pver do i = 1, ncol temp(i,k) = dse(i,k) enddo enddo ! Constant flux approximation at top interface topflx(:) = 0.0 !!$ call cfluxtop ( & !!$ lchnk ,ncol , & !!$ temp ,egwdffi ,pm ,rdpm ,dt , & !!$ pi ,t ,topflx) ! ! B.Foster 10/03: these routines not available in timegcm1. ! Decompose and solve the diffusion matrix ! call vd_lu_decomp( & ! lchnk ,ncol , & ! zero ,egwdffi ,tmpi2 ,rdpm ,dt ,& ! ca ,cc ,term ,ze ,ktop+1 ,& ! kbot+1 ) ! call vd_lu_solve( & ! lchnk ,ncol , & ! temp ,ca ,ze ,term , & ! ktop+1 ,kbot+1 ) ! Evaluated first part of tendency : diffusive term do k = 1, pver do i =1, ncol dttdf(i,k) = (temp(i,k) - dse(i,k)) / dt enddo enddo ! Evaluate second term: Conversion of kinetic energy into thermal do k = ktop+1, kbot do i = 1, ncol do l = -ngwv, ngwv dttke(i,k) = (c(l)-ubm(i,k)) * gwut(i,k,l) enddo enddo enddo return end subroutine gw_temp_tend !=============================================================================== ! ! B. Foster 10,03: sub lininterp copied from ! CAMROOT/models/atm/cam/src/control/lininterp.F90 ! subroutine lininterp (arrin, yin, nlev, nlatin, arrout, & yout, nlatout) !----------------------------------------------------------------------- ! ! Purpose: Do a linear interpolation from input mesh defined by yin to output ! mesh defined by yout. Where extrapolation is necessary, values will ! be copied from the extreme edge of the input grid. Vectorization is over ! the vertical (nlev) dimension. ! ! Method: Check validity of input, then determine weights, then do the N-S interpolation. ! ! Author: Jim Rosinski !----------------------------------------------------------------------- ! use shr_kind_mod, only: r8 => shr_kind_r8 !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! ! Arguments ! integer, intent(in) :: nlev ! number of vertical levels integer, intent(in) :: nlatin ! number of input latitudes integer, intent(in) :: nlatout ! number of output latitudes real(r8), intent(in) :: arrin(nlev,nlatin) ! input array of values to interpolate real(r8), intent(in) :: yin(nlatin) ! input mesh real(r8), intent(in) :: yout(nlatout) ! output mesh real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array ! ! Local workspace ! integer j, jj ! latitude indices integer js, jn, jjprev ! latitude indices integer k ! level index integer icount ! number of values real(r8) extrap ! percent grid non-overlap ! ! Dynamic ! integer :: jjm(nlatout) integer :: jjp(nlatout) real(r8) :: wgts(nlatout) real(r8) :: wgtn(nlatout) ! ! Check validity of input coordinate arrays: must be monotonically increasing, ! and have a total of at least 2 elements ! if (nlatin.lt.2) then write(6,*)'LININTERP: Must have at least 2 input points for interpolation' ! call endrun call shutdown('lininterp') end if icount = 0 do j=1,nlatin-1 if (yin(j).gt.yin(j+1)) icount = icount + 1 end do do j=1,nlatout-1 if (yout(j).gt.yout(j+1)) icount = icount + 1 end do if (icount.gt.0) then write(6,*)'LININTERP: Non-monotonic coordinate array(s) found' ! call endrun call shutdown('lininterp') end if ! ! Initialize index arrays for later checking ! do j=1,nlatout jjm(j) = 0 jjp(j) = 0 end do ! ! For values which extend beyond N and S boundaries, set weights ! such that values will just be copied. ! do js=1,nlatout if (yout(js).gt.yin(1)) goto 10 jjm(js) = 1 jjp(js) = 1 wgts(js) = 1. wgtn(js) = 0. end do 10 do jn=nlatout,1,-1 if (yout(jn).le.yin(nlatin)) goto 20 jjm(jn) = nlatin jjp(jn) = nlatin wgts(jn) = 1. wgtn(jn) = 0. end do ! ! Loop though output indices finding input indices and weights ! 20 jjprev = 1 do j=js,jn do jj=jjprev,nlatin-1 if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then jjm(j) = jj jjp(j) = jj + 1 wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj)) wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj)) goto 30 end if end do write(6,*)'LININTERP: Failed to find interp values' 30 jjprev = jj end do ! ! Check grid overlap ! extrap = 100.*((js - 1.) + float(nlatout - jn))/nlatout !+waccm+ !!$ if (extrap.gt.30.) then !!$ write(6,*)'LININTERP WARNING:',extrap,' % of output grid will have to be extrapolated' !!$ end if !-waccm- ! ! Check that interp/extrap points have been found for all outputs ! icount = 0 do j=1,nlatout if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1 end do if (icount.gt.0) then write(6,*)'LININTERP: Point found without interp indices' call shutdown('lininterp') end if ! ! Do the interpolation ! do j=1,nlatout do k=1,nlev arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j) end do end do return end subroutine lininterp end module gw_drag