CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
drag_suite.F90
1
4
6 module drag_suite
7
8 contains
9
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), optional :: 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), optional :: 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) :: 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 hmax(1:im) = 0.0
1920 komax(1:im) = 0
1921 kbmax(1:im) = 0
1922 kblk(1:im) = 0
1923!
1924 rd_inv = 1./rd
1925 do k = kts,km
1926 do i = its,im
1927 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
1928 vtk(i,k) = vtj(i,k) / prslk(i,k)
1929 ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k) ! density kg/m**3
1930 enddo
1931 enddo
1932!
1933! calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
1934!
1935 !zq=0.
1936 do k = kts,km
1937 do i = its,im
1938 !zq(i,k+1) = PHII(i,k+1)*g_inv
1939 !dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv
1940 zl(i,k) = phil(i,k)*g_inv
1941 enddo
1942 enddo
1943!
1944! determine reference level: maximum of 2*var and pbl heights
1945!
1946 do i = its,im
1947 if(vtype(i)==15) then
1948 zlowtop(i) = 1.0 * var_stoch(i) !!! reduce drag over land ice
1949 else
1950 zlowtop(i) = 2.0 * var_stoch(i)
1951 endif
1952 enddo
1953!
1954 do i = its,im
1955 flag(i) = .true.
1956 enddo
1957!
1958 do k = kts+1,km
1959 do i = its,im
1960 if(flag(i).and.zl(i,k).ge.zlowtop(i)) then
1961 klowtop(i) = k+1
1962 flag(i) = .false.
1963 endif
1964 enddo
1965 enddo
1966!
1967! determine the maximum height level
1968! note taht elvmax and zl are the heights from the model surface whereas
1969! oro (mean orography) is the height from the sea level
1970!
1971 do i = its,im
1972 flag(i) = .true.
1973 enddo
1974!
1975 do k = kts+1,km
1976 do i = its,im
1977 if(flag(i).and.zl(i,k).ge.elvmax(i)) then
1978 komax(i) = k+1
1979 flag(i) = .false.
1980 endif
1981 enddo
1982 enddo
1983!
1984! determine the launching level in determining blocking layer
1985!
1986 do i = its,im
1987 flag(i) = .true.
1988 enddo
1989!
1990 do k = kts+1,km
1991 do i = its,im
1992 if(flag(i).and.zl(i,k).ge.elvmax(i)+zlowtop(i)) then
1993 kbmax(i) = k+1
1994 flag(i) = .false.
1995 endif
1996 enddo
1997 enddo
1998!
1999! determing the reference level for gwd and blockding...
2000!
2001 do i = its,im
2002 hmax(i) = max(elvmax(i),zlowtop(i))
2003 enddo
2004!
2005 do i = its,im
2006!!! kbl(i) = max(kpbl(i), klowtop(i)) ! do not use pbl height for the time being...
2007 kbl(i) = max(komax(i), klowtop(i))
2008 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
2009 enddo
2010!
2011! compute low level averages below reference level
2012!
2013 do i = its,im
2014 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
2015 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
2016 enddo
2017 do k = kts,kpblmax
2018 do i = its,im
2019 if (k.lt.kbl(i)) then
2020 rcsks = rcs * del(i,k) * delks(i)
2021 rdelks = del(i,k) * delks(i)
2022 ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
2023 vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
2024 roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
2025 endif
2026 enddo
2027 enddo
2028!
2029! figure out low-level horizontal wind direction
2030!
2031! nwd 1 2 3 4 5 6 7 8
2032! wd w s sw nw e n ne se
2033!
2034 do i = its,im
2035 wdir = atan2(ubar(i),vbar(i)) + pi
2036 idir = mod(nint(fdir*wdir),mdir) + 1
2037 nwd = nwdir(idir)
2038 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
2039 ol(i) = max(ol4(i,mod(nwd-1,4)+1),olmin)
2040 oc(i) = min(max(oc1(i),oc_min),oc_max)
2041! if (var(i).le.var_min) then
2042! oc(i) = max(oc(i)*var(i)/var_min,oc_min)
2043! endif
2044 ! Repeat for small-scale gwd
2045 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
2046 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
2047
2048!
2049!----- compute orographic width along (ol) and perpendicular (olp)
2050!----- the direction of wind
2051!
2052 ol4p(1) = ol4(i,2)
2053 ol4p(2) = ol4(i,1)
2054 ol4p(3) = ol4(i,4)
2055 ol4p(4) = ol4(i,3)
2056 olp(i) = max(ol4p(mod(nwd-1,4)+1),olmin)
2057!
2058!----- compute orographic direction (horizontal orographic aspect ratio)
2059!
2060 od(i) = olp(i)/ol(i)
2061 od(i) = min(od(i),odmax)
2062 od(i) = max(od(i),odmin)
2063!
2064!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
2065!
2066 dxy(i) = dxy4(i,mod(nwd-1,4)+1)
2067 dxyp(i) = dxy4p(i,mod(nwd-1,4)+1)
2068 enddo
2069!
2070! END INITIALIZATION; BEGIN GWD CALCULATIONS:
2071!
2072IF ( (do_gsl_drag_ls_bl).and. &
2073 ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) ) then
2074
2075 g_cp = g/cp
2076 do i=its,im
2077
2078 if ( ls_taper(i).GT.1.e-02 ) then
2079
2080!
2081!--- saving richardson number in usqj for migwdi
2082!
2083 do k = kts,km-1
2084 ti = 2.0 / (t1(i,k)+t1(i,k+1))
2085 rdz = 1./(zl(i,k+1) - zl(i,k))
2086 tem1 = u1(i,k) - u1(i,k+1)
2087 tem2 = v1(i,k) - v1(i,k+1)
2088 dw2 = rcl*(tem1*tem1 + tem2*tem2)
2089 shr2 = max(dw2,dw2min) * rdz * rdz
2090 bvf2 = g*(g_cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
2091 usqj(i,k) = max(bvf2/shr2,rimin)
2092 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
2093 bnv2(i,k) = max( bnv2(i,k), bnv2min )
2094 enddo
2095!
2096!----compute the "low level" or 1/3 wind magnitude (m/s)
2097!
2098 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
2099 rulow(i) = 1./ulow(i)
2100!
2101 do k = kts,km-1
2102 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
2103 + (v1(i,k)+v1(i,k+1)) * vbar(i))
2104 velco(i,k) = velco(i,k) * rulow(i)
2105 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
2106 velco(i,k) = veleps
2107 endif
2108 enddo
2109!
2110! no drag when sub-oro is too small..
2111!
2112 ldrag(i) = hmax(i).le.hmt_min
2113!
2114! no drag when critical level in the base layer
2115!
2116 ldrag(i) = ldrag(i).or. velco(i,1).le.0.
2117!
2118! no drag when velco.lt.0
2119!
2120 do k = kpblmin,kpblmax
2121 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
2122 enddo
2123!
2124! no drag when bnv2.lt.0
2125!
2126 do k = kts,kpblmax
2127 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
2128 enddo
2129!
2130!-----the low level weighted average ri is stored in usqj(1,1; im)
2131!-----the low level weighted average n**2 is stored in bnv2(1,1; im)
2132!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
2133!---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
2134!
2135 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
2136 bnv2(i,1) = wtkbj * bnv2(i,1)
2137 usqj(i,1) = wtkbj * usqj(i,1)
2138!
2139 do k = kpblmin,kpblmax
2140 if (k .lt. kbl(i)) then
2141 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
2142 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
2143 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
2144 endif
2145 enddo
2146!
2147 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
2148 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
2149 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
2150 ldrag(i) = ldrag(i) .or. xland(i) .gt. 1.5
2151!
2152! set all ri low level values to the low level value
2153!
2154 do k = kpblmin,kpblmax
2155 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
2156 enddo
2157!
2158 if (.not.ldrag(i)) then
2159 bnv(i) = sqrt( bnv2(i,1) )
2160 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
2161 fr(i) = min(fr(i),frmax)
2162 xn(i) = ubar(i) * rulow(i)
2163 yn(i) = vbar(i) * rulow(i)
2164 endif
2165!
2166! compute the base level stress and store it in taub
2167! calculate enhancement factor, number of mountains & aspect
2168! ratio const. use simplified relationship between standard
2169! deviation & critical hgt
2170
2171 if (.not. ldrag(i)) then
2172 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
2173 efact = min( max(efact,efmin), efmax )
2174 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
2175 xlinv(i) = coefm(i) / cleff(i)
2176 tem = fr(i) * fr(i) * oc(i)
2177 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
2178 if ( gwd_opt_ls .NE. 0 ) then
2179 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
2180 * ulow(i) * gfobnv * efact
2181 else ! We've gotten what we need for the blocking scheme
2182 taub(i) = 0.0
2183 end if
2184 else
2185 taub(i) = 0.0
2186 xn(i) = 0.0
2187 yn(i) = 0.0
2188 endif
2189
2190 endif ! (ls_taper(i).GT.1.E-02)
2191
2192 enddo ! do i=its,im
2193
2194ENDIF ! (do_gsl_drag_ls_bl).and.((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1))
2195
2196!=========================================================
2197! add small-scale wavedrag for stable boundary layer
2198!=========================================================
2199 xnbv=0.
2200 tauwavex0=0.
2201 tauwavey0=0.
2202 density=1.2
2203 utendwave=0.
2204 vtendwave=0.
2205!
2206IF ( do_gsl_drag_ss ) THEN
2207
2208 do i=its,im
2209
2210 if ( ss_taper(i).GT.1.e-02 ) then
2211 !
2212 ! calculating potential temperature
2213 !
2214 do k = kts,km
2215 thx(i,k) = t1(i,k)/prslk(i,k)
2216 enddo
2217 !
2218 do k = kts,km
2219 tvcon = (1.+fv*q1(i,k))
2220 thvx(i,k) = thx(i,k)*tvcon
2221 enddo
2222
2223 hpbl2 = hpbl(i)+10.
2224 kpbl2 = kpbl(i)
2225 !kvar = MIN(kpbl, k-level of var)
2226 kvar = 1
2227 do k=kts+1,max(kpbl(i),kts+1)
2228! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then
2229 IF (zl(i,k)>300.) then
2230 kpbl2 = k
2231 IF (k == kpbl(i)) then
2232 hpbl2 = hpbl(i)+10.
2233 ELSE
2234 hpbl2 = zl(i,k)+10.
2235 ENDIF
2236 exit
2237 ENDIF
2238 enddo
2239 if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
2240 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
2241 coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
2242 xlinv(i) = coefm_ss(i) / cleff_ss(i)
2243 !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts)))
2244 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
2245 !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i))
2246 xnbv=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
2247!
2248 !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then
2249 if(abs(xnbv/u1(i,kpbl2)).gt.xlinv(i))then
2250 !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i))
2251 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2)
2252 !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
2253 ! Remove limit on varss_stoch
2254 var_temp = varss_stoch(i)
2255 ! Note: This is a semi-implicit treatment of the time differencing
2256 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
2257 tauwavex0=var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
2258 tauwavex0=tauwavex0*ss_taper(i)
2259 else
2260 tauwavex0=0.
2261 endif
2262!
2263 !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then
2264 if(abs(xnbv/v1(i,kpbl2)).gt.xlinv(i))then
2265 !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i))
2266 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2)
2267 !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
2268 ! Remove limit on varss_stoch
2269 var_temp = varss_stoch(i)
2270 ! Note: This is a semi-implicit treatment of the time differencing
2271 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
2272 tauwavey0=var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
2273 tauwavey0=tauwavey0*ss_taper(i)
2274 else
2275 tauwavey0=0.
2276 endif
2277
2278 do k=kts,kpbl(i) !MIN(kpbl2+1,km-1)
2279!original
2280 !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
2281 !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
2282!new
2283 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2284 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2285!mod-to be used in HRRRv3/RAPv4
2286 !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2
2287 !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2
2288 enddo
2289 endif
2290 endif
2291
2292 do k = kts,km
2293 dudt(i,k) = dudt(i,k) + utendwave(i,k)
2294 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
2295 dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k)
2296 dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k)
2297 enddo
2298 if(udtend>0) then
2299 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim
2300 endif
2301 if(vdtend>0) then
2302 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim
2303 endif
2304 if ( ldiag_ugwp ) then
2305 do k = kts,km
2306 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
2307 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
2308 dtaux2d_ss(i,k) = utendwave(i,k)
2309 dtauy2d_ss(i,k) = vtendwave(i,k)
2310 enddo
2311 endif
2312
2313 endif ! if (ss_taper(i).GT.1.E-02)
2314
2315 enddo ! i=its,im
2316
2317ENDIF ! if (do_gsl_drag_ss)
2318
2319!================================================================
2320! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16):
2321!================================================================
2322IF ( do_gsl_drag_tofd ) THEN
2323
2324 do i=its,im
2325
2326 if ( ss_taper(i).GT.1.e-02 ) then
2327
2328 utendform=0.
2329 vtendform=0.
2330
2331 IF ((xland(i)-1.5) .le. 0.) then
2332 !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
2333 ! Remove limit on varss_stoch
2334 var_temp = varss_stoch(i)
2335 !var_temp = MIN(var_temp, 250.)
2336 a1=0.00026615161*var_temp**2
2337! a1=0.00026615161*MIN(varss(i),varmax)**2
2338! a1=0.00026615161*(0.5*varss(i))**2
2339 ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
2340 a2=a1*0.005363
2341 ! Beljaars H_efold
2342 h_efold = 1500.
2343 DO k=kts,km
2344 wsp=sqrt(u1(i,k)**2 + v1(i,k)**2)
2345 ! Note: In Beljaars et al. (2004):
2346 ! alpha_fd*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
2347 ! lump beta*Cmd*Ccorr*2.109 into 1.*0.005*0.6*2.109 = coeff_fd ~ 6.325e-3_kind_phys
2348 var_temp = alpha_fd*coeff_fd*exp(-(zl(i,k)/h_efold)**1.5)*a2* &
2349 zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
2350 ! Note: This is a semi-implicit treatment of the time differencing
2351 ! per Beljaars et al. (2004, QJRMS)
2352 utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp)
2353 vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp)
2354 !IF(zl(i,k) > 4000.) exit
2355 ENDDO
2356 ENDIF
2357
2358 do k = kts,km
2359 dudt(i,k) = dudt(i,k) + utendform(i,k)
2360 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
2361 dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k)
2362 dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k)
2363 enddo
2364 if(udtend>0) then
2365 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim
2366 endif
2367 if(vdtend>0) then
2368 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim
2369 endif
2370 if ( ldiag_ugwp ) then
2371 do k = kts,km
2372 dtaux2d_fd(i,k) = utendform(i,k)
2373 dtauy2d_fd(i,k) = vtendform(i,k)
2374 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
2375 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
2376 enddo
2377 endif
2378
2379 endif ! if (ss_taper(i).GT.1.E-02)
2380
2381 enddo ! i=its,im
2382
2383ENDIF ! if (do_gsl_drag_tofd)
2384!=======================================================
2385! More for the large-scale gwd component
2386IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) ) THEN
2387
2388 do i=its,im
2389
2390 if ( ls_taper(i).GT.1.e-02 ) then
2391
2392!
2393! now compute vertical structure of the stress.
2394 do k = kts,kpblmax
2395 if (k .le. kbl(i)) taup(i,k) = taub(i)
2396 enddo
2397!
2398 do k = kpblmin, km-1 ! vertical level k loop!
2399 kp1 = k + 1
2400!
2401! unstablelayer if ri < ric
2402! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
2403! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
2404!
2405 if (k .ge. kbl(i)) then
2406 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
2407 .or. (velco(i,k) .le. 0.0)
2408 brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
2409 brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
2410 endif
2411!
2412 if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
2413 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
2414 temv = 1.0 / velco(i,k)
2415 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* &
2416 velco(i,k)*0.5
2417 hd = sqrt(taup(i,k) / tem1)
2418 fro = brvf(i) * hd * temv
2419!
2420! rim is the minimum-richardson number by shutts (1985)
2421 tem2 = sqrt(usqj(i,k))
2422 tem = 1. + tem2 * fro
2423 rim = usqj(i,k) * (1.-fro) / (tem * tem)
2424!
2425! check stability to employ the 'saturation hypothesis'
2426! of lindzen (1981) except at tropospheric downstream regions
2427!
2428 if (rim .le. ric) then ! saturation hypothesis!
2429 temc = 2.0 + 1.0 / tem2
2430 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
2431 taup(i,kp1) = tem1 * hd * hd
2432 else ! no wavebreaking!
2433 taup(i,kp1) = taup(i,k)
2434 endif
2435 endif
2436 endif
2437 enddo
2438!
2439 if(lcap.lt.km) then
2440 do klcap = lcapp1,km
2441 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
2442 enddo
2443 endif
2444
2445 endif ! if ( ls_taper(i).GT.1.E-02 )
2446
2447 enddo ! do i=its,im
2448
2449ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1)
2450!===============================================================
2451!COMPUTE BLOCKING COMPONENT
2452!===============================================================
2453IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) ) THEN
2454 do i = its,im
2455 flag(i) = .true.
2456 enddo
2457
2458 do i=its,im
2459
2460 if ( ls_taper(i).GT.1.e-02 ) then
2461
2462 if (.not.ldrag(i)) then
2463!
2464!------- determine the height of flow-blocking layer
2465!
2466 pe = 0.0
2467 ke = 0.0
2468 do k = km, kpblmin, -1
2469 if(flag(i).and. k.le.kbmax(i)) then
2470 pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* &
2471 del(i,k)*g_inv/ro(i,k)
2472 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
2473!
2474!---------- apply flow-blocking drag when pe >= ke
2475!
2476 if(pe.ge.ke.and.zl(i,k).le.hmax(i)) then
2477 kblk(i)= k
2478 zblk = zl(i,k)
2479 rdxzb(i) = real(k,kind=kind_phys)
2480 flag(i) = .false.
2481 endif
2482 endif
2483 enddo
2484 if(.not.flag(i)) then
2485!
2486!--------- compute flow-blocking stress
2487!
2488 cd = max(2.0-1.0/od(i),cdmin)
2489 taufb(i,kts) = 0.5 * roll(i) * coefm(i) / &
2490 max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * &
2491 olp(i) * zblk * ulow(i)**2
2492 tautem = taufb(i,kts)/float(kblk(i)-kts)
2493 do k = kts+1, kpblmax
2494 if (k .le. kblk(i)) taufb(i,k) = taufb(i,k-1) - tautem
2495 enddo
2496!
2497! reset gwd stress below blocking layer
2498!
2499 do k = kts,kpblmax
2500 if (k .le. kblk(i)) taup(i,k) = taup(i,kblk(i))
2501 enddo
2502! if(kblk(i).gt.5) print *,' gwd kbl komax kbmax kblk ',kbl(i),komax(i),kbmax(i),kblk(i)
2503! if(kblk(i).gt.5) print *,' gwd elvmax zlowtop zblk ',elvmax(i),zlowtop(i),zl(i,kblk(i))
2504 endif
2505
2506 endif ! if (.not.ldrag(i))
2507
2508 endif ! if ( ls_taper(i).GT.1.E-02 )
2509
2510 enddo ! do i=its,im
2511
2512ENDIF ! IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) )
2513!===========================================================
2514IF ( (do_gsl_drag_ls_bl) .and. &
2515 (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) ) THEN
2516
2517 do i=its,im
2518
2519 if ( ls_taper(i) .GT. 1.e-02 ) then
2520
2521!
2522! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
2523!
2524 do k = kts,km
2525 taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
2526 taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
2527 enddo
2528!
2529! limit de-acceleration (momentum deposition ) at top to 1/2 value
2530! the idea is some stuff must go out the 'top'
2531 do klcap = lcap,km
2532 taud_ls(i,klcap) = taud_ls(i,klcap) * factop
2533 taud_bl(i,klcap) = taud_bl(i,klcap) * factop
2534 enddo
2535!
2536! if the gravity wave drag would force a critical line
2537! in the lower ksmm1 layers during the next deltim timestep,
2538! then only apply drag until that critical line is reached.
2539!
2540 do k = kts,kpblmax-1
2541 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.) then
2542 dtfac(i,k) = min(dtfac(i,k),abs(velco(i,k) &
2543 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2544 endif
2545 enddo
2546! apply limiter to mesosphere drag, reduce the drag by density factor 10-3
2547! prevent wind reversal...
2548!
2549 do k = kpblmax,km-1
2550 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff) then
2551 denfac = min(ro(i,k)/pcutoff_den,1.)
2552 dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k) &
2553 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2554 endif
2555 enddo
2556!
2557 do k = kts,km
2558 taud_ls(i,k) = taud_ls(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2559 taud_bl(i,k) = taud_bl(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2560 dtaux = taud_ls(i,k) * xn(i)
2561 dtauy = taud_ls(i,k) * yn(i)
2562 dtauxb = taud_bl(i,k) * xn(i)
2563 dtauyb = taud_bl(i,k) * yn(i)
2564
2565 !add blocking and large-scale contributions to tendencies
2566 dudt(i,k) = dtaux + dtauxb + dudt(i,k)
2567 dvdt(i,k) = dtauy + dtauyb + dvdt(i,k)
2568
2569 if ( gsd_diss_ht_opt .EQ. 1 ) then
2570 ! Calculate dissipation heating
2571 ! Initial kinetic energy (at t0-dt)
2572 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
2573 ! Kinetic energy after wave-breaking/flow-blocking
2574 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
2575 (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
2576 ! Modify theta tendency
2577 dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim
2578 if ( tdtend>0 ) then
2579 dtend(i,k,tdtend) = dtend(i,k,tdtend) + max((eng0-eng1),0.0)/cp
2580 endif
2581 endif
2582
2583 dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + &
2584 taud_bl(i,k)*xn(i)*del(i,k)
2585 dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + &
2586 taud_bl(i,k)*yn(i)*del(i,k)
2587 if(udtend>0) then
2588 dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * &
2589 xn(i) + taud_bl(i,k) * xn(i)) * deltim
2590 endif
2591 if(vdtend>0) then
2592 dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * &
2593 yn(i) + taud_bl(i,k) * yn(i)) * deltim
2594 endif
2595
2596 enddo
2597
2598 ! Finalize dusfc and dvsfc diagnostics
2599 dusfc(i) = -(invgrcs) * dusfc(i)
2600 dvsfc(i) = -(invgrcs) * dvsfc(i)
2601
2602 if ( ldiag_ugwp ) then
2603 do k = kts,km
2604 dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i)
2605 dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i)
2606 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
2607 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
2608 dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k)
2609 dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k)
2610 dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k)
2611 dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k)
2612 enddo
2613 endif
2614
2615 endif ! if ( ls_taper(i) .GT. 1.E-02 )
2616
2617 enddo ! do i=its,im
2618
2619ENDIF ! (do_gsl_drag_ls_bl).and.(gwd_opt_ls.EQ.1 .OR. gwd_opt_bl.EQ.1)
2620
2621if ( ldiag_ugwp ) then
2622 ! Finalize dusfc and dvsfc diagnostics
2623 do i = its,im
2624 dusfc_ls(i) = -(invgrcs) * dusfc_ls(i)
2625 dvsfc_ls(i) = -(invgrcs) * dvsfc_ls(i)
2626 dusfc_bl(i) = -(invgrcs) * dusfc_bl(i)
2627 dvsfc_bl(i) = -(invgrcs) * dvsfc_bl(i)
2628 dusfc_ss(i) = -(invgrcs) * dusfc_ss(i)
2629 dvsfc_ss(i) = -(invgrcs) * dvsfc_ss(i)
2630 dusfc_fd(i) = -(invgrcs) * dusfc_fd(i)
2631 dvsfc_fd(i) = -(invgrcs) * dvsfc_fd(i)
2632 enddo
2633endif
2634!
2635 return
2636 end subroutine drag_suite_psl
2637!
2639
2640 end module drag_suite
drag_suite
This module contains the orographic drag scheme.
Definition drag_suite.F90:6