PageRenderTime 63ms CodeModel.GetById 28ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_bl_gwdo.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 700 lines | 454 code | 0 blank | 246 comment | 24 complexity | 43d380e6ed1a2ffc70951ad174f1fc1d MD5 | raw file
Possible License(s): AGPL-1.0
  1. ! WRf:model_layer:physics
  2. !
  3. !
  4. !
  5. !
  6. !
  7. module module_bl_gwdo
  8. contains
  9. !
  10. !-------------------------------------------------------------------
  11. !
  12. subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
  13. rublten,rvblten, &
  14. dtaux3d,dtauy3d,dusfcg,dvsfcg, &
  15. var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, &
  16. znu,znw,mut,p_top, &
  17. cp,g,rd,rv,ep1,pi, &
  18. dt,dx,kpbl2d,itimestep, &
  19. ids,ide, jds,jde, kds,kde, &
  20. ims,ime, jms,jme, kms,kme, &
  21. its,ite, jts,jte, kts,kte)
  22. !-------------------------------------------------------------------
  23. implicit none
  24. !------------------------------------------------------------------------------
  25. !
  26. !-- u3d 3d u-velocity interpolated to theta points (m/s)
  27. !-- v3d 3d v-velocity interpolated to theta points (m/s)
  28. !-- t3d temperature (k)
  29. !-- qv3d 3d water vapor mixing ratio (kg/kg)
  30. !-- p3d 3d pressure (pa)
  31. !-- p3di 3d pressure (pa) at interface level
  32. !-- pi3d 3d exner function (dimensionless)
  33. !-- rublten u tendency due to
  34. ! pbl parameterization (m/s/s)
  35. !-- rvblten v tendency due to
  36. !-- cp heat capacity at constant pressure for dry air (j/kg/k)
  37. !-- g acceleration due to gravity (m/s^2)
  38. !-- rd gas constant for dry air (j/kg/k)
  39. !-- z height above sea level (m)
  40. !-- rv gas constant for water vapor (j/kg/k)
  41. !-- dt time step (s)
  42. !-- dx model grid interval (m)
  43. !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
  44. !-- ids start index for i in domain
  45. !-- ide end index for i in domain
  46. !-- jds start index for j in domain
  47. !-- jde end index for j in domain
  48. !-- kds start index for k in domain
  49. !-- kde end index for k in domain
  50. !-- ims start index for i in memory
  51. !-- ime end index for i in memory
  52. !-- jms start index for j in memory
  53. !-- jme end index for j in memory
  54. !-- kms start index for k in memory
  55. !-- kme end index for k in memory
  56. !-- its start index for i in tile
  57. !-- ite end index for i in tile
  58. !-- jts start index for j in tile
  59. !-- jte end index for j in tile
  60. !-- kts start index for k in tile
  61. !-- kte end index for k in tile
  62. !-------------------------------------------------------------------
  63. !
  64. integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
  65. ims,ime, jms,jme, kms,kme, &
  66. its,ite, jts,jte, kts,kte
  67. integer, intent(in ) :: itimestep
  68. !
  69. real, intent(in ) :: dt,dx,cp,g,rd,rv,ep1,pi
  70. !
  71. real, dimension( ims:ime, kms:kme, jms:jme ) , &
  72. intent(in ) :: qv3d, &
  73. p3d, &
  74. pi3d, &
  75. t3d, &
  76. z
  77. real, dimension( ims:ime, kms:kme, jms:jme ) , &
  78. intent(in ) :: p3di
  79. !
  80. real, dimension( ims:ime, kms:kme, jms:jme ) , &
  81. intent(inout) :: rublten, &
  82. rvblten
  83. real, dimension( ims:ime, kms:kme, jms:jme ) , &
  84. intent(inout) :: dtaux3d, &
  85. dtauy3d
  86. !
  87. real, dimension( ims:ime, kms:kme, jms:jme ) , &
  88. intent(in ) :: u3d, &
  89. v3d
  90. !
  91. integer, dimension( ims:ime, jms:jme ) , &
  92. intent(in ) :: kpbl2d
  93. real, dimension( ims:ime, jms:jme ) , &
  94. intent(inout ) :: dusfcg, &
  95. dvsfcg
  96. !
  97. real, dimension( ims:ime, jms:jme ) , &
  98. intent(in ) :: var2d, &
  99. oc12d, &
  100. oa2d1,oa2d2,oa2d3,oa2d4, &
  101. ol2d1,ol2d2,ol2d3,ol2d4
  102. !
  103. real, dimension( ims:ime, jms:jme ) , &
  104. optional , &
  105. intent(in ) :: mut
  106. !
  107. real, dimension( kms:kme ) , &
  108. optional , &
  109. intent(in ) :: znu, &
  110. znw
  111. !
  112. real, optional, intent(in ) :: p_top
  113. !
  114. !local
  115. !
  116. real, dimension( its:ite, kts:kte ) :: delprsi, &
  117. pdh
  118. real, dimension( its:ite, kts:kte+1 ) :: pdhi
  119. real, dimension( its:ite, 4 ) :: oa4, &
  120. ol4
  121. integer :: i,j,k,kdt
  122. !
  123. do j = jts,jte
  124. if(present(mut))then
  125. ! For ARW we will replace p and p8w with dry hydrostatic pressure
  126. do k = kts,kte+1
  127. do i = its,ite
  128. if(k.le.kte)pdh(i,k) = mut(i,j)*znu(k) + p_top
  129. pdhi(i,k) = mut(i,j)*znw(k) + p_top
  130. enddo
  131. enddo
  132. else
  133. do k = kts,kte+1
  134. do i = its,ite
  135. if(k.le.kte)pdh(i,k) = p3d(i,k,j)
  136. pdhi(i,k) = p3di(i,k,j)
  137. enddo
  138. enddo
  139. endif
  140. !
  141. do k = kts,kte
  142. do i = its,ite
  143. delprsi(i,k) = pdhi(i,k)-pdhi(i,k+1)
  144. enddo
  145. enddo
  146. do i = its,ite
  147. oa4(i,1) = oa2d1(i,j)
  148. oa4(i,2) = oa2d2(i,j)
  149. oa4(i,3) = oa2d3(i,j)
  150. oa4(i,4) = oa2d4(i,j)
  151. ol4(i,1) = ol2d1(i,j)
  152. ol4(i,2) = ol2d2(i,j)
  153. ol4(i,3) = ol2d3(i,j)
  154. ol4(i,4) = ol2d4(i,j)
  155. enddo
  156. call gwdo2d(dudt=rublten(ims,kms,j),dvdt=rvblten(ims,kms,j) &
  157. ,dtaux2d=dtaux3d(ims,kms,j),dtauy2d=dtauy3d(ims,kms,j) &
  158. ,u1=u3d(ims,kms,j),v1=v3d(ims,kms,j) &
  159. ,t1=t3d(ims,kms,j),q1=qv3d(ims,kms,j) &
  160. ,prsi=pdhi(its,kts),del=delprsi(its,kts) &
  161. ,prsl=pdh(its,kts),prslk=pi3d(ims,kms,j) &
  162. ,zl=z(ims,kms,j),rcl=1.0 &
  163. ,dusfc=dusfcg(ims,j),dvsfc=dvsfcg(ims,j) &
  164. ,var=var2d(ims,j),oc1=oc12d(ims,j) &
  165. ,oa4=oa4,ol4=ol4 &
  166. ,g=g,cp=cp,rd=rd,rv=rv,fv=ep1,pi=pi &
  167. ,dxmeter=dx,deltim=dt &
  168. ,kpbl=kpbl2d(ims,j),kdt=itimestep,lat=j &
  169. ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
  170. ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
  171. ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
  172. enddo
  173. !
  174. !
  175. end subroutine gwdo
  176. !
  177. !-------------------------------------------------------------------
  178. !
  179. !
  180. !
  181. !
  182. subroutine gwdo2d(dudt,dvdt,dtaux2d,dtauy2d, &
  183. u1,v1,t1,q1, &
  184. prsi,del,prsl,prslk,zl,rcl, &
  185. var,oc1,oa4,ol4,dusfc,dvsfc, &
  186. g,cp,rd,rv,fv,pi,dxmeter,deltim,kpbl,kdt,lat, &
  187. ids,ide, jds,jde, kds,kde, &
  188. ims,ime, jms,jme, kms,kme, &
  189. its,ite, jts,jte, kts,kte)
  190. !-------------------------------------------------------------------
  191. !
  192. ! this code handles the time tendencies of u v due to the effect of mountain
  193. ! induced gravity wave drag from sub-grid scale orography. this routine
  194. ! not only treats the traditional upper-level wave breaking due to mountain
  195. ! variance (alpert 1988), but also the enhanced lower-tropospheric wave
  196. ! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
  197. ! thus, in addition to the terrain height data in a model grid gox,
  198. ! additional 10-2d topographic statistics files are needed, including
  199. ! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
  200. ! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
  201. ! hong (1999). the current scheme was implmented as in hong et al.(2008)
  202. !
  203. ! coded by song-you hong and young-joon kim and implemented by song-you hong
  204. !
  205. ! references:
  206. ! hong et al. (2008), wea. and forecasting
  207. ! kim and arakawa (1995), j. atmos. sci.
  208. ! alpet et al. (1988), NWP conference.
  209. ! hong (1999), NCEP office note 424.
  210. !
  211. ! notice : comparible or lower resolution orography files than model resolution
  212. ! are desirable in preprocess (wps) to prevent weakening of the drag
  213. !-------------------------------------------------------------------
  214. !
  215. ! input
  216. ! dudt (ims:ime,kms:kme) non-lin tendency for u wind component
  217. ! dvdt (ims:ime,kms:kme) non-lin tendency for v wind component
  218. ! u1(ims:ime,kms:kme) zonal wind / sqrt(rcl) m/sec at t0-dt
  219. ! v1(ims:ime,kms:kme) meridional wind / sqrt(rcl) m/sec at t0-dt
  220. ! t1(ims:ime,kms:kme) temperature deg k at t0-dt
  221. ! q1(ims:ime,kms:kme) specific humidity at t0-dt
  222. !
  223. ! rcl a scaling factor = reciprocal of square of cos(lat)
  224. ! for mrf gsm. rcl=1 if u1 and v1 are wind components.
  225. ! deltim time step secs
  226. ! del(kts:kte) positive increment of pressure across layer (pa)
  227. !
  228. ! output
  229. ! dudt, dvdt wind tendency due to gwdo
  230. !
  231. !-------------------------------------------------------------------
  232. implicit none
  233. !-------------------------------------------------------------------
  234. integer :: kdt,lat,latd,lond, &
  235. ids,ide, jds,jde, kds,kde, &
  236. ims,ime, jms,jme, kms,kme, &
  237. its,ite, jts,jte, kts,kte
  238. !
  239. real :: g,rd,rv,fv,cp,pi,dxmeter,deltim,rcl
  240. real :: dudt(ims:ime,kms:kme),dvdt(ims:ime,kms:kme), &
  241. dtaux2d(ims:ime,kms:kme),dtauy2d(ims:ime,kms:kme), &
  242. u1(ims:ime,kms:kme),v1(ims:ime,kms:kme), &
  243. t1(ims:ime,kms:kme),q1(ims:ime,kms:kme), &
  244. zl(ims:ime,kms:kme),prslk(ims:ime,kms:kme)
  245. real :: prsl(its:ite,kts:kte),prsi(its:ite,kts:kte+1), &
  246. del(its:ite,kts:kte)
  247. real :: oa4(its:ite,4),ol4(its:ite,4)
  248. !
  249. integer :: kpbl(ims:ime)
  250. real :: var(ims:ime),oc1(ims:ime), &
  251. dusfc(ims:ime),dvsfc(ims:ime)
  252. ! critical richardson number for wave breaking : ! larger drag with larger value
  253. !
  254. real,parameter :: ric = 0.25
  255. !
  256. real,parameter :: dw2min = 1.
  257. real,parameter :: rimin = -100.
  258. real,parameter :: bnv2min = 1.0e-5
  259. real,parameter :: efmin = 0.0
  260. real,parameter :: efmax = 10.0
  261. real,parameter :: xl = 4.0e4
  262. real,parameter :: critac = 1.0e-5
  263. real,parameter :: gmax = 1.
  264. real,parameter :: veleps = 1.0
  265. real,parameter :: factop = 0.5
  266. real,parameter :: frc = 1.0
  267. real,parameter :: ce = 0.8
  268. real,parameter :: cg = 0.5
  269. !
  270. ! local variables
  271. !
  272. integer :: i,k,lcap,lcapp1,nwd,idir,kpblmin,kpblmax, &
  273. klcap,kp1,ikount,kk
  274. !
  275. real :: rcs,rclcs,csg,fdir,cleff,cs,rcsks, &
  276. wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, &
  277. wtkbj,coefm,tem,gfobnv,hd,fro,rim,temc,tem1,efact, &
  278. temv,dtaux,dtauy
  279. !
  280. logical :: ldrag(its:ite),icrilv(its:ite), &
  281. flag(its:ite),kloop1(its:ite)
  282. !
  283. real :: taub(its:ite),taup(its:ite,kts:kte+1), &
  284. xn(its:ite),yn(its:ite), &
  285. ubar(its:ite),vbar(its:ite), &
  286. fr(its:ite),ulow(its:ite), &
  287. rulow(its:ite),bnv(its:ite), &
  288. oa(its:ite),ol(its:ite), &
  289. roll(its:ite),dtfac(its:ite), &
  290. brvf(its:ite),xlinv(its:ite), &
  291. delks(its:ite),delks1(its:ite), &
  292. bnv2(its:ite,kts:kte),usqj(its:ite,kts:kte), &
  293. taud(its:ite,kts:kte),ro(its:ite,kts:kte), &
  294. vtk(its:ite,kts:kte),vtj(its:ite,kts:kte), &
  295. zlowtop(its:ite),velco(its:ite,kts:kte-1)
  296. !
  297. integer :: kbl(its:ite),klowtop(its:ite), &
  298. lowlv(its:ite)
  299. !
  300. logical :: iope
  301. integer,parameter :: mdir=8
  302. integer :: nwdir(mdir)
  303. data nwdir/6,7,5,8,2,3,1,4/
  304. !
  305. ! initialize local variables
  306. !
  307. kbl=0 ; klowtop=0 ; lowlv=0
  308. !
  309. !---- constants
  310. !
  311. rcs = sqrt(rcl)
  312. cs = 1. / sqrt(rcl)
  313. csg = cs * g
  314. lcap = kte
  315. lcapp1 = lcap + 1
  316. fdir = mdir / (2.0*pi)
  317. !
  318. !
  319. !!!!!!! cleff (subgrid mountain scale ) is highly tunable parameter
  320. !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
  321. !
  322. cleff = max(dxmeter,50.e3)
  323. !
  324. ! initialize!!
  325. !
  326. dtaux = 0.0
  327. dtauy = 0.0
  328. do k = kts,kte
  329. do i = its,ite
  330. usqj(i,k) = 0.0
  331. bnv2(i,k) = 0.0
  332. vtj(i,k) = 0.0
  333. vtk(i,k) = 0.0
  334. taup(i,k) = 0.0
  335. taud(i,k) = 0.0
  336. dtaux2d(i,k)= 0.0
  337. dtauy2d(i,k)= 0.0
  338. enddo
  339. enddo
  340. do i = its,ite
  341. taup(i,kte+1) = 0.0
  342. xlinv(i) = 1.0/xl
  343. enddo
  344. !
  345. do k = kts,kte
  346. do i = its,ite
  347. vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
  348. vtk(i,k) = vtj(i,k) / prslk(i,k)
  349. ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3
  350. enddo
  351. enddo
  352. !
  353. do i = its,ite
  354. zlowtop(i) = 2. * var(i)
  355. enddo
  356. !
  357. !--- determine new reference level > 2*var
  358. !
  359. do i = its,ite
  360. kloop1(i) = .true.
  361. enddo
  362. do k = kts+1,kte
  363. do i = its,ite
  364. if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
  365. klowtop(i) = k+1
  366. kloop1(i) = .false.
  367. endif
  368. enddo
  369. enddo
  370. !
  371. kpblmax = 2
  372. do i = its,ite
  373. kbl(i) = max(2, kpbl(i))
  374. kbl(i) = max(kbl(i), klowtop(i))
  375. delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
  376. ubar (i) = 0.0
  377. vbar (i) = 0.0
  378. taup(i,1) = 0.0
  379. oa(i) = 0.0
  380. kpblmax = max(kpblmax,kbl(i))
  381. flag(i) = .true.
  382. lowlv(i) = 2
  383. enddo
  384. kpblmax = min(kpblmax+1,kte-1)
  385. !
  386. ! compute low level averages within pbl
  387. !
  388. do k = kts,kpblmax
  389. do i = its,ite
  390. if (k.lt.kbl(i)) then
  391. rcsks = rcs * del(i,k) * delks(i)
  392. ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
  393. vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
  394. endif
  395. enddo
  396. enddo
  397. !
  398. ! figure out low-level horizontal wind direction
  399. !
  400. ! nwd 1 2 3 4 5 6 7 8
  401. ! wd w s sw nw e n ne se
  402. !
  403. do i = its,ite
  404. wdir = atan2(ubar(i),vbar(i)) + pi
  405. idir = mod(nint(fdir*wdir),mdir) + 1
  406. nwd = nwdir(idir)
  407. oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
  408. ol(i) = ol4(i,mod(nwd-1,4)+1)
  409. enddo
  410. !
  411. kpblmin = kte
  412. do i = its,ite
  413. kpblmin = min(kpblmin, kbl(i))
  414. enddo
  415. !
  416. do i = its,ite
  417. if (oa(i).le.0.0) kbl(i) = kpbl(i) + 1
  418. enddo
  419. !
  420. do i = its,ite
  421. delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
  422. delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
  423. enddo
  424. !
  425. !--- saving richardson number in usqj for migwdi
  426. !
  427. do k = kts,kte-1
  428. do i = its,ite
  429. ti = 2.0 / (t1(i,k)+t1(i,k+1))
  430. rdz = 1./(zl(i,k+1) - zl(i,k))
  431. tem1 = u1(i,k) - u1(i,k+1)
  432. tem2 = v1(i,k) - v1(i,k+1)
  433. dw2 = rcl*(tem1*tem1 + tem2*tem2)
  434. shr2 = max(dw2,dw2min) * rdz * rdz
  435. bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
  436. usqj(i,k) = max(bvf2/shr2,rimin)
  437. bnv2(i,k) = 2*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
  438. bnv2(i,k) = max( bnv2(i,k), bnv2min )
  439. enddo
  440. enddo
  441. !
  442. !-----initialize arrays
  443. !
  444. do i = its,ite
  445. xn(i) = 0.0
  446. yn(i) = 0.0
  447. ubar (i) = 0.0
  448. vbar (i) = 0.0
  449. roll (i) = 0.0
  450. taub (i) = 0.0
  451. ulow (i) = 0.0
  452. dtfac(i) = 1.0
  453. ldrag(i) = .false.
  454. icrilv(i) = .false. ! initialize critical level control vector
  455. enddo
  456. !
  457. !---- compute low level averages
  458. !---- (u,v)*cos(lat) use uv=(u1,v1) which is wind at t0-1
  459. !---- use rcs=1/cos(lat) to get wind field
  460. !
  461. do k = 1,kpblmax
  462. do i = its,ite
  463. if (k .lt. kbl(i)) then
  464. rdelks = del(i,k) * delks(i)
  465. rcsks = rcs * rdelks
  466. ubar(i) = ubar(i) + rcsks * u1(i,k) ! u mean
  467. vbar(i) = vbar(i) + rcsks * v1(i,k) ! v mean
  468. roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
  469. endif
  470. enddo
  471. enddo
  472. !
  473. !----compute the "low level" or 1/3 wind magnitude (m/s)
  474. !
  475. do i = its,ite
  476. ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
  477. rulow(i) = 1./ulow(i)
  478. enddo
  479. !
  480. do k = kts,kte-1
  481. do i = its,ite
  482. velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
  483. + (v1(i,k)+v1(i,k+1)) * vbar(i))
  484. velco(i,k) = velco(i,k) * rulow(i)
  485. if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
  486. velco(i,k) = veleps
  487. endif
  488. enddo
  489. enddo
  490. !
  491. ! no drag when critical level in the base layer
  492. !
  493. do i = its,ite
  494. ldrag(i) = velco(i,1).le.0.
  495. enddo
  496. !
  497. do k = kts+1,kpblmax-1
  498. do i = its,ite
  499. if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
  500. enddo
  501. enddo
  502. !
  503. ! no drag when bnv2.lt.0
  504. !
  505. do k = kts,kpblmax-1
  506. do i = its,ite
  507. if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
  508. enddo
  509. enddo
  510. !
  511. !-----the low level weighted average ri is stored in usqj(1,1; im)
  512. !-----the low level weighted average n**2 is stored in bnv2(1,1; im)
  513. !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
  514. !---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
  515. !
  516. do i = its,ite
  517. wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
  518. bnv2(i,1) = wtkbj * bnv2(i,1)
  519. usqj(i,1) = wtkbj * usqj(i,1)
  520. enddo
  521. !
  522. do k = kts+1,kpblmax-1
  523. do i = its,ite
  524. if (k .lt. kbl(i)) then
  525. rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
  526. bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
  527. usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
  528. endif
  529. enddo
  530. enddo
  531. !
  532. do i = its,ite
  533. ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
  534. ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
  535. ldrag(i) = ldrag(i) .or. var(i) .le. 0.0
  536. enddo
  537. !
  538. ! ----- set all ri low level values to the low level value
  539. !
  540. do k = kts+1,kpblmax-1
  541. do i = its,ite
  542. if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
  543. enddo
  544. enddo
  545. !
  546. do i = its,ite
  547. if (.not.ldrag(i)) then
  548. bnv(i) = sqrt( bnv2(i,1) )
  549. fr(i) = bnv(i) * rulow(i) * var(i)
  550. xn(i) = ubar(i) * rulow(i)
  551. yn(i) = vbar(i) * rulow(i)
  552. endif
  553. enddo
  554. !
  555. ! compute the base level stress and store it in taub
  556. ! calculate enhancement factor, number of mountains & aspect
  557. ! ratio const. use simplified relationship between standard
  558. ! deviation & critical hgt
  559. !
  560. do i = its,ite
  561. if (.not. ldrag(i)) then
  562. efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
  563. efact = min( max(efact,efmin), efmax )
  564. coefm = (1. + ol(i)) ** (oa(i)+1.)
  565. xlinv(i) = coefm / cleff
  566. tem = fr(i) * fr(i) * oc1(i)
  567. gfobnv = gmax * tem / ((tem + cg)*bnv(i))
  568. taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
  569. * ulow(i) * gfobnv * efact
  570. else
  571. taub(i) = 0.0
  572. xn(i) = 0.0
  573. yn(i) = 0.0
  574. endif
  575. enddo
  576. !
  577. ! now compute vertical structure of the stress.
  578. !
  579. !----set up bottom values of stress
  580. !
  581. do k = kts,kpblmax
  582. do i = its,ite
  583. if (k .le. kbl(i)) taup(i,k) = taub(i)
  584. enddo
  585. enddo
  586. !
  587. do k = kpblmin, kte-1 ! vertical level k loop!
  588. kp1 = k + 1
  589. do i = its,ite
  590. !
  591. !-----unstablelayer if ri < ric
  592. !-----unstable layer if upper air vel comp along surf vel <=0 (crit lay)
  593. !---- at (u-c)=0. crit layer exists and bit vector should be set (.le.)
  594. !
  595. if (k .ge. kbl(i)) then
  596. icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
  597. .or. (velco(i,k) .le. 0.0)
  598. brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
  599. brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
  600. endif
  601. enddo
  602. !
  603. do i = its,ite
  604. if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
  605. if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
  606. temv = 1.0 / velco(i,k)
  607. tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5
  608. hd = sqrt(taup(i,k) / tem1)
  609. fro = brvf(i) * hd * temv
  610. !
  611. ! rim is the minimum-richardson number by shutts (1985)
  612. !
  613. tem2 = sqrt(usqj(i,k))
  614. tem = 1. + tem2 * fro
  615. rim = usqj(i,k) * (1.-fro) / (tem * tem)
  616. !
  617. ! check stability to employ the 'saturation hypothesis'
  618. ! of lindzen (1981) except at tropospheric downstream regions
  619. !
  620. if (rim .le. ric) then ! saturation hypothesis!
  621. if ((oa(i) .le. 0. .or. kp1 .ge. lowlv(i) )) then
  622. temc = 2.0 + 1.0 / tem2
  623. hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
  624. taup(i,kp1) = tem1 * hd * hd
  625. endif
  626. else ! no wavebreaking!
  627. taup(i,kp1) = taup(i,k)
  628. endif
  629. endif
  630. endif
  631. enddo
  632. enddo
  633. !
  634. if(lcap.lt.kte) then
  635. do klcap = lcapp1,kte
  636. do i = its,ite
  637. taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
  638. enddo
  639. enddo
  640. endif
  641. !
  642. ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
  643. !
  644. do k = kts,kte
  645. do i = its,ite
  646. taud(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
  647. enddo
  648. enddo
  649. !
  650. !------limit de-acceleration (momentum deposition ) at top to 1/2 value
  651. !------the idea is some stuff must go out the 'top'
  652. !
  653. do klcap = lcap,kte
  654. do i = its,ite
  655. taud(i,klcap) = taud(i,klcap) * factop
  656. enddo
  657. enddo
  658. !
  659. !------if the gravity wave drag would force a critical line
  660. !------in the lower ksmm1 layers during the next deltim timestep,
  661. !------then only apply drag until that critical line is reached.
  662. !
  663. do k = kts,kpblmax-1
  664. do i = its,ite
  665. if (k .le. kbl(i)) then
  666. if(taud(i,k).ne.0.) &
  667. dtfac(i) = min(dtfac(i),abs(velco(i,k) &
  668. /(deltim*rcs*taud(i,k))))
  669. endif
  670. enddo
  671. enddo
  672. !
  673. do i = its,ite
  674. dusfc(i) = 0.
  675. dvsfc(i) = 0.
  676. enddo
  677. !
  678. do k = kts,kte
  679. do i = its,ite
  680. taud(i,k) = taud(i,k) * dtfac(i)
  681. dtaux = taud(i,k) * xn(i)
  682. dtauy = taud(i,k) * yn(i)
  683. dtaux2d(i,k) = dtaux
  684. dtauy2d(i,k) = dtauy
  685. dudt(i,k) = dtaux + dudt(i,k)
  686. dvdt(i,k) = dtauy + dvdt(i,k)
  687. dusfc(i) = dusfc(i) + dtaux * del(i,k)
  688. dvsfc(i) = dvsfc(i) + dtauy * del(i,k)
  689. enddo
  690. enddo
  691. !
  692. do i = its,ite
  693. dusfc(i) = (-1./g*rcs) * dusfc(i)
  694. dvsfc(i) = (-1./g*rcs) * dvsfc(i)
  695. enddo
  696. !
  697. return
  698. end subroutine gwdo2d
  699. !-------------------------------------------------------------------
  700. end module module_bl_gwdo