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