CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
drag_suite.F90
1
4
6 module drag_suite
7
8 contains
9
15 subroutine drag_suite_init(gwd_opt, errmsg, errflg)
16
17 integer, intent(in) :: gwd_opt
18 character(len=*), intent(out) :: errmsg
19 integer, intent(out) :: errflg
20
21
22 ! Initialize CCPP error handling variables
23 errmsg = ''
24 errflg = 0
25
26 ! Consistency checks
27 if (gwd_opt/=3 .and. gwd_opt/=33) then
28 write(errmsg,'(*(a))') "Logic error: namelist choice of gravity wave &
29 & drag is different from drag_suite scheme"
30 errflg = 1
31 return
32 end if
33 end subroutine drag_suite_init
34
207 subroutine drag_suite_run( &
208 & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, &
209 & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, &
210 & var,oc1,oa4,ol4, &
211 & varss,oc1ss,oa4ss,ol4ss, &
212 & THETA,SIGMA,GAMMA,ELVMAX, &
213 & dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,dtauy2d_bl, &
214 & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, &
215 & dusfc,dvsfc, &
216 & dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl, &
217 & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
218 & slmsk,br1,hpbl, &
219 & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
220 & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
221 & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
222 & dtend, dtidx, index_of_process_orographic_gwd, &
223 & index_of_temperature, index_of_x_wind, &
224 & index_of_y_wind, ldiag3d, ldiag_ugwp, ugwp_seq_update, &
225 & spp_wts_gwd, spp_gwd, errmsg, errflg)
226
227! ********************************************************************
228! -----> I M P L E M E N T A T I O N V E R S I O N <----------
229!
230! ----- This code -----
231!begin WRF code
232
233! this code handles the time tendencies of u v due to the effect of mountain
234! induced gravity wave drag from sub-grid scale orography. this routine
235! not only treats the traditional upper-level wave breaking due to mountain
236! variance (alpert 1988), but also the enhanced lower-tropospheric wave
237! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
238! thus, in addition to the terrain height data in a model grid box,
239! additional 10-2d topographic statistics files are needed, including
240! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
241! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
242! hong (1999). the current scheme was implmented as in hong et al.(2008)
243!
244! Originally coded by song-you hong and young-joon kim and implemented by song-you hong
245!
246! program history log:
247! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle
248! with blocked height by dividing streamline theory
249! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale
250! orographic grabity wave drag:
251! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the
252! topographic form drag of Beljaars et al. (2004, QJRMS)
253! Activation of each component is done by specifying the integer-parameters
254! (defined below) to 0: inactive or 1: active
255! gwd_opt_ms = 0 or 1: mesoscale (changed to logical flag)
256! gwd_opt_bl = 0 or 1: blocking drag (changed to logical flag)
257! gwd_opt_ss = 0 or 1: small-scale gravity wave drag (removed)
258! gwd_opt_fd = 0 or 1: topographic form drag (removed)
259! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating (logical flags)
260! gsd_diss_ht_opt = .false. : dissipation heating off
261! gsd_diss_ht_opt = .true. : dissipation heating on
262! 2020-08-25 Michael Toy changed logic control for drag component selection
263! for CCPP.
264! Namelist options:
265! do_gsl_drag_ls_bl - logical flag for mesoscale GWD + blocking
266! do_gsl_drag_ss - logical flag for small-scale GWD
267! do_gsl_drag_tofd - logical flag for turbulent form drag
268! Compile-time options (changed from integer switches to logical flags):
269! gwd_opt_ms = : mesoscale GWD (active if == .true.)
270! gwd_opt_bl = : blocking drag (active if == .true.)
271!
272! References:
273! Hong et al. (2008), wea. and forecasting
274! Kim and Doyle (2005), Q. J. R. Meteor. Soc.
275! Kim and Arakawa (1995), j. atmos. sci.
276! Alpert et al. (1988), NWP conference.
277! Hong (1999), NCEP office note 424.
278! Steeneveld et al (2008), JAMC
279! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc.
280! Beljaars et al. (2004), Q. J. R. Meteor. Soc.
281!
282! notice : comparible or lower resolution orography files than model resolution
283! are desirable in preprocess (wps) to prevent weakening of the drag
284!-------------------------------------------------------------------------------
285!
286! input
287! dudt (im,km) non-lin tendency for u wind component
288! dvdt (im,km) non-lin tendency for v wind component
289! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt
290! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt
291! t1(im,km) temperature deg k at t0-dt
292! q1(im,km) specific humidity at t0-dt
293! deltim time step secs
294! del(km) positive increment of pressure across layer (pa)
295! KPBL(IM) is the index of the top layer of the PBL
296! ipr & lprnt for diagnostics
297!
298! output
299! dudt, dvdt wind tendency due to gwdo
300! dTdt
301!
302!-------------------------------------------------------------------------------
303
304!end wrf code
305!----------------------------------------------------------------------C
306! USE
307! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES)
308!
309! PURPOSE
310! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
311! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V
312! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
313! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
314! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
315! CRITICAL LEVELS
316!
317!
318! ********************************************************************
319 USE machine , ONLY : kind_phys
320 implicit none
321
322 ! Interface variables
323 integer, intent(in) :: im, km, imx, kdt, ipr, me, master
324 integer, intent(in) :: gwd_opt
325 logical, intent(in) :: lprnt
326 integer, intent(in) :: KPBL(:)
327 real(kind=kind_phys), intent(in) :: deltim, g, cp, rd, rv, &
328 & cdmbgwd(:), alpha_fd
329 real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
330 logical, intent(in) :: ldiag3d
331 integer, intent(in) :: dtidx(:,:), index_of_temperature, &
332 & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind
333
334 integer :: kpblmax
335 integer, parameter :: ims=1, kms=1, its=1, kts=1
336 real(kind=kind_phys), intent(in) :: fv, pi
337 real(kind=kind_phys) :: rcl, cdmb
338 real(kind=kind_phys) :: g_inv, rd_inv
339
340 real(kind=kind_phys), intent(inout) :: &
341 & dudt(:,:),dvdt(:,:), &
342 & dtdt(:,:)
343 real(kind=kind_phys), intent(inout) :: rdxzb(:)
344 real(kind=kind_phys), intent(in) :: &
345 & u1(:,:),v1(:,:), &
346 & t1(:,:),q1(:,:), &
347 & phii(:,:),prsl(:,:), &
348 & prslk(:,:),phil(:,:)
349 real(kind=kind_phys), intent(in) :: prsi(:,:), &
350 & del(:,:)
351 real(kind=kind_phys), intent(in) :: var(:),oc1(:), &
352 & oa4(:,:),ol4(:,:), &
353 & dx(:)
354 real(kind=kind_phys), intent(in), optional :: varss(:),oc1ss(:), &
355 & oa4ss(:,:),ol4ss(:,:)
356 real(kind=kind_phys), intent(in) :: theta(:),sigma(:), &
357 & gamma(:),elvmax(:)
358
359! added for small-scale orographic wave drag
360 real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx
361 real(kind=kind_phys), intent(in) :: br1(:), &
362 & hpbl(:), &
363 & slmsk(:)
364 real(kind=kind_phys), dimension(im) :: govrth,xland
365 !real(kind=kind_phys), dimension(im,km) :: dz2
366 real(kind=kind_phys) :: tauwavex0,tauwavey0, &
367 & xnbv,density,tvcon,hpbl2
368 integer :: kpbl2,kvar
369 !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g
370 real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g
371
372!SPP
373 real(kind=kind_phys), dimension(im) :: var_stoch, varss_stoch, &
374 varmax_fd_stoch
375 real(kind=kind_phys), intent(in), optional :: spp_wts_gwd(:,:)
376 integer, intent(in) :: spp_gwd
377
378 real(kind=kind_phys), dimension(im) :: rstoch
379
380!Output:
381 real(kind=kind_phys), intent(inout) :: &
382 & dusfc(:), dvsfc(:)
383!Output (optional):
384 real(kind=kind_phys), intent(inout), optional :: &
385 & dusfc_ms(:),dvsfc_ms(:), &
386 & dusfc_bl(:),dvsfc_bl(:), &
387 & dusfc_ss(:),dvsfc_ss(:), &
388 & dusfc_fd(:),dvsfc_fd(:)
389 real(kind=kind_phys), intent(inout), optional :: &
390 & dtaux2d_ms(:,:),dtauy2d_ms(:,:), &
391 & dtaux2d_bl(:,:),dtauy2d_bl(:,:), &
392 & dtaux2d_ss(:,:),dtauy2d_ss(:,:), &
393 & dtaux2d_fd(:,:),dtauy2d_fd(:,:)
394
395!Misc arrays
396 real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d
397
398!-------------------------------------------------------------------------
399! Flags to regulate the activation of specific components of drag suite:
400! Each component is tapered off automatically as a function of dx, so best to
401! keep them activated (.true.).
402 logical, intent(in) :: &
403 do_gsl_drag_ls_bl, & ! mesoscale gravity wave drag and blocking
404 do_gsl_drag_ss, & ! small-scale gravity wave drag (Steeneveld et al. 2008)
405 do_gsl_drag_tofd ! form drag (Beljaars et al. 2004, QJRMS)
406
407! Flag for diagnostic outputs
408 logical, intent(in) :: ldiag_ugwp
409
410! Flag for sequential update of u and v between
411! LSGWD + BLOCKING and SSGWD + TOFD calculations
412 logical, intent(in) :: ugwp_seq_update
413
414! More variables for sequential updating of winds
415 ! Updated winds
416 real(kind=kind_phys), dimension(im,km) :: uwnd1, vwnd1
417 real(kind=kind_phys) :: tmp1, tmp2 ! temporary variables
418
419! Additional flags
420 logical, parameter :: &
421 gwd_opt_ms = .true., & ! mesoscale gravity wave drag
422 gwd_opt_bl = .true., & ! blocking drag
423 gsd_diss_ht_opt = .false. ! dissipation heating
424
425! Parameters for bounding the scale-adaptive variability:
426! Small-scale GWD + turbulent form drag
427 real(kind=kind_phys), parameter :: dxmin_ss = 1000., &
428 & dxmax_ss = 12000. ! min,max range of tapering (m)
429! Mesoscale GWD + blocking
430 real(kind=kind_phys), parameter :: dxmin_ms = 3000., &
431 & dxmax_ms = 13000. ! min,max range of tapering (m)
432 real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and meso-scale tapering factors (-)
433!
434! Variables for limiting topographic standard deviation (var)
435 real(kind=kind_phys), parameter :: varmax_ss = 50., & ! varmax_ss not used
436 varmax_fd = 500., &
437 beta_fd = 0.2
438 real(kind=kind_phys) :: var_temp, var_temp2
439
440! added Beljaars orographic form drag
441 real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
442 real(kind=kind_phys) :: a1,a2,wsp
443 real(kind=kind_phys) :: h_efold
444 real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3
445
446! critical richardson number for wave breaking : ! larger drag with larger value
447 real(kind=kind_phys), parameter :: ric = 0.25
448 real(kind=kind_phys), parameter :: dw2min = 1.
449 real(kind=kind_phys), parameter :: rimin = -100.
450 real(kind=kind_phys), parameter :: bnv2min = 1.0e-5
451 real(kind=kind_phys), parameter :: efmin = 0.0
452 real(kind=kind_phys), parameter :: efmax = 10.0
453 real(kind=kind_phys), parameter :: xl = 4.0e4
454 real(kind=kind_phys), parameter :: critac = 1.0e-5
455 real(kind=kind_phys), parameter :: gmax = 1.
456 real(kind=kind_phys), parameter :: veleps = 1.0
457 real(kind=kind_phys), parameter :: factop = 0.5
458 real(kind=kind_phys), parameter :: frc = 1.0
459 real(kind=kind_phys), parameter :: ce = 0.8
460 real(kind=kind_phys), parameter :: cg = 0.5
461 real(kind=kind_phys), parameter :: sgmalolev = 0.5 ! max sigma lvl for dtfac
462 real(kind=kind_phys), parameter :: plolevmeso = 70.0 ! pres lvl for mesosphere OGWD reduction (Pa)
463 real(kind=kind_phys), parameter :: facmeso = 0.5 ! fractional velocity reduction for OGWD
464 integer,parameter :: kpblmin = 2
465
466!
467! local variables
468!
469 integer :: i,j,k,lcap,lcapp1,nwd,idir, &
470 klcap,kp1
471!
472 real(kind=kind_phys) :: rcs,csg,fdir,cleff,cleff_ss,cs, &
473 rcsks,wdir,ti,rdz,tem2,dw2,shr2, &
474 bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, &
475 rim,temc,tem1,efact,temv,dtaux,dtauy, &
476 dtauxb,dtauyb,eng0,eng1,ksmax,dtfac_meso
477!
478 logical :: ldrag(im),icrilv(im), &
479 flag(im),kloop1(im)
480 logical :: prop_test
481!
482 real(kind=kind_phys) :: onebgrcs
483!
484 real(kind=kind_phys) :: taub(im),taup(im,km+1), &
485 xn(im),yn(im), &
486 ubar(im),vbar(im), &
487 fr(im),ulow(im), &
488 rulow(im),bnv(im), &
489 oa(im),ol(im), &
490 oass(im),olss(im), &
491 roll(im),dtfac(im), &
492 brvf(im),xlinv(im), &
493 delks(im),delks1(im), &
494 bnv2(im,km),usqj(im,km), &
495 taud_ms(im,km),taud_bl(im,km), &
496 ro(im,km), &
497 vtk(im,km),vtj(im,km), &
498 zlowtop(im),velco(im,km-1), &
499 coefm(im),coefm_ss(im)
500!
501 integer :: kbl(im),klowtop(im)
502 integer,parameter :: mdir=8
503 !integer :: nwdir(mdir)
504 !data nwdir/6,7,5,8,2,3,1,4/
505 integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/)
506!
507! variables for flow-blocking drag
508!
509 real(kind=kind_phys),parameter :: frmax = 10.
510 real(kind=kind_phys),parameter :: olmin = 1.0e-5
511 real(kind=kind_phys),parameter :: odmin = 0.1
512 real(kind=kind_phys),parameter :: odmax = 10.
513 integer :: komax(im)
514 integer :: kblk
515 real(kind=kind_phys) :: cd
516 real(kind=kind_phys) :: zblk,tautem
517 real(kind=kind_phys) :: pe,ke
518 real(kind=kind_phys) :: delx,dely
519 real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4)
520 real(kind=kind_phys) :: dxy(im),dxyp(im)
521 real(kind=kind_phys) :: ol4p(4),olp(im),od(im)
522 real(kind=kind_phys) :: taufb(im,km+1)
523
524 character(len=*), intent(out) :: errmsg
525 integer, intent(out) :: errflg
526
527 integer :: udtend, vdtend, Tdtend
528
529 ! Calculate inverse of gravitational acceleration
530 g_inv = 1./g
531
532 ! Initialize CCPP error handling variables
533 errmsg = ''
534 errflg = 0
535
536 ! Initialize local variables
537 var_temp2 = 0.
538 udtend = -1
539 vdtend = -1
540 tdtend = -1
541
542 if(ldiag3d) then
543 udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd)
544 vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd)
545 tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd)
546 endif
547
548
549 ! Initialize winds for sequential updating.
550 ! These are for optional sequential updating of the wind between the
551 ! LSGWD+BLOCKING and SSGWD+TOFD steps. They are placeholders
552 ! for the u1,v1 winds that are updated within the scheme if
553 ! ugwp_seq_update == T, otherwise, they retain the values
554 ! passed in to the subroutine.
555 uwnd1(:,:) = u1(:,:)
556 vwnd1(:,:) = v1(:,:)
557
558!--------------------------------------------------------------------
559! SCALE-ADPTIVE PARAMETER FROM GFS GWD SCHEME
560!--------------------------------------------------------------------
561! parameter (cdmb = 1.0) ! non-dim sub grid mtn drag Amp (*j*)
562! non-dim sub grid mtn drag Amp (*j*)
563! cdmb = 1.0/float(IMX/192)
564! cdmb = 192.0/float(IMX)
565 ! New cdmbgwd addition for GSL blocking drag
566 cdmb = 1.0
567 if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
568
570 kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2
571!
572! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
573!
574 if (imx > 0) then
575! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0)
576! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
577! cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
578! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192)
579! cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
580 cleff = 0.5e-5 / sqrt(float(imx)/192.0) ! this is inverse of CLEFF!
581! hmhj for ndsl
582! jw cleff = 0.1E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
583! cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
584! cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
585 endif
586 if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
587!--------------------------------------------------------------------
588! END SCALE-ADPTIVE PARAMETER SECTION
589!--------------------------------------------------------------------
590!
591!---- constants
592!
593 rcl = 1.
594 rcs = sqrt(rcl)
595 cs = 1. / sqrt(rcl)
596 csg = cs * g
597 lcap = km
598 lcapp1 = lcap + 1
599 fdir = mdir / (2.0*pi)
600 onebgrcs = 1._kind_phys/g*rcs
601
602 do i=1,im
603 if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3
604 xland(i)=1.0 !but land/water = (1/2) in this module
605 else
606 xland(i)=2.0
607 endif
608 enddo
609
610!--- calculate scale-aware tapering factors
611do i=1,im
612 if ( dx(i) .ge. dxmax_ms ) then
613 ls_taper(i) = 1.
614 else
615 if ( dx(i) .le. dxmin_ms) then
616 ls_taper(i) = 0.
617 else
618 ls_taper(i) = 0.5 * ( sin(pi*(dx(i)-0.5*(dxmax_ms+dxmin_ms))/ &
619 (dxmax_ms-dxmin_ms)) + 1. )
620 endif
621 endif
622enddo
623
624! Remove ss_tapering
625ss_taper(:) = 1.
626
627! SPP, if spp_gwd is 0, no perturbations are applied.
628if ( spp_gwd==1 ) then
629 do i = its,im
630 var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1)
631 varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1)
632 varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1)
633 enddo
634else
635 do i = its,im
636 var_stoch(i) = var(i)
637 varss_stoch(i) = varss(i)
638 varmax_fd_stoch(i) = varmax_fd
639 enddo
640endif
641
642!--- calculate length of grid for flow-blocking drag
643!
644do i=1,im
645 delx = dx(i)
646 dely = dx(i)
647 dxy4(i,1) = delx
648 dxy4(i,2) = dely
649 dxy4(i,3) = sqrt(delx*delx + dely*dely)
650 dxy4(i,4) = dxy4(i,3)
651 dxy4p(i,1) = dxy4(i,2)
652 dxy4p(i,2) = dxy4(i,1)
653 dxy4p(i,3) = dxy4(i,4)
654 dxy4p(i,4) = dxy4(i,3)
655enddo
656!
657!-----initialize arrays
658!
659 dtaux = 0.0
660 dtauy = 0.0
661 do i = its,im
662 klowtop(i) = 0
663 kbl(i) = 0
664 enddo
665!
666 do i = its,im
667 xn(i) = 0.0
668 yn(i) = 0.0
669 ubar(i) = 0.0
670 vbar(i) = 0.0
671 roll(i) = 0.0
672 taub(i) = 0.0
673 oa(i) = 0.0
674 ol(i) = 0.0
675 oass(i) = 0.0
676 olss(i) = 0.0
677 ulow(i) = 0.0
678 dtfac(i) = 1.0
679 rstoch(i) = 0.0
680 ldrag(i) = .false.
681 icrilv(i) = .false.
682 flag(i) = .true.
683 enddo
684
685 do k = kts,km
686 do i = its,im
687 usqj(i,k) = 0.0
688 bnv2(i,k) = 0.0
689 vtj(i,k) = 0.0
690 vtk(i,k) = 0.0
691 taup(i,k) = 0.0
692 taud_ms(i,k) = 0.0
693 taud_bl(i,k) = 0.0
694 dtaux2d(i,k) = 0.0
695 dtauy2d(i,k) = 0.0
696 enddo
697 enddo
698!
699 do i = its,im
700 xlinv(i) = 1.0/xl
701 enddo
702!
703! initialize array for flow-blocking drag
704!
705 taufb(1:im,1:km+1) = 0.0
706 komax(1:im) = 0
707!
708 rd_inv = 1./rd
709 do k = kts,km
710 do i = its,im
711 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
712 vtk(i,k) = vtj(i,k) / prslk(i,k)
713 ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3
714 enddo
715 enddo
716!
717! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
718!
719 !zq=0.
720 do k = kts,km
721 do i = its,im
722 !zq(i,k+1) = PHII(i,k+1)*g_inv
723 !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv
724 zl(i,k) = phil(i,k)*g_inv
725 enddo
726 enddo
727!
728! determine reference level: maximum of 2*var and pbl heights
729!
730 do i = its,im
731 zlowtop(i) = 2. * var_stoch(i)
732 enddo
733!
734 do i = its,im
735 kloop1(i) = .true.
736 enddo
737!
738 do k = kts+1,km
739 do i = its,im
740 if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
741 klowtop(i) = k+1
742 kloop1(i) = .false.
743 endif
744 enddo
745 enddo
746!
747 do i = its,im
748 kbl(i) = max(kpbl(i), klowtop(i))
749 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
750 enddo
751!
752! determine the level of maximum orographic height
753!
754 ! komax(:) = kbl(:)
755 komax(:) = klowtop(:) - 1 ! modification by NOAA/GSD March 2018
756!
757 do i = its,im
758 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
759 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
760 enddo
761!
762! compute low level averages within pbl
763!
764 do k = kts,kpblmax
765 do i = its,im
766 if (k.lt.kbl(i)) then
767 rcsks = rcs * del(i,k) * delks(i)
768 rdelks = del(i,k) * delks(i)
769 ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
770 vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
771 roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
772 endif
773 enddo
774 enddo
775!
776! figure out low-level horizontal wind direction
777!
778! nwd 1 2 3 4 5 6 7 8
779! wd w s sw nw e n ne se
780!
781 do i = its,im
782 wdir = atan2(ubar(i),vbar(i)) + pi
783 idir = mod(nint(fdir*wdir),mdir) + 1
784 nwd = nwdir(idir)
785 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
786 ol(i) = ol4(i,mod(nwd-1,4)+1)
787 ! Repeat for small-scale gwd
788 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
789 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
790
791!
792!----- compute orographic width along (ol) and perpendicular (olp)
793!----- the direction of wind
794!
795 ol4p(1) = ol4(i,2)
796 ol4p(2) = ol4(i,1)
797 ol4p(3) = ol4(i,4)
798 ol4p(4) = ol4(i,3)
799 olp(i) = ol4p(mod(nwd-1,4)+1)
800!
801!----- compute orographic direction (horizontal orographic aspect ratio)
802!
803 od(i) = olp(i)/max(ol(i),olmin)
804 od(i) = min(od(i),odmax)
805 od(i) = max(od(i),odmin)
806!
807!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
808!
809 dxy(i) = dxy4(i,mod(nwd-1,4)+1)
810 dxyp(i) = dxy4p(i,mod(nwd-1,4)+1)
811 enddo
812!
813! END INITIALIZATION; BEGIN GWD CALCULATIONS:
814!
815
816IF ( (do_gsl_drag_ls_bl).and. &
817 (gwd_opt_ms.or.gwd_opt_bl) ) then
818
819 do i=its,im
820
821 rdxzb(i) = 0.0
822
823 if ( ls_taper(i).GT.1.e-02 ) then
824
825!
826!--- saving richardson number in usqj for migwdi
827!
828 do k = kts,km-1
829 ti = 2.0 / (t1(i,k)+t1(i,k+1))
830 rdz = 1./(zl(i,k+1) - zl(i,k))
831 tem1 = u1(i,k) - u1(i,k+1)
832 tem2 = v1(i,k) - v1(i,k+1)
833 dw2 = rcl*(tem1*tem1 + tem2*tem2)
834 shr2 = max(dw2,dw2min) * rdz * rdz
835 bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
836 usqj(i,k) = max(bvf2/shr2,rimin)
837 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
838 bnv2(i,k) = max( bnv2(i,k), bnv2min )
839 enddo
840!
841!----compute the "low level" or 1/3 wind magnitude (m/s)
842!
843 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
844 rulow(i) = 1./ulow(i)
845!
846 do k = kts,km-1
847 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
848 + (v1(i,k)+v1(i,k+1)) * vbar(i))
849 velco(i,k) = velco(i,k) * rulow(i)
850 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
851 velco(i,k) = veleps
852 endif
853 enddo
854!
855! no drag when critical level in the base layer
856!
857 ldrag(i) = velco(i,1).le.0.
858!
859! no drag when velco.lt.0
860!
861 do k = kpblmin,kpblmax
862 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
863 enddo
864!
865! no drag when bnv2.lt.0
866!
867 do k = kts,kpblmax
868 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
869 enddo
870!
871!-----the low level weighted average ri is stored in usqj(1,1; im)
872!-----the low level weighted average n**2 is stored in bnv2(1,1; im)
873!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
874!---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
875!
876 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
877 bnv2(i,1) = wtkbj * bnv2(i,1)
878 usqj(i,1) = wtkbj * usqj(i,1)
879!
880 do k = kpblmin,kpblmax
881 if (k .lt. kbl(i)) then
882 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
883 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
884 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
885 endif
886 enddo
887!
888 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
889 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
890 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
891! Check if mesoscale gravity waves will propagate vertically or be evanescent
892! and not impart a drag force -- consider the maximum sub-grid horizontal
893! topographic wavelength to be one-half the horizontal grid spacing -- calculate
894! ksmax accordingly
895 ksmax = 4.0*pi/dx(i) ! based on wavelength = 0.5*dx(i)
896 if ( bnv2(i,1).gt.0.0 ) then
897 ldrag(i) = ldrag(i) .or. sqrt(bnv2(i,1))*rulow(i).lt.ksmax
898 endif
899!
900! set all ri low level values to the low level value
901!
902 do k = kpblmin,kpblmax
903 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
904 enddo
905!
906 if (.not.ldrag(i)) then
907 bnv(i) = sqrt( bnv2(i,1) )
908 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
909 fr(i) = min(fr(i),frmax)
910 xn(i) = ubar(i) * rulow(i)
911 yn(i) = vbar(i) * rulow(i)
912 endif
913!
914! compute the base level stress and store it in taub
915! calculate enhancement factor, number of mountains & aspect
916! ratio const. use simplified relationship between standard
917! deviation & critical hgt
918
919 if (.not. ldrag(i)) then
920 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
921 efact = min( max(efact,efmin), efmax )
922!!!!!!! cleff (effective grid length) is highly tunable parameter
923!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
924!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.)
925!WRF cleff = 3. * max(dx(i),cleff)
926 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
927!WRF xlinv(i) = coefm(i) / cleff
928 xlinv(i) = coefm(i) * cleff
929 tem = fr(i) * fr(i) * oc1(i)
930 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
931 if ( gwd_opt_ms ) then
932 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
933 * ulow(i) * gfobnv * efact
934 else ! We've gotten what we need for the blocking scheme
935 taub(i) = 0.0
936 end if
937 else
938 taub(i) = 0.0
939 xn(i) = 0.0
940 yn(i) = 0.0
941 endif
942
943 endif ! (ls_taper(i).GT.1.E-02)
944
945 enddo ! do i=its,im
946
947ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ms.or.gwd_opt_bl)
948
949
950
951
952!=======================================================
953! Mesoscale GWD + blocking
954!=======================================================
955IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ms) ) THEN
956
957 do i=its,im
958
959 if ( ls_taper(i).GT.1.e-02 ) then
960
961!
962! now compute vertical structure of the stress.
963 do k = kts,kpblmax
964 if (k .le. kbl(i)) taup(i,k) = taub(i)
965 enddo
966!
967 do k = kpblmin, km-1 ! vertical level k loop!
968 kp1 = k + 1
969!
970! unstablelayer if ri < ric
971! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
972! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
973!
974 if (k .ge. kbl(i)) then
975 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
976 .or. (velco(i,k) .le. 0.0)
977 brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
978 brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
979 endif
980!
981 if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
982 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
983 temv = 1.0 / velco(i,k)
984 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* &
985 velco(i,k)*0.5
986 hd = sqrt(taup(i,k) / tem1)
987 fro = brvf(i) * hd * temv
988!
989! rim is the minimum-richardson number by shutts (1985)
990 tem2 = sqrt(usqj(i,k))
991 tem = 1. + tem2 * fro
992 rim = usqj(i,k) * (1.-fro) / (tem * tem)
993!
994! check stability to employ the 'saturation hypothesis'
995! of lindzen (1981) except at tropospheric downstream regions
996!
997 if (rim .le. ric) then ! saturation hypothesis!
998 if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then
999 temc = 2.0 + 1.0 / tem2
1000 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
1001 taup(i,kp1) = tem1 * hd * hd
1002 endif
1003 else ! no wavebreaking!
1004 taup(i,kp1) = taup(i,k)
1005 endif
1006 endif
1007 endif
1008 enddo
1009!
1010 if(lcap.lt.km) then
1011 do klcap = lcapp1,km
1012 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
1013 enddo
1014 endif
1015
1016 endif ! if ( ls_taper(i).GT.1.E-02 )
1017
1018 enddo ! do i=its,im
1019
1020ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ms)
1021!===============================================================
1022!COMPUTE BLOCKING COMPONENT
1023!===============================================================
1024IF ( do_gsl_drag_ls_bl .and. gwd_opt_bl ) THEN
1025
1026 do i=its,im
1027
1028 if ( ls_taper(i).GT.1.e-02 ) then
1029
1030 if (.not.ldrag(i)) then
1031!
1032!------- determine the height of flow-blocking layer
1033!
1034 kblk = 0
1035 pe = 0.0
1036 do k = km, kpblmin, -1
1037 if(kblk.eq.0 .and. k.le.komax(i)) then
1038 pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))* &
1039 del(i,k)/g/ro(i,k)
1040 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
1041!
1042!---------- apply flow-blocking drag when pe >= ke
1043!
1044 if(pe.ge.ke) then
1045 kblk = k
1046 kblk = min(kblk,kbl(i))
1047 zblk = zl(i,kblk)-zl(i,kts)
1048 rdxzb(i) = real(k,kind=kind_phys)
1049 endif
1050 endif
1051 enddo
1052 if(kblk.ne.0) then
1053!
1054!--------- compute flow-blocking stress
1055!
1056 cd = max(2.0-1.0/od(i),0.0)
1057 ! New cdmbgwd addition for GSL blocking drag
1058 taufb(i,kts) = cdmb * 0.5 * roll(i) * coefm(i) / &
1059 max(dxmax_ms,dxy(i))**2 * cd * dxyp(i) * &
1060 olp(i) * zblk * ulow(i)**2
1061 tautem = taufb(i,kts)/float(kblk-kts)
1062 do k = kts+1, kblk
1063 taufb(i,k) = taufb(i,k-1) - tautem
1064 enddo
1065!
1066!----------sum orographic GW stress and flow-blocking stress
1067!
1068 ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now
1069 endif
1070
1071 endif ! if (.not.ldrag(i))
1072
1073 endif ! if ( ls_taper(i).GT.1.E-02 )
1074
1075 enddo ! do i=its,im
1076
1077ENDIF ! IF ( do_gsl_drag_ls_bl .and. gwd_opt_bl )
1078!===========================================================
1079IF ( (do_gsl_drag_ls_bl) .and. &
1080 (gwd_opt_ms .OR. gwd_opt_bl) ) THEN
1081
1082 do i=its,im
1083
1084 if ( ls_taper(i) .GT. 1.e-02 ) then
1085
1086!
1087! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
1088!
1089! First, set taup (momentum flux) at model top equal to that of the layer
1090! interface just below the top, i.e., taup(km)
1091! The idea is to allow the momentum flux to go out the 'top'. This
1092! ensures there is no GWD force at the top layer.
1093!
1094 taup(i,km+1) = taup(i,km)
1095 do k = kts,km
1096 taud_ms(i,k) = (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
1097 taud_bl(i,k) = (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
1098 enddo
1099!
1100!
1101! if the gravity wave drag + blocking would force a critical line
1102! in the layers below pressure-based 'sigma' level = sgmalolev during the next deltim
1103! timestep, then only apply drag until that critical line is reached, i.e.,
1104! reduce drag to limit resulting wind components to zero
1105! Note: 'sigma' = prsi(k)/prsi(k=1), where prsi(k=1) is the surface pressure
1106!
1107 do k = kts,kpblmax-1
1108 if (prsi(i,k).ge.sgmalolev*prsi(i,1)) then
1109 if ((taud_ms(i,k)+taud_bl(i,k)).ne.0.) &
1110 dtfac(i) = min(dtfac(i),abs(velco(i,k) &
1111 /(deltim*rcs*(taud_ms(i,k)+taud_bl(i,k)))))
1112 else
1113 exit
1114 endif
1115 enddo
1116!
1117 do k = kts,km
1118
1119 ! Check if well into mesosphere -- if so, perform similar reduction of
1120 ! velocity tendency due to mesoscale GWD to prevent sudden reversal of
1121 ! wind direction (similar to above)
1122 dtfac_meso = 1.0
1123 if (prsl(i,k).le.plolevmeso) then
1124 if (taud_ms(i,k).ne.0.) &
1125 dtfac_meso = min(dtfac_meso,facmeso*abs(velco(i,k) &
1126 /(deltim*rcs*taud_ms(i,k))))
1127 end if
1128
1129 taud_ms(i,k) = taud_ms(i,k)*dtfac(i)*dtfac_meso* &
1130 ls_taper(i) *(1.-rstoch(i))
1131 taud_bl(i,k) = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i))
1132
1133 dtaux = taud_ms(i,k) * xn(i)
1134 dtauy = taud_ms(i,k) * yn(i)
1135 dtauxb = taud_bl(i,k) * xn(i)
1136 dtauyb = taud_bl(i,k) * yn(i)
1137
1138 !add blocking and mesoscale contributions to tendencies
1139 tmp1 = dtaux + dtauxb
1140 tmp2 = dtauy + dtauyb
1141 dudt(i,k) = tmp1 + dudt(i,k)
1142 dvdt(i,k) = tmp2 + dvdt(i,k)
1143
1144 ! Update winds if sequential updating is selected
1145 ! and SSGWD and TOFD will be calculated
1146 ! Note: uwnd1 and vwnd1 replace u1 and u2,respectively,
1147 ! for the SSGWD and TOFD calculations
1148 if ( ugwp_seq_update .and. (do_gsl_drag_ss.or.do_gsl_drag_tofd) ) then
1149 uwnd1(i,k) = uwnd1(i,k) + tmp1*deltim
1150 vwnd1(i,k) = vwnd1(i,k) + tmp2*deltim
1151 endif
1152
1153 if ( gsd_diss_ht_opt ) then
1154 ! Calculate dissipation heating
1155 ! Initial kinetic energy (at t0-dt)
1156 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
1157 ! Kinetic energy after wave-breaking/flow-blocking
1158 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
1159 (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
1160 ! Modify theta tendency
1161 dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim
1162 if ( tdtend>0 ) then
1163 dtend(i,k,tdtend) = dtend(i,k,tdtend) + max((eng0-eng1),0.0)/cp
1164 endif
1165 endif
1166
1167 dusfc(i) = dusfc(i) - onebgrcs * ( taud_ms(i,k)*xn(i)*del(i,k) + &
1168 taud_bl(i,k)*xn(i)*del(i,k) )
1169 dvsfc(i) = dvsfc(i) - onebgrcs * ( taud_ms(i,k)*yn(i)*del(i,k) + &
1170 taud_bl(i,k)*yn(i)*del(i,k) )
1171 if(udtend>0) then
1172 dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ms(i,k) * &
1173 xn(i) + taud_bl(i,k) * xn(i)) * deltim
1174 endif
1175 if(vdtend>0) then
1176 dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ms(i,k) * &
1177 yn(i) + taud_bl(i,k) * yn(i)) * deltim
1178 endif
1179
1180 enddo
1181
1182 if ( ldiag_ugwp ) then
1183 do k = kts,km
1184 dtaux2d_ms(i,k) = taud_ms(i,k) * xn(i)
1185 dtauy2d_ms(i,k) = taud_ms(i,k) * yn(i)
1186 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
1187 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
1188 dusfc_ms(i) = dusfc_ms(i) - onebgrcs * dtaux2d_ms(i,k) * del(i,k)
1189 dvsfc_ms(i) = dvsfc_ms(i) - onebgrcs * dtauy2d_ms(i,k) * del(i,k)
1190 dusfc_bl(i) = dusfc_bl(i) - onebgrcs * dtaux2d_bl(i,k) * del(i,k)
1191 dvsfc_bl(i) = dvsfc_bl(i) - onebgrcs * dtauy2d_bl(i,k) * del(i,k)
1192 enddo
1193 endif
1194
1195 endif ! if ( ls_taper(i) .GT. 1.E-02 )
1196
1197 enddo ! do i=its,im
1198
1199ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ms .OR. gwd_opt_bl)
1200
1201
1202!====================================================================
1203! Calculate small-scale gravity wave drag for stable boundary layer
1204!====================================================================
1205 xnbv=0.
1206 tauwavex0=0.
1207 tauwavey0=0.
1208 density=1.2
1209 utendwave=0.
1210 vtendwave=0.
1211!
1212IF ( do_gsl_drag_ss ) THEN
1213
1214 do i=its,im
1215
1216 if ( ss_taper(i).GT.1.e-02 ) then
1217 !
1218 ! calculating potential temperature
1219 !
1220 do k = kts,km
1221 thx(i,k) = t1(i,k)/prslk(i,k)
1222 enddo
1223 !
1224 do k = kts,km
1225 tvcon = (1.+fv*q1(i,k))
1226 thvx(i,k) = thx(i,k)*tvcon
1227 enddo
1228
1229 hpbl2 = hpbl(i)+10.
1230 kpbl2 = kpbl(i)
1231 !kvar = MIN(kpbl, k-level of var)
1232 kvar = 1
1233 do k=kts+1,max(kpbl(i),kts+1)
1234! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then
1235 IF (zl(i,k)>300.) then
1236 kpbl2 = k
1237 IF (k == kpbl(i)) then
1238 hpbl2 = hpbl(i)+10.
1239 ELSE
1240 hpbl2 = zl(i,k)+10.
1241 ENDIF
1242 exit
1243 ENDIF
1244 enddo
1245 if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
1246 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
1247 ! Modify xlinv to represent wave number of "typical" small-scale topography
1248! cleff_ss = 3. * max(dx(i),cleff_ss)
1249! cleff_ss = 10. * max(dxmax_ss,cleff_ss)
1250! cleff_ss = 0.1 * 12000.
1251 xlinv(i) = 0.001*pi ! 2km horizontal wavelength
1252 !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts)))
1253 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
1254 !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i))
1255 xnbv=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
1256!
1257 ! check for possibility of vertical wave propagation
1258 ! (avoids division by zero if uwnd1(i,kpbl2).eq.0.)
1259 if (uwnd1(i,kpbl2).eq.0.) then
1260 prop_test = .true.
1261 elseif (abs(xnbv/uwnd1(i,kpbl2)).gt.xlinv(i)) then
1262 prop_test = .true.
1263 else
1264 prop_test = .false.
1265 endif
1266 if (prop_test) then
1267 ! Remove limit on varss_stoch
1268 var_temp = varss_stoch(i)
1269 ! Note: This is a semi-implicit treatment of the time differencing
1270 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
1271 tauwavex0=var_temp2*uwnd1(i,kvar)/(1.+var_temp2*deltim)
1272 tauwavex0=tauwavex0*ss_taper(i)
1273 else
1274 tauwavex0=0.
1275 endif
1276!
1277 ! check for possibility of vertical wave propagation
1278 ! (avoids division by zero if vwnd1(i,kpbl2).eq.0.)
1279 if (vwnd1(i,kpbl2).eq.0.) then
1280 prop_test = .true.
1281 elseif (abs(xnbv/vwnd1(i,kpbl2)).gt.xlinv(i)) then
1282 prop_test = .true.
1283 else
1284 prop_test = .false.
1285 endif
1286 if (prop_test) then
1287 ! Remove limit on varss_stoch
1288 var_temp = varss_stoch(i)
1289 ! Note: This is a semi-implicit treatment of the time differencing
1290 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
1291 tauwavey0=var_temp2*vwnd1(i,kvar)/(1.+var_temp2*deltim)
1292 tauwavey0=tauwavey0*ss_taper(i)
1293 else
1294 tauwavey0=0.
1295 endif
1296
1297 do k=kts,kpbl(i) !MIN(kpbl2+1,km-1)
1298!original
1299 !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
1300 !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
1301!new
1302 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
1303 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
1304!mod-to be used in HRRRv3/RAPv4
1305 !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2
1306 !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2
1307 enddo
1308 endif
1309 endif
1310
1311 do k = kts,km
1312 dudt(i,k) = dudt(i,k) + utendwave(i,k)
1313 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
1314 dusfc(i) = dusfc(i) - onebgrcs * utendwave(i,k) * del(i,k)
1315 dvsfc(i) = dvsfc(i) - onebgrcs * vtendwave(i,k) * del(i,k)
1316 enddo
1317 if(udtend>0) then
1318 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim
1319 endif
1320 if(vdtend>0) then
1321 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim
1322 endif
1323 if ( ldiag_ugwp ) then
1324 do k = kts,km
1325 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
1326 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
1327 dtaux2d_ss(i,k) = utendwave(i,k)
1328 dtauy2d_ss(i,k) = vtendwave(i,k)
1329 enddo
1330 endif
1331
1332 endif ! if (ss_taper(i).GT.1.E-02)
1333
1334 enddo ! i=its,im
1335
1336ENDIF ! if (do_gsl_drag_ss)
1337
1338
1339!===================================================================
1340! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16):
1341!===================================================================
1342IF ( do_gsl_drag_tofd ) THEN
1343
1344 do i=its,im
1345
1346 if ( ss_taper(i).GT.1.e-02 ) then
1347
1348 utendform=0.
1349 vtendform=0.
1350
1351 IF ((xland(i)-1.5) .le. 0.) then
1352 !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
1353 var_temp = min(varss_stoch(i),varmax_fd_stoch(i)) + &
1354 max(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i)))
1355 a1=0.00026615161*var_temp**2
1356! a1=0.00026615161*MIN(varss(i),varmax)**2
1357! a1=0.00026615161*(0.5*varss(i))**2
1358 ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
1359 a2=a1*0.005363
1360 ! Beljaars H_efold
1361 h_efold = 1500.
1362 DO k=kts,km
1363 wsp=sqrt(uwnd1(i,k)**2 + vwnd1(i,k)**2)
1364 ! Note: In Beljaars et al. (2004):
1365 ! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
1366 ! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys
1367 var_temp = alpha_fd*coeff_fd*exp(-(zl(i,k)/h_efold)**1.5)*a2* &
1368 zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
1369 ! Note: This is a semi-implicit treatment of the time differencing
1370 ! per Beljaars et al. (2004, QJRMS)
1371 utendform(i,k) = - var_temp*wsp*uwnd1(i,k)/(1. + var_temp*deltim*wsp)
1372 vtendform(i,k) = - var_temp*wsp*vwnd1(i,k)/(1. + var_temp*deltim*wsp)
1373 !IF(zl(i,k) > 4000.) exit
1374 ENDDO
1375 ENDIF
1376
1377 do k = kts,km
1378 dudt(i,k) = dudt(i,k) + utendform(i,k)
1379 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
1380 dusfc(i) = dusfc(i) - onebgrcs * utendform(i,k) * del(i,k)
1381 dvsfc(i) = dvsfc(i) - onebgrcs * vtendform(i,k) * del(i,k)
1382 enddo
1383 if(udtend>0) then
1384 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim
1385 endif
1386 if(vdtend>0) then
1387 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim
1388 endif
1389 if ( ldiag_ugwp ) then
1390 do k = kts,km
1391 dtaux2d_fd(i,k) = utendform(i,k)
1392 dtauy2d_fd(i,k) = vtendform(i,k)
1393 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
1394 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
1395 enddo
1396 endif
1397
1398 endif ! if (ss_taper(i).GT.1.E-02)
1399
1400 enddo ! i=its,im
1401
1402ENDIF ! if (do_gsl_drag_tofd)
1403
1404
1405
1406if ( ldiag_ugwp ) then
1407 ! Finalize dusfc and dvsfc diagnostics for gsl small-scale drag components
1408 dusfc_ss(:) = -onebgrcs * dusfc_ss(:)
1409 dvsfc_ss(:) = -onebgrcs * dvsfc_ss(:)
1410 dusfc_fd(:) = -onebgrcs * dusfc_fd(:)
1411 dvsfc_fd(:) = -onebgrcs * dvsfc_fd(:)
1412endif
1413!
1414 return
1415 end subroutine drag_suite_run
1416!-------------------------------------------------------------------
1417
1418 subroutine drag_suite_psl( &
1419 & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, &
1420 & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, &
1421 & var,oc1,oa4,ol4, &
1422 & varss,oc1ss,oa4ss,ol4ss, &
1423 & THETA,SIGMA,GAMMA,ELVMAX, &
1424 & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, &
1425 & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, &
1426 & dusfc,dvsfc, &
1427 & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
1428 & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
1429 & slmsk,br1,hpbl,vtype, &
1430 & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
1431 & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
1432 & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
1433 & psl_gwd_dx_factor, &
1434 & dtend, dtidx, index_of_process_orographic_gwd, &
1435 & index_of_temperature, index_of_x_wind, &
1436 & index_of_y_wind, ldiag3d, ldiag_ugwp, ugwp_seq_update, &
1437 & spp_wts_gwd, spp_gwd, errmsg, errflg)
1438
1439! ********************************************************************
1440! -----> I M P L E M E N T A T I O N V E R S I O N <----------
1441!
1442! ----- This code -----
1443!begin WRF code
1444
1445! this code handles the time tendencies of u v due to the effect of mountain
1446! induced gravity wave drag from sub-grid scale orography. this routine
1447! not only treats the traditional upper-level wave breaking due to mountain
1448! variance (alpert 1988), but also the enhanced lower-tropospheric wave
1449! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
1450! thus, in addition to the terrain height data in a model grid box,
1451! additional 10-2d topographic statistics files are needed, including
1452! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
1453! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
1454! hong (1999). the current scheme was implmented as in hong et al.(2008)
1455!
1456! Originally coded by song-you hong and young-joon kim and implemented by song-you hong
1457!
1458! program history log:
1459! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle
1460! with blocked height by dividing streamline theory
1461! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale
1462! orographic grabity wave drag:
1463! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the
1464! topographic form drag of Beljaars et al. (2004, QJRMS)
1465! Activation of each component is done by specifying the integer-parameters
1466! (defined below) to 0: inactive or 1: active
1467! gwd_opt_ls = 0 or 1: large-scale
1468! gwd_opt_bl = 0 or 1: blocking drag
1469! gwd_opt_ss = 0 or 1: small-scale gravity wave drag
1470! gwd_opt_fd = 0 or 1: topographic form drag
1471! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating
1472! gsd_diss_ht_opt = 0: dissipation heating off
1473! gsd_diss_ht_opt = 1: dissipation heating on
1474! 2020-08-25 Michael Toy changed logic control for drag component selection
1475! for CCPP.
1476! Namelist options:
1477! do_gsl_drag_ls_bl - logical flag for large-scale GWD + blocking
1478! do_gsl_drag_ss - logical flag for small-scale GWD
1479! do_gsl_drag_tofd - logical flag for turbulent form drag
1480! Compile-time options (same as before):
1481! gwd_opt_ls = 0 or 1: large-scale GWD
1482! gwd_opt_bl = 0 or 1: blocking drag
1483!
1484! References:
1485! Choi and Hong (2015) J. Geophys. Res.
1486! Hong et al. (2008), wea. and forecasting
1487! Kim and Doyle (2005), Q. J. R. Meteor. Soc.
1488! Kim and Arakawa (1995), j. atmos. sci.
1489! Alpert et al. (1988), NWP conference.
1490! Hong (1999), NCEP office note 424.
1491! Steeneveld et al (2008), JAMC
1492! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc.
1493! Beljaars et al. (2004), Q. J. R. Meteor. Soc.
1494!
1495! notice : comparible or lower resolution orography files than model resolution
1496! are desirable in preprocess (wps) to prevent weakening of the drag
1497!-------------------------------------------------------------------------------
1498!
1499! input
1500! dudt (im,km) non-lin tendency for u wind component
1501! dvdt (im,km) non-lin tendency for v wind component
1502! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt
1503! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt
1504! t1(im,km) temperature deg k at t0-dt
1505! q1(im,km) specific humidity at t0-dt
1506! deltim time step secs
1507! del(km) positive increment of pressure across layer (pa)
1508! KPBL(IM) is the index of the top layer of the PBL
1509! ipr & lprnt for diagnostics
1510!
1511! output
1512! dudt, dvdt wind tendency due to gwdo
1513! dTdt
1514!
1515!-------------------------------------------------------------------------------
1516
1517!end wrf code
1518!----------------------------------------------------------------------C
1519! USE
1520! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES)
1521!
1522! PURPOSE
1523! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
1524! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V
1525! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
1526! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
1527! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
1528! CRITICAL LEVELS
1529!
1530!
1531! ********************************************************************
1532 USE machine , ONLY : kind_phys
1533 implicit none
1534
1535 ! Interface variables
1536 integer, intent(in) :: im, km, imx, kdt, ipr, me, master
1537 integer, intent(in) :: gwd_opt
1538 logical, intent(in) :: lprnt
1539 integer, intent(in) :: KPBL(:)
1540 real(kind=kind_phys), intent(in) :: deltim, g, cp, rd, rv, &
1541 & cdmbgwd(:), alpha_fd
1542 real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
1543 logical, intent(in) :: ldiag3d
1544 integer, intent(in) :: dtidx(:,:), index_of_temperature, &
1545 & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind
1546
1547 integer :: kpblmax
1548 integer, parameter :: ims=1, kms=1, its=1, kts=1
1549 real(kind=kind_phys), intent(in) :: fv, pi
1550 real(kind=kind_phys) :: rcl, cdmb
1551 real(kind=kind_phys) :: g_inv, g_cp, rd_inv
1552
1553 real(kind=kind_phys), intent(inout) :: &
1554 & dudt(:,:),dvdt(:,:), &
1555 & dtdt(:,:)
1556 real(kind=kind_phys), intent(out) :: rdxzb(:)
1557 real(kind=kind_phys), intent(in) :: &
1558 & u1(:,:),v1(:,:), &
1559 & t1(:,:),q1(:,:), &
1560 & phii(:,:),prsl(:,:), &
1561 & prslk(:,:),phil(:,:)
1562 real(kind=kind_phys), intent(in) :: prsi(:,:), &
1563 & del(:,:)
1564 real(kind=kind_phys), intent(in) :: var(:),oc1(:), &
1565 & oa4(:,:),ol4(:,:), &
1566 & dx(:)
1567 real(kind=kind_phys), intent(in), optional :: varss(:),oc1ss(:), &
1568 & oa4ss(:,:),ol4ss(:,:)
1569 real(kind=kind_phys), intent(in) :: theta(:),sigma(:), &
1570 & gamma(:),elvmax(:)
1571
1572! added for small-scale orographic wave drag
1573 real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx
1574 integer, intent(in) :: vtype(:)
1575 real(kind=kind_phys), intent(in) :: br1(:), &
1576 & hpbl(:), &
1577 & slmsk(:)
1578 real(kind=kind_phys), dimension(im) :: govrth,xland
1579 !real(kind=kind_phys), dimension(im,km) :: dz2
1580 real(kind=kind_phys) :: tauwavex0,tauwavey0, &
1581 & xnbv,density,tvcon,hpbl2
1582 integer :: kpbl2,kvar
1583 !real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g
1584 real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g
1585
1586!SPP
1587 real(kind=kind_phys), dimension(im) :: var_stoch, varss_stoch, &
1588 varmax_ss_stoch, varmax_fd_stoch
1589 real(kind=kind_phys), intent(in), optional :: spp_wts_gwd(:,:)
1590 integer, intent(in) :: spp_gwd
1591
1592 real(kind=kind_phys), dimension(im) :: rstoch
1593
1594!Output:
1595 real(kind=kind_phys), intent(out) :: &
1596 & dusfc(:), dvsfc(:)
1597!Output (optional):
1598 real(kind=kind_phys), intent(out), optional :: &
1599 & dusfc_ls(:),dvsfc_ls(:), &
1600 & dusfc_bl(:),dvsfc_bl(:), &
1601 & dusfc_ss(:),dvsfc_ss(:), &
1602 & dusfc_fd(:),dvsfc_fd(:)
1603 real(kind=kind_phys), intent(out), optional :: &
1604 & dtaux2d_ls(:,:),dtauy2d_ls(:,:), &
1605 & dtaux2d_bl(:,:),dtauy2d_bl(:,:), &
1606 & dtaux2d_ss(:,:),dtauy2d_ss(:,:), &
1607 & dtaux2d_fd(:,:),dtauy2d_fd(:,:)
1608
1609!Misc arrays
1610 real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d
1611
1612!-------------------------------------------------------------------------
1613! Flags to regulate the activation of specific components of drag suite:
1614! Each component is tapered off automatically as a function of dx, so best to
1615! keep them activated (.true.).
1616 logical, intent(in) :: &
1617 do_gsl_drag_ls_bl, & ! large-scale gravity wave drag and blocking
1618 do_gsl_drag_ss, & ! small-scale gravity wave drag (Steeneveld et al. 2008)
1619 do_gsl_drag_tofd ! form drag (Beljaars et al. 2004, QJRMS)
1620! Flag for diagnostic outputs
1621 logical, intent(in) :: ldiag_ugwp
1622
1623! Flag for sequential update of u and v between
1624! LSGWD + BLOCKING and SSGWD + TOFD calculations
1625 logical, intent(in) :: ugwp_seq_update
1626!
1627! Additional flags
1628 integer, parameter :: &
1629 gwd_opt_ls = 1, & ! large-scale gravity wave drag
1630 gwd_opt_bl = 1, & ! blocking drag
1631 gsd_diss_ht_opt = 0
1632
1633! Parameters for bounding the scale-adaptive variability:
1634! Small-scale GWD + turbulent form drag
1635 real(kind=kind_phys), parameter :: dxmin_ss = 1000., &
1636 & dxmax_ss = 12000. ! min,max range of tapering (m)
1637! Large-scale GWD + blocking
1638 real(kind=kind_phys), parameter :: dxmin_ls = 3000., &
1639 & dxmax_ls = 13000. ! min,max range of tapering (m)
1640 real(kind=kind_phys), dimension(im) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-)
1641!
1642! Variables for limiting topographic standard deviation (var)
1643 real(kind=kind_phys), parameter :: varmax_ss = 50., &
1644 varmax_fd = 150., &
1645 beta_ss = 0.1, &
1646 beta_fd = 0.2
1647 real(kind=kind_phys) :: var_temp, var_temp2
1648
1649! added Beljaars orographic form drag
1650 real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
1651 real(kind=kind_phys) :: a1,a2,wsp
1652 real(kind=kind_phys) :: h_efold
1653 real(kind=kind_phys), parameter :: coeff_fd = 6.325e-3
1654
1655! multification factor of standard deviation : ! larger drag with larger value
1656!!! real(kind=kind_phys), parameter :: psl_gwd_dx_factor = 6.0
1657 real(kind=kind_phys), intent(in) :: psl_gwd_dx_factor
1658
1659! critical richardson number for wave breaking : ! larger drag with larger value
1660 real(kind=kind_phys), parameter :: ric = 0.25
1661 real(kind=kind_phys), parameter :: dw2min = 1.
1662 real(kind=kind_phys), parameter :: rimin = -100.
1663 real(kind=kind_phys), parameter :: bnv2min = 1.0e-5
1664 real(kind=kind_phys), parameter :: efmin = 0.0
1665 real(kind=kind_phys), parameter :: efmax = 10.0
1666 real(kind=kind_phys), parameter :: xl = 4.0e4
1667 real(kind=kind_phys), parameter :: critac = 1.0e-5
1668 real(kind=kind_phys), parameter :: gmax = 1.
1669 real(kind=kind_phys), parameter :: veleps = 1.0
1670 real(kind=kind_phys), parameter :: factop = 0.5
1671 real(kind=kind_phys), parameter :: frc = 1.0
1672 real(kind=kind_phys), parameter :: ce = 0.8
1673 real(kind=kind_phys), parameter :: cg = 1.0
1674! real(kind=kind_phys), parameter :: var_min = 100.0
1675 real(kind=kind_phys), parameter :: var_min = 10.0
1676 real(kind=kind_phys), parameter :: hmt_min = 50.
1677 real(kind=kind_phys), parameter :: oc_min = 1.0
1678 real(kind=kind_phys), parameter :: oc_max = 10.0
1679! 7.5 mb -- 33 km ... 0.01 kgm-3 reduce gwd drag above cutoff level
1680 real(kind=kind_phys), parameter :: pcutoff = 7.5e2
1681! 0.76 mb -- 50 km ...0.001 kgm-3 --- 0.1 mb 65 km 0.0001 kgm-3
1682 real(kind=kind_phys), parameter :: pcutoff_den = 0.01 !
1683
1684 integer,parameter :: kpblmin = 2
1685
1686!
1687! local variables
1688!
1689 integer :: i,j,k,lcap,lcapp1,nwd,idir, &
1690 klcap,kp1
1691!
1692 real(kind=kind_phys) :: rcs,csg,fdir,cs, &
1693 rcsks,wdir,ti,rdz,tem2,dw2,shr2, &
1694 bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, &
1695 rim,temc,tem1,efact,temv,dtaux,dtauy, &
1696 dtauxb,dtauyb,eng0,eng1
1697 real(kind=kind_phys) :: denfac
1698!
1699 logical :: ldrag(im),icrilv(im), &
1700 flag(im)
1701!
1702 real(kind=kind_phys) :: invgrcs
1703!
1704 real(kind=kind_phys) :: taub(im),taup(im,km+1), &
1705 xn(im),yn(im), &
1706 ubar(im),vbar(im), &
1707 fr(im),ulow(im), &
1708 rulow(im),bnv(im), &
1709 oa(im),ol(im),oc(im), &
1710 oass(im),olss(im), &
1711 roll(im),dtfac(im,km), &
1712 brvf(im),xlinv(im), &
1713 delks(im),delks1(im), &
1714 bnv2(im,km),usqj(im,km), &
1715 taud_ls(im,km),taud_bl(im,km), &
1716 ro(im,km), &
1717 vtk(im,km),vtj(im,km), &
1718 zlowtop(im),velco(im,km-1), &
1719 coefm(im),coefm_ss(im)
1720 real(kind=kind_phys) :: cleff(im),cleff_ss(im)
1721!
1722 integer :: kbl(im),klowtop(im)
1723 integer,parameter :: mdir=8
1724 !integer :: nwdir(mdir)
1725 !data nwdir/6,7,5,8,2,3,1,4/
1726 integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/)
1727!
1728! variables for flow-blocking drag
1729!
1730 real(kind=kind_phys),parameter :: frmax = 10.
1731 real(kind=kind_phys),parameter :: olmin = 1.e-5
1732 real(kind=kind_phys),parameter :: odmin = 0.1
1733 real(kind=kind_phys),parameter :: odmax = 10.
1734 real(kind=kind_phys),parameter :: cdmin = 0.0
1735 integer :: komax(im),kbmax(im),kblk(im)
1736 real(kind=kind_phys) :: hmax(im)
1737 real(kind=kind_phys) :: cd
1738 real(kind=kind_phys) :: zblk,tautem
1739 real(kind=kind_phys) :: pe,ke
1740 real(kind=kind_phys) :: delx,dely
1741 real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4)
1742 real(kind=kind_phys) :: dxy(im),dxyp(im)
1743 real(kind=kind_phys) :: ol4p(4),olp(im),od(im)
1744 real(kind=kind_phys) :: taufb(im,km+1)
1745
1746 character(len=*), intent(out) :: errmsg
1747 integer, intent(out) :: errflg
1748
1749 integer :: udtend, vdtend, Tdtend
1750
1751 ! Calculate inverse of gravitational acceleration
1752 g_inv = 1./g
1753
1754 ! Initialize CCPP error handling variables
1755 errmsg = ''
1756 errflg = 0
1757
1758 ! Initialize local variables
1759 var_temp2 = 0.
1760 udtend = -1
1761 vdtend = -1
1762 tdtend = -1
1763
1764 if(ldiag3d) then
1765 udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd)
1766 vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd)
1767 tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd)
1768 endif
1769!
1770!---- constants
1771!
1772 rcl = 1.
1773 rcs = sqrt(rcl)
1774 cs = 1. / sqrt(rcl)
1775 csg = cs * g
1776 lcap = km
1777 lcapp1 = lcap + 1
1778 fdir = mdir / (2.0*pi)
1779 invgrcs = 1._kind_phys/g*rcs
1780 kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2
1781 denfac = 1.0
1782
1783 do i=1,im
1784 if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3
1785 xland(i)=1.0 !but land/water = (1/2) in this module
1786 else
1787 xland(i)=2.0
1788 endif
1789 rdxzb(i) = 0.0
1790 enddo
1791
1792!--- calculate scale-aware tapering factors
1793 do i=1,im
1794 if ( dx(i) .ge. dxmax_ls ) then
1795 ls_taper(i) = 1.
1796 else
1797 if ( dx(i) .le. dxmin_ls) then
1798 ls_taper(i) = 0.
1799 else
1800 ls_taper(i) = 0.5 * ( sin(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ &
1801 (dxmax_ls-dxmin_ls)) + 1. )
1802 endif
1803 endif
1804 enddo
1805
1806 ! Remove ss_tapering
1807 ss_taper(:) = 1.
1808
1809 ! SPP, if spp_gwd is 0, no perturbations are applied.
1810 if ( spp_gwd==1 ) then
1811 do i = its,im
1812 var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1)
1813 varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1)
1814 varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1)
1815 varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1)
1816 enddo
1817 else
1818 do i = its,im
1819 var_stoch(i) = var(i)
1820 varss_stoch(i) = varss(i)
1821 varmax_ss_stoch(i) = varmax_ss
1822 varmax_fd_stoch(i) = varmax_fd
1823 enddo
1824 endif
1825
1826 !--- calculate length of grid for flow-blocking drag
1827 !
1828 do i=1,im
1829 delx = dx(i)
1830 dely = dx(i)
1831 dxy4(i,1) = delx
1832 dxy4(i,2) = dely
1833 dxy4(i,3) = sqrt(delx*delx + dely*dely)
1834 dxy4(i,4) = dxy4(i,3)
1835 dxy4p(i,1) = dxy4(i,2)
1836 dxy4p(i,2) = dxy4(i,1)
1837 dxy4p(i,3) = dxy4(i,4)
1838 dxy4p(i,4) = dxy4(i,3)
1839 cleff(i) = psl_gwd_dx_factor*(delx+dely)*0.5
1840 cleff_ss(i) = 0.1 * max(dxmax_ss,dxy4(i,3))
1841 ! cleff_ss(i) = cleff(i) ! consider .....
1842 enddo
1843!
1844!-----initialize arrays
1845!
1846 dtaux = 0.0
1847 dtauy = 0.0
1848 do i = its,im
1849 klowtop(i) = 0
1850 kbl(i) = 0
1851 enddo
1852!
1853 do i = its,im
1854 xn(i) = 0.0
1855 yn(i) = 0.0
1856 ubar(i) = 0.0
1857 vbar(i) = 0.0
1858 roll(i) = 0.0
1859 taub(i) = 0.0
1860 oa(i) = 0.0
1861 ol(i) = 0.0
1862 oc(i) = 0.0
1863 oass(i) = 0.0
1864 olss(i) = 0.0
1865 ulow(i) = 0.0
1866 rstoch(i) = 0.0
1867 ldrag(i) = .false.
1868 icrilv(i) = .false.
1869 enddo
1870
1871 do k = kts,km
1872 do i = its,im
1873 usqj(i,k) = 0.0
1874 bnv2(i,k) = 0.0
1875 vtj(i,k) = 0.0
1876 vtk(i,k) = 0.0
1877 taup(i,k) = 0.0
1878 taud_ls(i,k) = 0.0
1879 taud_bl(i,k) = 0.0
1880 dtaux2d(i,k) = 0.0
1881 dtauy2d(i,k) = 0.0
1882 dtfac(i,k) = 1.0
1883 enddo
1884 enddo
1885!
1886 if ( ldiag_ugwp ) then
1887 do i = its,im
1888 dusfc_ls(i) = 0.0
1889 dvsfc_ls(i) = 0.0
1890 dusfc_bl(i) = 0.0
1891 dvsfc_bl(i) = 0.0
1892 dusfc_ss(i) = 0.0
1893 dvsfc_ss(i) = 0.0
1894 dusfc_fd(i) = 0.0
1895 dvsfc_fd(i) = 0.0
1896 enddo
1897 do k = kts,km
1898 do i = its,im
1899 dtaux2d_ls(i,k)= 0.0
1900 dtauy2d_ls(i,k)= 0.0
1901 dtaux2d_bl(i,k)= 0.0
1902 dtauy2d_bl(i,k)= 0.0
1903 dtaux2d_ss(i,k)= 0.0
1904 dtauy2d_ss(i,k)= 0.0
1905 dtaux2d_fd(i,k)= 0.0
1906 dtauy2d_fd(i,k)= 0.0
1907 enddo
1908 enddo
1909 endif
1910
1911 do i = its,im
1912 taup(i,km+1) = 0.0
1913 xlinv(i) = 1.0/xl
1914 dusfc(i) = 0.0
1915 dvsfc(i) = 0.0
1916 enddo
1917!
1918! initialize array for flow-blocking drag
1919!
1920 taufb(1:im,1:km+1) = 0.0
1921 hmax(1:im) = 0.0
1922 komax(1:im) = 0
1923 kbmax(1:im) = 0
1924 kblk(1:im) = 0
1925!
1926 rd_inv = 1./rd
1927 do k = kts,km
1928 do i = its,im
1929 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
1930 vtk(i,k) = vtj(i,k) / prslk(i,k)
1931 ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3
1932 enddo
1933 enddo
1934!
1935! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
1936!
1937 !zq=0.
1938 do k = kts,km
1939 do i = its,im
1940 !zq(i,k+1) = PHII(i,k+1)*g_inv
1941 !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv
1942 zl(i,k) = phil(i,k)*g_inv
1943 enddo
1944 enddo
1945!
1946! determine reference level: maximum of 2*var and pbl heights
1947!
1948 do i = its,im
1949 if(vtype(i)==15) then
1950 zlowtop(i) = 1.0 * var_stoch(i) !!! reduce drag over land ice
1951 else
1952 zlowtop(i) = 2.0 * var_stoch(i)
1953 endif
1954 enddo
1955!
1956 do i = its,im
1957 flag(i) = .true.
1958 enddo
1959!
1960 do k = kts+1,km
1961 do i = its,im
1962 if(flag(i).and.zl(i,k).ge.zlowtop(i)) then
1963 klowtop(i) = k+1
1964 flag(i) = .false.
1965 endif
1966 enddo
1967 enddo
1968!
1969! determine the maximum height level
1970! note taht elvmax and zl are the heights from the model surface whereas
1971! oro (mean orography) is the height from the sea level
1972!
1973 do i = its,im
1974 flag(i) = .true.
1975 enddo
1976!
1977 do k = kts+1,km
1978 do i = its,im
1979 if(flag(i).and.zl(i,k).ge.elvmax(i)) then
1980 komax(i) = k+1
1981 flag(i) = .false.
1982 endif
1983 enddo
1984 enddo
1985!
1986! determine the launching level in determining blocking layer
1987!
1988 do i = its,im
1989 flag(i) = .true.
1990 enddo
1991!
1992 do k = kts+1,km
1993 do i = its,im
1994 if(flag(i).and.zl(i,k).ge.elvmax(i)+zlowtop(i)) then
1995 kbmax(i) = k+1
1996 flag(i) = .false.
1997 endif
1998 enddo
1999 enddo
2000!
2001! determing the reference level for gwd and blockding...
2002!
2003 do i = its,im
2004 hmax(i) = max(elvmax(i),zlowtop(i))
2005 enddo
2006!
2007 do i = its,im
2008!!! kbl(i) = max(kpbl(i), klowtop(i)) ! do not use pbl height for the time being...
2009 kbl(i) = max(komax(i), klowtop(i))
2010 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
2011 enddo
2012!
2013! compute low level averages below reference level
2014!
2015 do i = its,im
2016 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
2017 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
2018 enddo
2019 do k = kts,kpblmax
2020 do i = its,im
2021 if (k.lt.kbl(i)) then
2022 rcsks = rcs * del(i,k) * delks(i)
2023 rdelks = del(i,k) * delks(i)
2024 ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
2025 vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
2026 roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
2027 endif
2028 enddo
2029 enddo
2030!
2031! figure out low-level horizontal wind direction
2032!
2033! nwd 1 2 3 4 5 6 7 8
2034! wd w s sw nw e n ne se
2035!
2036 do i = its,im
2037 wdir = atan2(ubar(i),vbar(i)) + pi
2038 idir = mod(nint(fdir*wdir),mdir) + 1
2039 nwd = nwdir(idir)
2040 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
2041 ol(i) = max(ol4(i,mod(nwd-1,4)+1),olmin)
2042 oc(i) = min(max(oc1(i),oc_min),oc_max)
2043! if (var(i).le.var_min) then
2044! oc(i) = max(oc(i)*var(i)/var_min,oc_min)
2045! endif
2046 ! Repeat for small-scale gwd
2047 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
2048 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
2049
2050!
2051!----- compute orographic width along (ol) and perpendicular (olp)
2052!----- the direction of wind
2053!
2054 ol4p(1) = ol4(i,2)
2055 ol4p(2) = ol4(i,1)
2056 ol4p(3) = ol4(i,4)
2057 ol4p(4) = ol4(i,3)
2058 olp(i) = max(ol4p(mod(nwd-1,4)+1),olmin)
2059!
2060!----- compute orographic direction (horizontal orographic aspect ratio)
2061!
2062 od(i) = olp(i)/ol(i)
2063 od(i) = min(od(i),odmax)
2064 od(i) = max(od(i),odmin)
2065!
2066!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
2067!
2068 dxy(i) = dxy4(i,mod(nwd-1,4)+1)
2069 dxyp(i) = dxy4p(i,mod(nwd-1,4)+1)
2070 enddo
2071!
2072! END INITIALIZATION; BEGIN GWD CALCULATIONS:
2073!
2074IF ( (do_gsl_drag_ls_bl).and. &
2075 ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then
2076
2077 g_cp = g/cp
2078 do i=its,im
2079
2080 if ( ls_taper(i).GT.1.e-02 ) then
2081
2082!
2083!--- saving richardson number in usqj for migwdi
2084!
2085 do k = kts,km-1
2086 ti = 2.0 / (t1(i,k)+t1(i,k+1))
2087 rdz = 1./(zl(i,k+1) - zl(i,k))
2088 tem1 = u1(i,k) - u1(i,k+1)
2089 tem2 = v1(i,k) - v1(i,k+1)
2090 dw2 = rcl*(tem1*tem1 + tem2*tem2)
2091 shr2 = max(dw2,dw2min) * rdz * rdz
2092 bvf2 = g*(g_cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
2093 usqj(i,k) = max(bvf2/shr2,rimin)
2094 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
2095 bnv2(i,k) = max( bnv2(i,k), bnv2min )
2096 enddo
2097!
2098!----compute the "low level" or 1/3 wind magnitude (m/s)
2099!
2100 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
2101 rulow(i) = 1./ulow(i)
2102!
2103 do k = kts,km-1
2104 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
2105 + (v1(i,k)+v1(i,k+1)) * vbar(i))
2106 velco(i,k) = velco(i,k) * rulow(i)
2107 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
2108 velco(i,k) = veleps
2109 endif
2110 enddo
2111!
2112! no drag when sub-oro is too small..
2113!
2114 ldrag(i) = hmax(i).le.hmt_min
2115!
2116! no drag when critical level in the base layer
2117!
2118 ldrag(i) = ldrag(i).or. velco(i,1).le.0.
2119!
2120! no drag when velco.lt.0
2121!
2122 do k = kpblmin,kpblmax
2123 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
2124 enddo
2125!
2126! no drag when bnv2.lt.0
2127!
2128 do k = kts,kpblmax
2129 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
2130 enddo
2131!
2132!-----the low level weighted average ri is stored in usqj(1,1; im)
2133!-----the low level weighted average n**2 is stored in bnv2(1,1; im)
2134!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
2135!---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
2136!
2137 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
2138 bnv2(i,1) = wtkbj * bnv2(i,1)
2139 usqj(i,1) = wtkbj * usqj(i,1)
2140!
2141 do k = kpblmin,kpblmax
2142 if (k .lt. kbl(i)) then
2143 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
2144 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
2145 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
2146 endif
2147 enddo
2148!
2149 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
2150 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
2151 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
2152 ldrag(i) = ldrag(i) .or. xland(i) .gt. 1.5
2153!
2154! set all ri low level values to the low level value
2155!
2156 do k = kpblmin,kpblmax
2157 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
2158 enddo
2159!
2160 if (.not.ldrag(i)) then
2161 bnv(i) = sqrt( bnv2(i,1) )
2162 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
2163 fr(i) = min(fr(i),frmax)
2164 xn(i) = ubar(i) * rulow(i)
2165 yn(i) = vbar(i) * rulow(i)
2166 endif
2167!
2168! compute the base level stress and store it in taub
2169! calculate enhancement factor, number of mountains & aspect
2170! ratio const. use simplified relationship between standard
2171! deviation & critical hgt
2172
2173 if (.not. ldrag(i)) then
2174 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
2175 efact = min( max(efact,efmin), efmax )
2176 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
2177 xlinv(i) = coefm(i) / cleff(i)
2178 tem = fr(i) * fr(i) * oc(i)
2179 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
2180 if ( gwd_opt_ls .NE. 0 ) then
2181 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
2182 * ulow(i) * gfobnv * efact
2183 else ! We've gotten what we need for the blocking scheme
2184 taub(i) = 0.0
2185 end if
2186 else
2187 taub(i) = 0.0
2188 xn(i) = 0.0
2189 yn(i) = 0.0
2190 endif
2191
2192 endif ! (ls_taper(i).GT.1.E-02)
2193
2194 enddo ! do i=its,im
2195
2196ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1))
2197
2198!=========================================================
2199! add small-scale wavedrag for stable boundary layer
2200!=========================================================
2201 xnbv=0.
2202 tauwavex0=0.
2203 tauwavey0=0.
2204 density=1.2
2205 utendwave=0.
2206 vtendwave=0.
2207!
2208IF ( do_gsl_drag_ss ) THEN
2209
2210 do i=its,im
2211
2212 if ( ss_taper(i).GT.1.e-02 ) then
2213 !
2214 ! calculating potential temperature
2215 !
2216 do k = kts,km
2217 thx(i,k) = t1(i,k)/prslk(i,k)
2218 enddo
2219 !
2220 do k = kts,km
2221 tvcon = (1.+fv*q1(i,k))
2222 thvx(i,k) = thx(i,k)*tvcon
2223 enddo
2224
2225 hpbl2 = hpbl(i)+10.
2226 kpbl2 = kpbl(i)
2227 !kvar = MIN(kpbl, k-level of var)
2228 kvar = 1
2229 do k=kts+1,max(kpbl(i),kts+1)
2230! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then
2231 IF (zl(i,k)>300.) then
2232 kpbl2 = k
2233 IF (k == kpbl(i)) then
2234 hpbl2 = hpbl(i)+10.
2235 ELSE
2236 hpbl2 = zl(i,k)+10.
2237 ENDIF
2238 exit
2239 ENDIF
2240 enddo
2241 if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
2242 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
2243 coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
2244 xlinv(i) = coefm_ss(i) / cleff_ss(i)
2245 !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts)))
2246 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
2247 !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i))
2248 xnbv=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
2249!
2250 !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then
2251 if(abs(xnbv/u1(i,kpbl2)).gt.xlinv(i))then
2252 !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i))
2253 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2)
2254 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
2255 ! Remove limit on varss_stoch
2256 var_temp = varss_stoch(i)
2257 ! Note: This is a semi-implicit treatment of the time differencing
2258 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
2259 tauwavex0=var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
2260 tauwavex0=tauwavex0*ss_taper(i)
2261 else
2262 tauwavex0=0.
2263 endif
2264!
2265 !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then
2266 if(abs(xnbv/v1(i,kpbl2)).gt.xlinv(i))then
2267 !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i))
2268 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2)
2269 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
2270 ! Remove limit on varss_stoch
2271 var_temp = varss_stoch(i)
2272 ! Note: This is a semi-implicit treatment of the time differencing
2273 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
2274 tauwavey0=var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
2275 tauwavey0=tauwavey0*ss_taper(i)
2276 else
2277 tauwavey0=0.
2278 endif
2279
2280 do k=kts,kpbl(i) !MIN(kpbl2+1,km-1)
2281!original
2282 !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
2283 !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
2284!new
2285 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2286 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2287!mod-to be used in HRRRv3/RAPv4
2288 !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2
2289 !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2
2290 enddo
2291 endif
2292 endif
2293
2294 do k = kts,km
2295 dudt(i,k) = dudt(i,k) + utendwave(i,k)
2296 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
2297 dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k)
2298 dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k)
2299 enddo
2300 if(udtend>0) then
2301 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim
2302 endif
2303 if(vdtend>0) then
2304 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim
2305 endif
2306 if ( ldiag_ugwp ) then
2307 do k = kts,km
2308 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
2309 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
2310 dtaux2d_ss(i,k) = utendwave(i,k)
2311 dtauy2d_ss(i,k) = vtendwave(i,k)
2312 enddo
2313 endif
2314
2315 endif ! if (ss_taper(i).GT.1.E-02)
2316
2317 enddo ! i=its,im
2318
2319ENDIF ! if (do_gsl_drag_ss)
2320
2321!================================================================
2322! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16):
2323!================================================================
2324IF ( do_gsl_drag_tofd ) THEN
2325
2326 do i=its,im
2327
2328 if ( ss_taper(i).GT.1.e-02 ) then
2329
2330 utendform=0.
2331 vtendform=0.
2332
2333 IF ((xland(i)-1.5) .le. 0.) then
2334 !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
2335 ! Remove limit on varss_stoch
2336 var_temp = varss_stoch(i)
2337 !var_temp = MIN(var_temp, 250.)
2338 a1=0.00026615161*var_temp**2
2339! a1=0.00026615161*MIN(varss(i),varmax)**2
2340! a1=0.00026615161*(0.5*varss(i))**2
2341 ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
2342 a2=a1*0.005363
2343 ! Beljaars H_efold
2344 h_efold = 1500.
2345 DO k=kts,km
2346 wsp=sqrt(u1(i,k)**2 + v1(i,k)**2)
2347 ! Note: In Beljaars et al. (2004):
2348 ! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
2349 ! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys
2350 var_temp = alpha_fd*coeff_fd*exp(-(zl(i,k)/h_efold)**1.5)*a2* &
2351 zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
2352 ! Note: This is a semi-implicit treatment of the time differencing
2353 ! per Beljaars et al. (2004, QJRMS)
2354 utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp)
2355 vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp)
2356 !IF(zl(i,k) > 4000.) exit
2357 ENDDO
2358 ENDIF
2359
2360 do k = kts,km
2361 dudt(i,k) = dudt(i,k) + utendform(i,k)
2362 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
2363 dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k)
2364 dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k)
2365 enddo
2366 if(udtend>0) then
2367 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim
2368 endif
2369 if(vdtend>0) then
2370 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim
2371 endif
2372 if ( ldiag_ugwp ) then
2373 do k = kts,km
2374 dtaux2d_fd(i,k) = utendform(i,k)
2375 dtauy2d_fd(i,k) = vtendform(i,k)
2376 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
2377 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
2378 enddo
2379 endif
2380
2381 endif ! if (ss_taper(i).GT.1.E-02)
2382
2383 enddo ! i=its,im
2384
2385ENDIF ! if (do_gsl_drag_tofd)
2386!=======================================================
2387! More for the large-scale gwd component
2388IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN
2389
2390 do i=its,im
2391
2392 if ( ls_taper(i).GT.1.e-02 ) then
2393
2394!
2395! now compute vertical structure of the stress.
2396 do k = kts,kpblmax
2397 if (k .le. kbl(i)) taup(i,k) = taub(i)
2398 enddo
2399!
2400 do k = kpblmin, km-1 ! vertical level k loop!
2401 kp1 = k + 1
2402!
2403! unstablelayer if ri < ric
2404! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
2405! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
2406!
2407 if (k .ge. kbl(i)) then
2408 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
2409 .or. (velco(i,k) .le. 0.0)
2410 brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
2411 brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
2412 endif
2413!
2414 if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
2415 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
2416 temv = 1.0 / velco(i,k)
2417 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* &
2418 velco(i,k)*0.5
2419 hd = sqrt(taup(i,k) / tem1)
2420 fro = brvf(i) * hd * temv
2421!
2422! rim is the minimum-richardson number by shutts (1985)
2423 tem2 = sqrt(usqj(i,k))
2424 tem = 1. + tem2 * fro
2425 rim = usqj(i,k) * (1.-fro) / (tem * tem)
2426!
2427! check stability to employ the 'saturation hypothesis'
2428! of lindzen (1981) except at tropospheric downstream regions
2429!
2430 if (rim .le. ric) then ! saturation hypothesis!
2431 temc = 2.0 + 1.0 / tem2
2432 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
2433 taup(i,kp1) = tem1 * hd * hd
2434 else ! no wavebreaking!
2435 taup(i,kp1) = taup(i,k)
2436 endif
2437 endif
2438 endif
2439 enddo
2440!
2441 if(lcap.lt.km) then
2442 do klcap = lcapp1,km
2443 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
2444 enddo
2445 endif
2446
2447 endif ! if ( ls_taper(i).GT.1.E-02 )
2448
2449 enddo ! do i=its,im
2450
2451ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1)
2452!===============================================================
2453!COMPUTE BLOCKING COMPONENT
2454!===============================================================
2455IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN
2456 do i = its,im
2457 flag(i) = .true.
2458 enddo
2459
2460 do i=its,im
2461
2462 if ( ls_taper(i).GT.1.e-02 ) then
2463
2464 if (.not.ldrag(i)) then
2465!
2466!------- determine the height of flow-blocking layer
2467!
2468 pe = 0.0
2469 ke = 0.0
2470 do k = km, kpblmin, -1
2471 if(flag(i).and. k.le.kbmax(i)) then
2472 pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* &
2473 del(i,k)*g_inv/ro(i,k)
2474 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
2475!
2476!---------- apply flow-blocking drag when pe >= ke
2477!
2478 if(pe.ge.ke.and.zl(i,k).le.hmax(i)) then
2479 kblk(i)= k
2480 zblk = zl(i,k)
2481 rdxzb(i) = real(k,kind=kind_phys)
2482 flag(i) = .false.
2483 endif
2484 endif
2485 enddo
2486 if(.not.flag(i)) then
2487!
2488!--------- compute flow-blocking stress
2489!
2490 cd = max(2.0-1.0/od(i),cdmin)
2491 taufb(i,kts) = 0.5 * roll(i) * coefm(i) / &
2492 max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * &
2493 olp(i) * zblk * ulow(i)**2
2494 tautem = taufb(i,kts)/float(kblk(i)-kts)
2495 do k = kts+1, kpblmax
2496 if (k .le. kblk(i)) taufb(i,k) = taufb(i,k-1) - tautem
2497 enddo
2498!
2499! reset gwd stress below blocking layer
2500!
2501 do k = kts,kpblmax
2502 if (k .le. kblk(i)) taup(i,k) = taup(i,kblk(i))
2503 enddo
2504! if(kblk(i).gt.5) print *,' gwd kbl komax kbmax kblk ',kbl(i),komax(i),kbmax(i),kblk(i)
2505! if(kblk(i).gt.5) print *,' gwd elvmax zlowtop zblk ',elvmax(i),zlowtop(i),zl(i,kblk(i))
2506 endif
2507
2508 endif ! if (.not.ldrag(i))
2509
2510 endif ! if ( ls_taper(i).GT.1.E-02 )
2511
2512 enddo ! do i=its,im
2513
2514ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) )
2515!===========================================================
2516IF ( (do_gsl_drag_ls_bl) .and. &
2517 (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN
2518
2519 do i=its,im
2520
2521 if ( ls_taper(i) .GT. 1.e-02 ) then
2522
2523!
2524! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
2525!
2526 do k = kts,km
2527 taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
2528 taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
2529 enddo
2530!
2531! limit de-acceleration (momentum deposition ) at top to 1/2 value
2532! the idea is some stuff must go out the 'top'
2533 do klcap = lcap,km
2534 taud_ls(i,klcap) = taud_ls(i,klcap) * factop
2535 taud_bl(i,klcap) = taud_bl(i,klcap) * factop
2536 enddo
2537!
2538! if the gravity wave drag would force a critical line
2539! in the lower ksmm1 layers during the next deltim timestep,
2540! then only apply drag until that critical line is reached.
2541!
2542 do k = kts,kpblmax-1
2543 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) then
2544 dtfac(i,k) = min(dtfac(i,k),abs(velco(i,k) &
2545 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2546 endif
2547 enddo
2548! apply limiter to mesosphere drag, reduce the drag by density factor 10-3
2549! prevent wind reversal...
2550!
2551 do k = kpblmax,km-1
2552 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff) then
2553 denfac = min(ro(i,k)/pcutoff_den,1.)
2554 dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k) &
2555 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2556 endif
2557 enddo
2558!
2559 do k = kts,km
2560 taud_ls(i,k) = taud_ls(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2561 taud_bl(i,k) = taud_bl(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2562 dtaux = taud_ls(i,k) * xn(i)
2563 dtauy = taud_ls(i,k) * yn(i)
2564 dtauxb = taud_bl(i,k) * xn(i)
2565 dtauyb = taud_bl(i,k) * yn(i)
2566
2567 !add blocking and large-scale contributions to tendencies
2568 dudt(i,k) = dtaux + dtauxb + dudt(i,k)
2569 dvdt(i,k) = dtauy + dtauyb + dvdt(i,k)
2570
2571 if ( gsd_diss_ht_opt .EQ. 1 ) then
2572 ! Calculate dissipation heating
2573 ! Initial kinetic energy (at t0-dt)
2574 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
2575 ! Kinetic energy after wave-breaking/flow-blocking
2576 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
2577 (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
2578 ! Modify theta tendency
2579 dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim
2580 if ( tdtend>0 ) then
2581 dtend(i,k,tdtend) = dtend(i,k,tdtend) + max((eng0-eng1),0.0)/cp
2582 endif
2583 endif
2584
2585 dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + &
2586 taud_bl(i,k)*xn(i)*del(i,k)
2587 dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + &
2588 taud_bl(i,k)*yn(i)*del(i,k)
2589 if(udtend>0) then
2590 dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * &
2591 xn(i) + taud_bl(i,k) * xn(i)) * deltim
2592 endif
2593 if(vdtend>0) then
2594 dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * &
2595 yn(i) + taud_bl(i,k) * yn(i)) * deltim
2596 endif
2597
2598 enddo
2599
2600 ! Finalize dusfc and dvsfc diagnostics
2601 dusfc(i) = -(invgrcs) * dusfc(i)
2602 dvsfc(i) = -(invgrcs) * dvsfc(i)
2603
2604 if ( ldiag_ugwp ) then
2605 do k = kts,km
2606 dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i)
2607 dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i)
2608 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
2609 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
2610 dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k)
2611 dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k)
2612 dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k)
2613 dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k)
2614 enddo
2615 endif
2616
2617 endif ! if ( ls_taper(i) .GT. 1.E-02 )
2618
2619 enddo ! do i=its,im
2620
2621ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1)
2622
2623if ( ldiag_ugwp ) then
2624 ! Finalize dusfc and dvsfc diagnostics
2625 do i = its,im
2626 dusfc_ls(i) = -(invgrcs) * dusfc_ls(i)
2627 dvsfc_ls(i) = -(invgrcs) * dvsfc_ls(i)
2628 dusfc_bl(i) = -(invgrcs) * dusfc_bl(i)
2629 dvsfc_bl(i) = -(invgrcs) * dvsfc_bl(i)
2630 dusfc_ss(i) = -(invgrcs) * dusfc_ss(i)
2631 dvsfc_ss(i) = -(invgrcs) * dvsfc_ss(i)
2632 dusfc_fd(i) = -(invgrcs) * dusfc_fd(i)
2633 dvsfc_fd(i) = -(invgrcs) * dvsfc_fd(i)
2634 enddo
2635endif
2636!
2637 return
2638 end subroutine drag_suite_psl
2639!
2641
2642 end module drag_suite
This module contains the orographic gravity wave drag scheme.
Definition drag_suite.F90:6