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