207 subroutine drag_suite_run( &
208 & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, &
209 & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, &
211 & varss,oc1ss,oa4ss,ol4ss, &
212 & THETA,SIGMA,GAMMA,ELVMAX, &
213 & dtaux2d_ms,dtauy2d_ms,dtaux2d_bl,dtauy2d_bl, &
214 & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, &
216 & dusfc_ms,dvsfc_ms,dusfc_bl,dvsfc_bl, &
217 & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
219 & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
220 & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
221 & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
222 & dtend, dtidx, index_of_process_orographic_gwd, &
223 & index_of_temperature, index_of_x_wind, &
224 & index_of_y_wind, ldiag3d, ldiag_ugwp, ugwp_seq_update, &
225 & spp_wts_gwd, spp_gwd, errmsg, errflg)
319 USE machine ,
ONLY : kind_phys
323 integer,
intent(in) :: im, km, imx, kdt, ipr, me, master
324 integer,
intent(in) :: gwd_opt
325 logical,
intent(in) :: lprnt
326 integer,
intent(in) :: KPBL(:)
327 real(kind=kind_phys),
intent(in) :: deltim, g, cp, rd, rv, &
328 & cdmbgwd(:), alpha_fd
329 real(kind=kind_phys),
intent(inout),
optional :: dtend(:,:,:)
330 logical,
intent(in) :: ldiag3d
331 integer,
intent(in) :: dtidx(:,:), index_of_temperature, &
332 & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind
335 integer,
parameter :: ims=1, kms=1, its=1, kts=1
336 real(kind=kind_phys),
intent(in) :: fv, pi
337 real(kind=kind_phys) :: rcl, cdmb
338 real(kind=kind_phys) :: g_inv, rd_inv
340 real(kind=kind_phys),
intent(inout) :: &
341 & dudt(:,:),dvdt(:,:), &
343 real(kind=kind_phys),
intent(inout) :: rdxzb(:)
344 real(kind=kind_phys),
intent(in) :: &
347 & phii(:,:),prsl(:,:), &
348 & prslk(:,:),phil(:,:)
349 real(kind=kind_phys),
intent(in) :: prsi(:,:), &
351 real(kind=kind_phys),
intent(in) :: var(:),oc1(:), &
352 & oa4(:,:),ol4(:,:), &
354 real(kind=kind_phys),
intent(in),
optional :: varss(:),oc1ss(:), &
355 & oa4ss(:,:),ol4ss(:,:)
356 real(kind=kind_phys),
intent(in) :: theta(:),sigma(:), &
360 real(kind=kind_phys),
dimension(im,km) :: utendwave,vtendwave,thx,thvx
361 real(kind=kind_phys),
intent(in) :: br1(:), &
364 real(kind=kind_phys),
dimension(im) :: govrth,xland
366 real(kind=kind_phys) :: tauwavex0,tauwavey0, &
367 & xnbv,density,tvcon,hpbl2
368 integer :: kpbl2,kvar
370 real(kind=kind_phys),
dimension(im,km) :: zl
373 real(kind=kind_phys),
dimension(im) :: var_stoch, varss_stoch, &
375 real(kind=kind_phys),
intent(in),
optional :: spp_wts_gwd(:,:)
376 integer,
intent(in) :: spp_gwd
378 real(kind=kind_phys),
dimension(im) :: rstoch
381 real(kind=kind_phys),
intent(inout) :: &
384 real(kind=kind_phys),
intent(inout),
optional :: &
385 & dusfc_ms(:),dvsfc_ms(:), &
386 & dusfc_bl(:),dvsfc_bl(:), &
387 & dusfc_ss(:),dvsfc_ss(:), &
388 & dusfc_fd(:),dvsfc_fd(:)
389 real(kind=kind_phys),
intent(inout),
optional :: &
390 & dtaux2d_ms(:,:),dtauy2d_ms(:,:), &
391 & dtaux2d_bl(:,:),dtauy2d_bl(:,:), &
392 & dtaux2d_ss(:,:),dtauy2d_ss(:,:), &
393 & dtaux2d_fd(:,:),dtauy2d_fd(:,:)
396 real(kind=kind_phys),
dimension(im,km) :: dtaux2d, dtauy2d
402 logical,
intent(in) :: &
408 logical,
intent(in) :: ldiag_ugwp
412 logical,
intent(in) :: ugwp_seq_update
416 real(kind=kind_phys),
dimension(im,km) :: uwnd1, vwnd1
417 real(kind=kind_phys) :: tmp1, tmp2
420 logical,
parameter :: &
421 gwd_opt_ms = .true., &
422 gwd_opt_bl = .true., &
423 gsd_diss_ht_opt = .false.
427 real(kind=kind_phys),
parameter :: dxmin_ss = 1000., &
430 real(kind=kind_phys),
parameter :: dxmin_ms = 3000., &
432 real(kind=kind_phys),
dimension(im) :: ss_taper, ls_taper
435 real(kind=kind_phys),
parameter :: varmax_ss = 50., &
438 real(kind=kind_phys) :: var_temp, var_temp2
441 real(kind=kind_phys),
dimension(im,km) :: utendform,vtendform
442 real(kind=kind_phys) :: a1,a2,wsp
443 real(kind=kind_phys) :: h_efold
444 real(kind=kind_phys),
parameter :: coeff_fd = 6.325e-3
447 real(kind=kind_phys),
parameter :: ric = 0.25
448 real(kind=kind_phys),
parameter :: dw2min = 1.
449 real(kind=kind_phys),
parameter :: rimin = -100.
450 real(kind=kind_phys),
parameter :: bnv2min = 1.0e-5
451 real(kind=kind_phys),
parameter :: efmin = 0.0
452 real(kind=kind_phys),
parameter :: efmax = 10.0
453 real(kind=kind_phys),
parameter :: xl = 4.0e4
454 real(kind=kind_phys),
parameter :: critac = 1.0e-5
455 real(kind=kind_phys),
parameter :: gmax = 1.
456 real(kind=kind_phys),
parameter :: veleps = 1.0
457 real(kind=kind_phys),
parameter :: factop = 0.5
458 real(kind=kind_phys),
parameter :: frc = 1.0
459 real(kind=kind_phys),
parameter :: ce = 0.8
460 real(kind=kind_phys),
parameter :: cg = 0.5
461 real(kind=kind_phys),
parameter :: sgmalolev = 0.5
462 real(kind=kind_phys),
parameter :: plolevmeso = 70.0
463 real(kind=kind_phys),
parameter :: facmeso = 0.5
464 integer,
parameter :: kpblmin = 2
469 integer :: i,j,k,lcap,lcapp1,nwd,idir, &
472 real(kind=kind_phys) :: rcs,csg,fdir,cleff,cleff_ss,cs, &
473 rcsks,wdir,ti,rdz,tem2,dw2,shr2, &
474 bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, &
475 rim,temc,tem1,efact,temv,dtaux,dtauy, &
476 dtauxb,dtauyb,eng0,eng1,ksmax,dtfac_meso
478 logical :: ldrag(im),icrilv(im), &
482 real(kind=kind_phys) :: onebgrcs
484 real(kind=kind_phys) :: taub(im),taup(im,km+1), &
491 roll(im),dtfac(im), &
492 brvf(im),xlinv(im), &
493 delks(im),delks1(im), &
494 bnv2(im,km),usqj(im,km), &
495 taud_ms(im,km),taud_bl(im,km), &
497 vtk(im,km),vtj(im,km), &
498 zlowtop(im),velco(im,km-1), &
499 coefm(im),coefm_ss(im)
501 integer :: kbl(im),klowtop(im)
502 integer,
parameter :: mdir=8
505 integer,
parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/)
509 real(kind=kind_phys),
parameter :: frmax = 10.
510 real(kind=kind_phys),
parameter :: olmin = 1.0e-5
511 real(kind=kind_phys),
parameter :: odmin = 0.1
512 real(kind=kind_phys),
parameter :: odmax = 10.
515 real(kind=kind_phys) :: cd
516 real(kind=kind_phys) :: zblk,tautem
517 real(kind=kind_phys) :: pe,ke
518 real(kind=kind_phys) :: delx,dely
519 real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4)
520 real(kind=kind_phys) :: dxy(im),dxyp(im)
521 real(kind=kind_phys) :: ol4p(4),olp(im),od(im)
522 real(kind=kind_phys) :: taufb(im,km+1)
524 character(len=*),
intent(out) :: errmsg
525 integer,
intent(out) :: errflg
527 integer :: udtend, vdtend, Tdtend
543 udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd)
544 vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd)
545 tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd)
567 if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
580 cleff = 0.5e-5 / sqrt(float(imx)/192.0)
586 if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
599 fdir = mdir / (2.0*pi)
600 onebgrcs = 1._kind_phys/g*rcs
603 if (slmsk(i)==1. .or. slmsk(i)==2.)
then
612 if ( dx(i) .ge. dxmax_ms )
then
615 if ( dx(i) .le. dxmin_ms)
then
618 ls_taper(i) = 0.5 * ( sin(pi*(dx(i)-0.5*(dxmax_ms+dxmin_ms))/ &
619 (dxmax_ms-dxmin_ms)) + 1. )
628if ( spp_gwd==1 )
then
630 var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1)
631 varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1)
632 varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1)
636 var_stoch(i) = var(i)
637 varss_stoch(i) = varss(i)
638 varmax_fd_stoch(i) = varmax_fd
649 dxy4(i,3) = sqrt(delx*delx + dely*dely)
650 dxy4(i,4) = dxy4(i,3)
651 dxy4p(i,1) = dxy4(i,2)
652 dxy4p(i,2) = dxy4(i,1)
653 dxy4p(i,3) = dxy4(i,4)
654 dxy4p(i,4) = dxy4(i,3)
705 taufb(1:im,1:km+1) = 0.0
711 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
712 vtk(i,k) = vtj(i,k) / prslk(i,k)
713 ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k)
724 zl(i,k) = phil(i,k)*g_inv
731 zlowtop(i) = 2. * var_stoch(i)
740 if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i))
then
748 kbl(i) = max(kpbl(i), klowtop(i))
749 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
755 komax(:) = klowtop(:) - 1
758 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
759 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
766 if (k.lt.kbl(i))
then
767 rcsks = rcs * del(i,k) * delks(i)
768 rdelks = del(i,k) * delks(i)
769 ubar(i) = ubar(i) + rcsks * u1(i,k)
770 vbar(i) = vbar(i) + rcsks * v1(i,k)
771 roll(i) = roll(i) + rdelks * ro(i,k)
782 wdir = atan2(ubar(i),vbar(i)) + pi
783 idir = mod(nint(fdir*wdir),mdir) + 1
785 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
786 ol(i) = ol4(i,mod(nwd-1,4)+1)
788 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
789 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
799 olp(i) = ol4p(mod(nwd-1,4)+1)
803 od(i) = olp(i)/max(ol(i),olmin)
804 od(i) = min(od(i),odmax)
805 od(i) = max(od(i),odmin)
809 dxy(i) = dxy4(i,mod(nwd-1,4)+1)
810 dxyp(i) = dxy4p(i,mod(nwd-1,4)+1)
816IF ( (do_gsl_drag_ls_bl).and. &
817 (gwd_opt_ms.or.gwd_opt_bl) )
then
823 if ( ls_taper(i).GT.1.e-02 )
then
829 ti = 2.0 / (t1(i,k)+t1(i,k+1))
830 rdz = 1./(zl(i,k+1) - zl(i,k))
831 tem1 = u1(i,k) - u1(i,k+1)
832 tem2 = v1(i,k) - v1(i,k+1)
833 dw2 = rcl*(tem1*tem1 + tem2*tem2)
834 shr2 = max(dw2,dw2min) * rdz * rdz
835 bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
836 usqj(i,k) = max(bvf2/shr2,rimin)
837 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
838 bnv2(i,k) = max( bnv2(i,k), bnv2min )
843 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
844 rulow(i) = 1./ulow(i)
847 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
848 + (v1(i,k)+v1(i,k+1)) * vbar(i))
849 velco(i,k) = velco(i,k) * rulow(i)
850 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.))
then
857 ldrag(i) = velco(i,1).le.0.
861 do k = kpblmin,kpblmax
862 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
868 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
876 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
877 bnv2(i,1) = wtkbj * bnv2(i,1)
878 usqj(i,1) = wtkbj * usqj(i,1)
880 do k = kpblmin,kpblmax
881 if (k .lt. kbl(i))
then
882 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
883 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
884 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
888 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
889 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
890 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
896 if ( bnv2(i,1).gt.0.0 )
then
897 ldrag(i) = ldrag(i) .or. sqrt(bnv2(i,1))*rulow(i).lt.ksmax
902 do k = kpblmin,kpblmax
903 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
906 if (.not.ldrag(i))
then
907 bnv(i) = sqrt( bnv2(i,1) )
908 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
909 fr(i) = min(fr(i),frmax)
910 xn(i) = ubar(i) * rulow(i)
911 yn(i) = vbar(i) * rulow(i)
919 if (.not. ldrag(i))
then
920 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
921 efact = min( max(efact,efmin), efmax )
926 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
928 xlinv(i) = coefm(i) * cleff
929 tem = fr(i) * fr(i) * oc1(i)
930 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
931 if ( gwd_opt_ms )
then
932 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
933 * ulow(i) * gfobnv * efact
955IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ms) )
THEN
959 if ( ls_taper(i).GT.1.e-02 )
then
964 if (k .le. kbl(i)) taup(i,k) = taub(i)
974 if (k .ge. kbl(i))
then
975 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
976 .or. (velco(i,k) .le. 0.0)
977 brvf(i) = max(bnv2(i,k),bnv2min)
978 brvf(i) = sqrt(brvf(i))
981 if (k .ge. kbl(i) .and. (.not. ldrag(i)))
then
982 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 )
then
983 temv = 1.0 / velco(i,k)
984 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* &
986 hd = sqrt(taup(i,k) / tem1)
987 fro = brvf(i) * hd * temv
990 tem2 = sqrt(usqj(i,k))
991 tem = 1. + tem2 * fro
992 rim = usqj(i,k) * (1.-fro) / (tem * tem)
997 if (rim .le. ric)
then
998 if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin ))
then
999 temc = 2.0 + 1.0 / tem2
1000 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
1001 taup(i,kp1) = tem1 * hd * hd
1004 taup(i,kp1) = taup(i,k)
1011 do klcap = lcapp1,km
1012 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
1024IF ( do_gsl_drag_ls_bl .and. gwd_opt_bl )
THEN
1028 if ( ls_taper(i).GT.1.e-02 )
then
1030 if (.not.ldrag(i))
then
1036 do k = km, kpblmin, -1
1037 if(kblk.eq.0 .and. k.le.komax(i))
then
1038 pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))* &
1040 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
1046 kblk = min(kblk,kbl(i))
1047 zblk = zl(i,kblk)-zl(i,kts)
1048 rdxzb(i) = real(k,kind=kind_phys)
1056 cd = max(2.0-1.0/od(i),0.0)
1058 taufb(i,kts) = cdmb * 0.5 * roll(i) * coefm(i) / &
1059 max(dxmax_ms,dxy(i))**2 * cd * dxyp(i) * &
1060 olp(i) * zblk * ulow(i)**2
1061 tautem = taufb(i,kts)/float(kblk-kts)
1063 taufb(i,k) = taufb(i,k-1) - tautem
1079IF ( (do_gsl_drag_ls_bl) .and. &
1080 (gwd_opt_ms .OR. gwd_opt_bl) )
THEN
1084 if ( ls_taper(i) .GT. 1.e-02 )
then
1094 taup(i,km+1) = taup(i,km)
1096 taud_ms(i,k) = (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
1097 taud_bl(i,k) = (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
1107 do k = kts,kpblmax-1
1108 if (prsi(i,k).ge.sgmalolev*prsi(i,1))
then
1109 if ((taud_ms(i,k)+taud_bl(i,k)).ne.0.) &
1110 dtfac(i) = min(dtfac(i),abs(velco(i,k) &
1111 /(deltim*rcs*(taud_ms(i,k)+taud_bl(i,k)))))
1123 if (prsl(i,k).le.plolevmeso)
then
1124 if (taud_ms(i,k).ne.0.) &
1125 dtfac_meso = min(dtfac_meso,facmeso*abs(velco(i,k) &
1126 /(deltim*rcs*taud_ms(i,k))))
1129 taud_ms(i,k) = taud_ms(i,k)*dtfac(i)*dtfac_meso* &
1130 ls_taper(i) *(1.-rstoch(i))
1131 taud_bl(i,k) = taud_bl(i,k)*dtfac(i)* ls_taper(i) *(1.-rstoch(i))
1133 dtaux = taud_ms(i,k) * xn(i)
1134 dtauy = taud_ms(i,k) * yn(i)
1135 dtauxb = taud_bl(i,k) * xn(i)
1136 dtauyb = taud_bl(i,k) * yn(i)
1139 tmp1 = dtaux + dtauxb
1140 tmp2 = dtauy + dtauyb
1141 dudt(i,k) = tmp1 + dudt(i,k)
1142 dvdt(i,k) = tmp2 + dvdt(i,k)
1148 if ( ugwp_seq_update .and. (do_gsl_drag_ss.or.do_gsl_drag_tofd) )
then
1149 uwnd1(i,k) = uwnd1(i,k) + tmp1*deltim
1150 vwnd1(i,k) = vwnd1(i,k) + tmp2*deltim
1153 if ( gsd_diss_ht_opt )
then
1156 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
1158 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
1159 (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
1161 dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim
1162 if ( tdtend>0 )
then
1163 dtend(i,k,tdtend) = dtend(i,k,tdtend) + max((eng0-eng1),0.0)/cp
1167 dusfc(i) = dusfc(i) - onebgrcs * ( taud_ms(i,k)*xn(i)*del(i,k) + &
1168 taud_bl(i,k)*xn(i)*del(i,k) )
1169 dvsfc(i) = dvsfc(i) - onebgrcs * ( taud_ms(i,k)*yn(i)*del(i,k) + &
1170 taud_bl(i,k)*yn(i)*del(i,k) )
1172 dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ms(i,k) * &
1173 xn(i) + taud_bl(i,k) * xn(i)) * deltim
1176 dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ms(i,k) * &
1177 yn(i) + taud_bl(i,k) * yn(i)) * deltim
1182 if ( ldiag_ugwp )
then
1184 dtaux2d_ms(i,k) = taud_ms(i,k) * xn(i)
1185 dtauy2d_ms(i,k) = taud_ms(i,k) * yn(i)
1186 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
1187 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
1188 dusfc_ms(i) = dusfc_ms(i) - onebgrcs * dtaux2d_ms(i,k) * del(i,k)
1189 dvsfc_ms(i) = dvsfc_ms(i) - onebgrcs * dtauy2d_ms(i,k) * del(i,k)
1190 dusfc_bl(i) = dusfc_bl(i) - onebgrcs * dtaux2d_bl(i,k) * del(i,k)
1191 dvsfc_bl(i) = dvsfc_bl(i) - onebgrcs * dtauy2d_bl(i,k) * del(i,k)
1212IF ( do_gsl_drag_ss )
THEN
1216 if ( ss_taper(i).GT.1.e-02 )
then
1221 thx(i,k) = t1(i,k)/prslk(i,k)
1225 tvcon = (1.+fv*q1(i,k))
1226 thvx(i,k) = thx(i,k)*tvcon
1233 do k=kts+1,max(kpbl(i),kts+1)
1235 IF (zl(i,k)>300.)
then
1237 IF (k == kpbl(i))
then
1245 if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))
then
1246 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)
then
1253 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
1255 xnbv=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
1259 if (uwnd1(i,kpbl2).eq.0.)
then
1261 elseif (abs(xnbv/uwnd1(i,kpbl2)).gt.xlinv(i))
then
1268 var_temp = varss_stoch(i)
1270 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)
1271 tauwavex0=var_temp2*uwnd1(i,kvar)/(1.+var_temp2*deltim)
1272 tauwavex0=tauwavex0*ss_taper(i)
1279 if (vwnd1(i,kpbl2).eq.0.)
then
1281 elseif (abs(xnbv/vwnd1(i,kpbl2)).gt.xlinv(i))
then
1288 var_temp = varss_stoch(i)
1290 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)
1291 tauwavey0=var_temp2*vwnd1(i,kvar)/(1.+var_temp2*deltim)
1292 tauwavey0=tauwavey0*ss_taper(i)
1302 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
1303 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
1312 dudt(i,k) = dudt(i,k) + utendwave(i,k)
1313 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
1314 dusfc(i) = dusfc(i) - onebgrcs * utendwave(i,k) * del(i,k)
1315 dvsfc(i) = dvsfc(i) - onebgrcs * vtendwave(i,k) * del(i,k)
1318 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim
1321 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim
1323 if ( ldiag_ugwp )
then
1325 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
1326 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
1327 dtaux2d_ss(i,k) = utendwave(i,k)
1328 dtauy2d_ss(i,k) = vtendwave(i,k)
1342IF ( do_gsl_drag_tofd )
THEN
1346 if ( ss_taper(i).GT.1.e-02 )
then
1351 IF ((xland(i)-1.5) .le. 0.)
then
1353 var_temp = min(varss_stoch(i),varmax_fd_stoch(i)) + &
1354 max(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i)))
1355 a1=0.00026615161*var_temp**2
1363 wsp=sqrt(uwnd1(i,k)**2 + vwnd1(i,k)**2)
1367 var_temp = alpha_fd*coeff_fd*exp(-(zl(i,k)/h_efold)**1.5)*a2* &
1368 zl(i,k)**(-1.2)*ss_taper(i)
1371 utendform(i,k) = - var_temp*wsp*uwnd1(i,k)/(1. + var_temp*deltim*wsp)
1372 vtendform(i,k) = - var_temp*wsp*vwnd1(i,k)/(1. + var_temp*deltim*wsp)
1378 dudt(i,k) = dudt(i,k) + utendform(i,k)
1379 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
1380 dusfc(i) = dusfc(i) - onebgrcs * utendform(i,k) * del(i,k)
1381 dvsfc(i) = dvsfc(i) - onebgrcs * vtendform(i,k) * del(i,k)
1384 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim
1387 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim
1389 if ( ldiag_ugwp )
then
1391 dtaux2d_fd(i,k) = utendform(i,k)
1392 dtauy2d_fd(i,k) = vtendform(i,k)
1393 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
1394 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
1406if ( ldiag_ugwp )
then
1408 dusfc_ss(:) = -onebgrcs * dusfc_ss(:)
1409 dvsfc_ss(:) = -onebgrcs * dvsfc_ss(:)
1410 dusfc_fd(:) = -onebgrcs * dusfc_fd(:)
1411 dvsfc_fd(:) = -onebgrcs * dvsfc_fd(:)
1418 subroutine drag_suite_psl( &
1419 & IM,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, &
1420 & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, &
1421 & var,oc1,oa4,ol4, &
1422 & varss,oc1ss,oa4ss,ol4ss, &
1423 & THETA,SIGMA,GAMMA,ELVMAX, &
1424 & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, &
1425 & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, &
1427 & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
1428 & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
1429 & slmsk,br1,hpbl,vtype, &
1430 & g, cp, rd, rv, fv, pi, imx, cdmbgwd, alpha_fd, &
1431 & me, master, lprnt, ipr, rdxzb, dx, gwd_opt, &
1432 & do_gsl_drag_ls_bl, do_gsl_drag_ss, do_gsl_drag_tofd, &
1433 & psl_gwd_dx_factor, &
1434 & dtend, dtidx, index_of_process_orographic_gwd, &
1435 & index_of_temperature, index_of_x_wind, &
1436 & index_of_y_wind, ldiag3d, ldiag_ugwp, ugwp_seq_update, &
1437 & spp_wts_gwd, spp_gwd, errmsg, errflg)
1532 USE machine ,
ONLY : kind_phys
1536 integer,
intent(in) :: im, km, imx, kdt, ipr, me, master
1537 integer,
intent(in) :: gwd_opt
1538 logical,
intent(in) :: lprnt
1539 integer,
intent(in) :: KPBL(:)
1540 real(kind=kind_phys),
intent(in) :: deltim, g, cp, rd, rv, &
1541 & cdmbgwd(:), alpha_fd
1542 real(kind=kind_phys),
intent(inout),
optional :: dtend(:,:,:)
1543 logical,
intent(in) :: ldiag3d
1544 integer,
intent(in) :: dtidx(:,:), index_of_temperature, &
1545 & index_of_process_orographic_gwd, index_of_x_wind, index_of_y_wind
1548 integer,
parameter :: ims=1, kms=1, its=1, kts=1
1549 real(kind=kind_phys),
intent(in) :: fv, pi
1550 real(kind=kind_phys) :: rcl, cdmb
1551 real(kind=kind_phys) :: g_inv, g_cp, rd_inv
1553 real(kind=kind_phys),
intent(inout) :: &
1554 & dudt(:,:),dvdt(:,:), &
1556 real(kind=kind_phys),
intent(out) :: rdxzb(:)
1557 real(kind=kind_phys),
intent(in) :: &
1558 & u1(:,:),v1(:,:), &
1559 & t1(:,:),q1(:,:), &
1560 & phii(:,:),prsl(:,:), &
1561 & prslk(:,:),phil(:,:)
1562 real(kind=kind_phys),
intent(in) :: prsi(:,:), &
1564 real(kind=kind_phys),
intent(in) :: var(:),oc1(:), &
1565 & oa4(:,:),ol4(:,:), &
1567 real(kind=kind_phys),
intent(in),
optional :: varss(:),oc1ss(:), &
1568 & oa4ss(:,:),ol4ss(:,:)
1569 real(kind=kind_phys),
intent(in) :: theta(:),sigma(:), &
1570 & gamma(:),elvmax(:)
1573 real(kind=kind_phys),
dimension(im,km) :: utendwave,vtendwave,thx,thvx
1574 integer,
intent(in) :: vtype(:)
1575 real(kind=kind_phys),
intent(in) :: br1(:), &
1578 real(kind=kind_phys),
dimension(im) :: govrth,xland
1580 real(kind=kind_phys) :: tauwavex0,tauwavey0, &
1581 & xnbv,density,tvcon,hpbl2
1582 integer :: kpbl2,kvar
1584 real(kind=kind_phys),
dimension(im,km) :: zl
1587 real(kind=kind_phys),
dimension(im) :: var_stoch, varss_stoch, &
1588 varmax_ss_stoch, varmax_fd_stoch
1589 real(kind=kind_phys),
intent(in),
optional :: spp_wts_gwd(:,:)
1590 integer,
intent(in) :: spp_gwd
1592 real(kind=kind_phys),
dimension(im) :: rstoch
1595 real(kind=kind_phys),
intent(out) :: &
1596 & dusfc(:), dvsfc(:)
1598 real(kind=kind_phys),
intent(out),
optional :: &
1599 & dusfc_ls(:),dvsfc_ls(:), &
1600 & dusfc_bl(:),dvsfc_bl(:), &
1601 & dusfc_ss(:),dvsfc_ss(:), &
1602 & dusfc_fd(:),dvsfc_fd(:)
1603 real(kind=kind_phys),
intent(out),
optional :: &
1604 & dtaux2d_ls(:,:),dtauy2d_ls(:,:), &
1605 & dtaux2d_bl(:,:),dtauy2d_bl(:,:), &
1606 & dtaux2d_ss(:,:),dtauy2d_ss(:,:), &
1607 & dtaux2d_fd(:,:),dtauy2d_fd(:,:)
1610 real(kind=kind_phys),
dimension(im,km) :: dtaux2d, dtauy2d
1616 logical,
intent(in) :: &
1617 do_gsl_drag_ls_bl, &
1621 logical,
intent(in) :: ldiag_ugwp
1625 logical,
intent(in) :: ugwp_seq_update
1628 integer,
parameter :: &
1635 real(kind=kind_phys),
parameter :: dxmin_ss = 1000., &
1638 real(kind=kind_phys),
parameter :: dxmin_ls = 3000., &
1640 real(kind=kind_phys),
dimension(im) :: ss_taper, ls_taper
1643 real(kind=kind_phys),
parameter :: varmax_ss = 50., &
1647 real(kind=kind_phys) :: var_temp, var_temp2
1650 real(kind=kind_phys),
dimension(im,km) :: utendform,vtendform
1651 real(kind=kind_phys) :: a1,a2,wsp
1652 real(kind=kind_phys) :: h_efold
1653 real(kind=kind_phys),
parameter :: coeff_fd = 6.325e-3
1657 real(kind=kind_phys),
intent(in) :: psl_gwd_dx_factor
1660 real(kind=kind_phys),
parameter :: ric = 0.25
1661 real(kind=kind_phys),
parameter :: dw2min = 1.
1662 real(kind=kind_phys),
parameter :: rimin = -100.
1663 real(kind=kind_phys),
parameter :: bnv2min = 1.0e-5
1664 real(kind=kind_phys),
parameter :: efmin = 0.0
1665 real(kind=kind_phys),
parameter :: efmax = 10.0
1666 real(kind=kind_phys),
parameter :: xl = 4.0e4
1667 real(kind=kind_phys),
parameter :: critac = 1.0e-5
1668 real(kind=kind_phys),
parameter :: gmax = 1.
1669 real(kind=kind_phys),
parameter :: veleps = 1.0
1670 real(kind=kind_phys),
parameter :: factop = 0.5
1671 real(kind=kind_phys),
parameter :: frc = 1.0
1672 real(kind=kind_phys),
parameter :: ce = 0.8
1673 real(kind=kind_phys),
parameter :: cg = 1.0
1675 real(kind=kind_phys),
parameter :: var_min = 10.0
1676 real(kind=kind_phys),
parameter :: hmt_min = 50.
1677 real(kind=kind_phys),
parameter :: oc_min = 1.0
1678 real(kind=kind_phys),
parameter :: oc_max = 10.0
1680 real(kind=kind_phys),
parameter :: pcutoff = 7.5e2
1682 real(kind=kind_phys),
parameter :: pcutoff_den = 0.01
1684 integer,
parameter :: kpblmin = 2
1689 integer :: i,j,k,lcap,lcapp1,nwd,idir, &
1692 real(kind=kind_phys) :: rcs,csg,fdir,cs, &
1693 rcsks,wdir,ti,rdz,tem2,dw2,shr2, &
1694 bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, &
1695 rim,temc,tem1,efact,temv,dtaux,dtauy, &
1696 dtauxb,dtauyb,eng0,eng1
1697 real(kind=kind_phys) :: denfac
1699 logical :: ldrag(im),icrilv(im), &
1702 real(kind=kind_phys) :: invgrcs
1704 real(kind=kind_phys) :: taub(im),taup(im,km+1), &
1706 ubar(im),vbar(im), &
1708 rulow(im),bnv(im), &
1709 oa(im),ol(im),oc(im), &
1710 oass(im),olss(im), &
1711 roll(im),dtfac(im,km), &
1712 brvf(im),xlinv(im), &
1713 delks(im),delks1(im), &
1714 bnv2(im,km),usqj(im,km), &
1715 taud_ls(im,km),taud_bl(im,km), &
1717 vtk(im,km),vtj(im,km), &
1718 zlowtop(im),velco(im,km-1), &
1719 coefm(im),coefm_ss(im)
1720 real(kind=kind_phys) :: cleff(im),cleff_ss(im)
1722 integer :: kbl(im),klowtop(im)
1723 integer,
parameter :: mdir=8
1726 integer,
parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/)
1730 real(kind=kind_phys),
parameter :: frmax = 10.
1731 real(kind=kind_phys),
parameter :: olmin = 1.e-5
1732 real(kind=kind_phys),
parameter :: odmin = 0.1
1733 real(kind=kind_phys),
parameter :: odmax = 10.
1734 real(kind=kind_phys),
parameter :: cdmin = 0.0
1735 integer :: komax(im),kbmax(im),kblk(im)
1736 real(kind=kind_phys) :: hmax(im)
1737 real(kind=kind_phys) :: cd
1738 real(kind=kind_phys) :: zblk,tautem
1739 real(kind=kind_phys) :: pe,ke
1740 real(kind=kind_phys) :: delx,dely
1741 real(kind=kind_phys) :: dxy4(im,4),dxy4p(im,4)
1742 real(kind=kind_phys) :: dxy(im),dxyp(im)
1743 real(kind=kind_phys) :: ol4p(4),olp(im),od(im)
1744 real(kind=kind_phys) :: taufb(im,km+1)
1746 character(len=*),
intent(out) :: errmsg
1747 integer,
intent(out) :: errflg
1749 integer :: udtend, vdtend, Tdtend
1765 udtend = dtidx(index_of_x_wind,index_of_process_orographic_gwd)
1766 vdtend = dtidx(index_of_y_wind,index_of_process_orographic_gwd)
1767 tdtend = dtidx(index_of_temperature,index_of_process_orographic_gwd)
1778 fdir = mdir / (2.0*pi)
1779 invgrcs = 1._kind_phys/g*rcs
1784 if (slmsk(i)==1. .or. slmsk(i)==2.)
then
1794 if ( dx(i) .ge. dxmax_ls )
then
1797 if ( dx(i) .le. dxmin_ls)
then
1800 ls_taper(i) = 0.5 * ( sin(pi*(dx(i)-0.5*(dxmax_ls+dxmin_ls))/ &
1801 (dxmax_ls-dxmin_ls)) + 1. )
1810 if ( spp_gwd==1 )
then
1812 var_stoch(i) = var(i) + var(i)*0.75*spp_wts_gwd(i,1)
1813 varss_stoch(i) = varss(i) + varss(i)*0.75*spp_wts_gwd(i,1)
1814 varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.75*spp_wts_gwd(i,1)
1815 varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.75*spp_wts_gwd(i,1)
1819 var_stoch(i) = var(i)
1820 varss_stoch(i) = varss(i)
1821 varmax_ss_stoch(i) = varmax_ss
1822 varmax_fd_stoch(i) = varmax_fd
1833 dxy4(i,3) = sqrt(delx*delx + dely*dely)
1834 dxy4(i,4) = dxy4(i,3)
1835 dxy4p(i,1) = dxy4(i,2)
1836 dxy4p(i,2) = dxy4(i,1)
1837 dxy4p(i,3) = dxy4(i,4)
1838 dxy4p(i,4) = dxy4(i,3)
1839 cleff(i) = psl_gwd_dx_factor*(delx+dely)*0.5
1840 cleff_ss(i) = 0.1 * max(dxmax_ss,dxy4(i,3))
1886 if ( ldiag_ugwp )
then
1899 dtaux2d_ls(i,k)= 0.0
1900 dtauy2d_ls(i,k)= 0.0
1901 dtaux2d_bl(i,k)= 0.0
1902 dtauy2d_bl(i,k)= 0.0
1903 dtaux2d_ss(i,k)= 0.0
1904 dtauy2d_ss(i,k)= 0.0
1905 dtaux2d_fd(i,k)= 0.0
1906 dtauy2d_fd(i,k)= 0.0
1920 taufb(1:im,1:km+1) = 0.0
1929 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
1930 vtk(i,k) = vtj(i,k) / prslk(i,k)
1931 ro(i,k) = rd_inv * prsl(i,k) / vtj(i,k)
1942 zl(i,k) = phil(i,k)*g_inv
1949 if(vtype(i)==15)
then
1950 zlowtop(i) = 1.0 * var_stoch(i)
1952 zlowtop(i) = 2.0 * var_stoch(i)
1962 if(flag(i).and.zl(i,k).ge.zlowtop(i))
then
1979 if(flag(i).and.zl(i,k).ge.elvmax(i))
then
1994 if(flag(i).and.zl(i,k).ge.elvmax(i)+zlowtop(i))
then
2004 hmax(i) = max(elvmax(i),zlowtop(i))
2009 kbl(i) = max(komax(i), klowtop(i))
2010 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
2016 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
2017 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
2021 if (k.lt.kbl(i))
then
2022 rcsks = rcs * del(i,k) * delks(i)
2023 rdelks = del(i,k) * delks(i)
2024 ubar(i) = ubar(i) + rcsks * u1(i,k)
2025 vbar(i) = vbar(i) + rcsks * v1(i,k)
2026 roll(i) = roll(i) + rdelks * ro(i,k)
2037 wdir = atan2(ubar(i),vbar(i)) + pi
2038 idir = mod(nint(fdir*wdir),mdir) + 1
2040 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
2041 ol(i) = max(ol4(i,mod(nwd-1,4)+1),olmin)
2042 oc(i) = min(max(oc1(i),oc_min),oc_max)
2047 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
2048 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
2058 olp(i) = max(ol4p(mod(nwd-1,4)+1),olmin)
2062 od(i) = olp(i)/ol(i)
2063 od(i) = min(od(i),odmax)
2064 od(i) = max(od(i),odmin)
2068 dxy(i) = dxy4(i,mod(nwd-1,4)+1)
2069 dxyp(i) = dxy4p(i,mod(nwd-1,4)+1)
2074IF ( (do_gsl_drag_ls_bl).and. &
2075 ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)) )
then
2080 if ( ls_taper(i).GT.1.e-02 )
then
2086 ti = 2.0 / (t1(i,k)+t1(i,k+1))
2087 rdz = 1./(zl(i,k+1) - zl(i,k))
2088 tem1 = u1(i,k) - u1(i,k+1)
2089 tem2 = v1(i,k) - v1(i,k+1)
2090 dw2 = rcl*(tem1*tem1 + tem2*tem2)
2091 shr2 = max(dw2,dw2min) * rdz * rdz
2092 bvf2 = g*(g_cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
2093 usqj(i,k) = max(bvf2/shr2,rimin)
2094 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
2095 bnv2(i,k) = max( bnv2(i,k), bnv2min )
2100 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
2101 rulow(i) = 1./ulow(i)
2104 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
2105 + (v1(i,k)+v1(i,k+1)) * vbar(i))
2106 velco(i,k) = velco(i,k) * rulow(i)
2107 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.))
then
2114 ldrag(i) = hmax(i).le.hmt_min
2118 ldrag(i) = ldrag(i).or. velco(i,1).le.0.
2122 do k = kpblmin,kpblmax
2123 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
2129 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
2137 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
2138 bnv2(i,1) = wtkbj * bnv2(i,1)
2139 usqj(i,1) = wtkbj * usqj(i,1)
2141 do k = kpblmin,kpblmax
2142 if (k .lt. kbl(i))
then
2143 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
2144 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
2145 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
2149 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
2150 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
2151 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
2152 ldrag(i) = ldrag(i) .or. xland(i) .gt. 1.5
2156 do k = kpblmin,kpblmax
2157 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
2160 if (.not.ldrag(i))
then
2161 bnv(i) = sqrt( bnv2(i,1) )
2162 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
2163 fr(i) = min(fr(i),frmax)
2164 xn(i) = ubar(i) * rulow(i)
2165 yn(i) = vbar(i) * rulow(i)
2173 if (.not. ldrag(i))
then
2174 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
2175 efact = min( max(efact,efmin), efmax )
2176 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
2177 xlinv(i) = coefm(i) / cleff(i)
2178 tem = fr(i) * fr(i) * oc(i)
2179 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
2180 if ( gwd_opt_ls .NE. 0 )
then
2181 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
2182 * ulow(i) * gfobnv * efact
2208IF ( do_gsl_drag_ss )
THEN
2212 if ( ss_taper(i).GT.1.e-02 )
then
2217 thx(i,k) = t1(i,k)/prslk(i,k)
2221 tvcon = (1.+fv*q1(i,k))
2222 thvx(i,k) = thx(i,k)*tvcon
2229 do k=kts+1,max(kpbl(i),kts+1)
2231 IF (zl(i,k)>300.)
then
2233 IF (k == kpbl(i))
then
2241 if((xland(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))
then
2242 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)
then
2243 coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
2244 xlinv(i) = coefm_ss(i) / cleff_ss(i)
2246 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
2248 xnbv=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
2251 if(abs(xnbv/u1(i,kpbl2)).gt.xlinv(i))
then
2256 var_temp = varss_stoch(i)
2258 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)
2259 tauwavex0=var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
2260 tauwavex0=tauwavex0*ss_taper(i)
2266 if(abs(xnbv/v1(i,kpbl2)).gt.xlinv(i))
then
2271 var_temp = varss_stoch(i)
2273 var_temp2 = 0.5*xnbv*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)
2274 tauwavey0=var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
2275 tauwavey0=tauwavey0*ss_taper(i)
2285 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2286 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
2295 dudt(i,k) = dudt(i,k) + utendwave(i,k)
2296 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
2297 dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k)
2298 dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k)
2301 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendwave(i,kts:km)*deltim
2304 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendwave(i,kts:km)*deltim
2306 if ( ldiag_ugwp )
then
2308 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
2309 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
2310 dtaux2d_ss(i,k) = utendwave(i,k)
2311 dtauy2d_ss(i,k) = vtendwave(i,k)
2324IF ( do_gsl_drag_tofd )
THEN
2328 if ( ss_taper(i).GT.1.e-02 )
then
2333 IF ((xland(i)-1.5) .le. 0.)
then
2336 var_temp = varss_stoch(i)
2338 a1=0.00026615161*var_temp**2
2346 wsp=sqrt(u1(i,k)**2 + v1(i,k)**2)
2350 var_temp = alpha_fd*coeff_fd*exp(-(zl(i,k)/h_efold)**1.5)*a2* &
2351 zl(i,k)**(-1.2)*ss_taper(i)
2354 utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp)
2355 vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp)
2361 dudt(i,k) = dudt(i,k) + utendform(i,k)
2362 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
2363 dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k)
2364 dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k)
2367 dtend(i,kts:km,udtend) = dtend(i,kts:km,udtend) + utendform(i,kts:km)*deltim
2370 dtend(i,kts:km,vdtend) = dtend(i,kts:km,vdtend) + vtendform(i,kts:km)*deltim
2372 if ( ldiag_ugwp )
then
2374 dtaux2d_fd(i,k) = utendform(i,k)
2375 dtauy2d_fd(i,k) = vtendform(i,k)
2376 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
2377 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
2388IF ( (do_gsl_drag_ls_bl).and.(gwd_opt_ls .EQ. 1) )
THEN
2392 if ( ls_taper(i).GT.1.e-02 )
then
2397 if (k .le. kbl(i)) taup(i,k) = taub(i)
2400 do k = kpblmin, km-1
2407 if (k .ge. kbl(i))
then
2408 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
2409 .or. (velco(i,k) .le. 0.0)
2410 brvf(i) = max(bnv2(i,k),bnv2min)
2411 brvf(i) = sqrt(brvf(i))
2414 if (k .ge. kbl(i) .and. (.not. ldrag(i)))
then
2415 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 )
then
2416 temv = 1.0 / velco(i,k)
2417 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)* &
2419 hd = sqrt(taup(i,k) / tem1)
2420 fro = brvf(i) * hd * temv
2423 tem2 = sqrt(usqj(i,k))
2424 tem = 1. + tem2 * fro
2425 rim = usqj(i,k) * (1.-fro) / (tem * tem)
2430 if (rim .le. ric)
then
2431 temc = 2.0 + 1.0 / tem2
2432 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
2433 taup(i,kp1) = tem1 * hd * hd
2435 taup(i,kp1) = taup(i,k)
2442 do klcap = lcapp1,km
2443 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
2455IF ( (do_gsl_drag_ls_bl) .and. (gwd_opt_bl .EQ. 1) )
THEN
2462 if ( ls_taper(i).GT.1.e-02 )
then
2464 if (.not.ldrag(i))
then
2470 do k = km, kpblmin, -1
2471 if(flag(i).and. k.le.kbmax(i))
then
2472 pe = pe + bnv2(i,k)*(zl(i,kbmax(i))-zl(i,k))* &
2473 del(i,k)*g_inv/ro(i,k)
2474 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
2478 if(pe.ge.ke.and.zl(i,k).le.hmax(i))
then
2481 rdxzb(i) = real(k,kind=kind_phys)
2486 if(.not.flag(i))
then
2490 cd = max(2.0-1.0/od(i),cdmin)
2491 taufb(i,kts) = 0.5 * roll(i) * coefm(i) / &
2492 max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) * &
2493 olp(i) * zblk * ulow(i)**2
2494 tautem = taufb(i,kts)/float(kblk(i)-kts)
2495 do k = kts+1, kpblmax
2496 if (k .le. kblk(i)) taufb(i,k) = taufb(i,k-1) - tautem
2502 if (k .le. kblk(i)) taup(i,k) = taup(i,kblk(i))
2516IF ( (do_gsl_drag_ls_bl) .and. &
2517 (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) )
THEN
2521 if ( ls_taper(i) .GT. 1.e-02 )
then
2527 taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
2528 taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
2534 taud_ls(i,klcap) = taud_ls(i,klcap) * factop
2535 taud_bl(i,klcap) = taud_bl(i,klcap) * factop
2542 do k = kts,kpblmax-1
2543 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0.)
then
2544 dtfac(i,k) = min(dtfac(i,k),abs(velco(i,k) &
2545 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2552 if ((taud_ls(i,k)+taud_bl(i,k)).ne.0..and.prsl(i,k).le.pcutoff)
then
2553 denfac = min(ro(i,k)/pcutoff_den,1.)
2554 dtfac(i,k) = min(dtfac(i,k),denfac*abs(velco(i,k) &
2555 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
2560 taud_ls(i,k) = taud_ls(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2561 taud_bl(i,k) = taud_bl(i,k)*dtfac(i,k)* ls_taper(i) *(1.-rstoch(i))
2562 dtaux = taud_ls(i,k) * xn(i)
2563 dtauy = taud_ls(i,k) * yn(i)
2564 dtauxb = taud_bl(i,k) * xn(i)
2565 dtauyb = taud_bl(i,k) * yn(i)
2568 dudt(i,k) = dtaux + dtauxb + dudt(i,k)
2569 dvdt(i,k) = dtauy + dtauyb + dvdt(i,k)
2571 if ( gsd_diss_ht_opt .EQ. 1 )
then
2574 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
2576 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
2577 (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
2579 dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim
2580 if ( tdtend>0 )
then
2581 dtend(i,k,tdtend) = dtend(i,k,tdtend) + max((eng0-eng1),0.0)/cp
2585 dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + &
2586 taud_bl(i,k)*xn(i)*del(i,k)
2587 dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + &
2588 taud_bl(i,k)*yn(i)*del(i,k)
2590 dtend(i,k,udtend) = dtend(i,k,udtend) + (taud_ls(i,k) * &
2591 xn(i) + taud_bl(i,k) * xn(i)) * deltim
2594 dtend(i,k,vdtend) = dtend(i,k,vdtend) + (taud_ls(i,k) * &
2595 yn(i) + taud_bl(i,k) * yn(i)) * deltim
2601 dusfc(i) = -(invgrcs) * dusfc(i)
2602 dvsfc(i) = -(invgrcs) * dvsfc(i)
2604 if ( ldiag_ugwp )
then
2606 dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i)
2607 dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i)
2608 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
2609 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
2610 dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k)
2611 dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k)
2612 dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k)
2613 dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k)
2623if ( ldiag_ugwp )
then
2626 dusfc_ls(i) = -(invgrcs) * dusfc_ls(i)
2627 dvsfc_ls(i) = -(invgrcs) * dvsfc_ls(i)
2628 dusfc_bl(i) = -(invgrcs) * dusfc_bl(i)
2629 dvsfc_bl(i) = -(invgrcs) * dvsfc_bl(i)
2630 dusfc_ss(i) = -(invgrcs) * dusfc_ss(i)
2631 dvsfc_ss(i) = -(invgrcs) * dvsfc_ss(i)
2632 dusfc_fd(i) = -(invgrcs) * dusfc_fd(i)
2633 dvsfc_fd(i) = -(invgrcs) * dvsfc_fd(i)