!WRF:model_layer:physics ! ! ! ! module a_module_bl_gwdo contains ! ! Generated by TAPENADE (INRIA, Tropics team) ! Tapenade 3.10 (r5498) - 20 Jan 2015 09:48 ! ! Differentiation of gwdo in reverse (adjoint) mode: ! gradient of useful results: rublten p3di u3d dusfcg z dvsfcg ! dtauy3d rvblten t3d qv3d pi3d v3d dtaux3d ! with respect to varying inputs: rublten p3di u3d dusfcg z dvsfcg ! dtauy3d rvblten t3d qv3d pi3d v3d dtaux3d ! RW status of diff variables: rublten:in-out p3di:incr u3d:incr ! dusfcg:in-out z:incr dvsfcg:in-out dtauy3d:in-out ! rvblten:in-out t3d:incr qv3d:incr pi3d:incr v3d:incr ! dtaux3d:in-out !------------------------------------------------------------------------------- SUBROUTINE GWDO_B(u3d, u3db, v3d, v3db, t3d, t3db, qv3d, qv3db, p3d, & & p3di, p3dib, pi3d, pi3db, z, zb, rublten, rubltenb, rvblten, rvbltenb& & , dtaux3d, dtaux3db, dtauy3d, dtauy3db, dusfcg, dusfcgb, dvsfcg, & & dvsfcgb, var2d, oc12d, oa2d1, oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3& & , ol2d4, znu, znw, p_top, cp, g, rd, rv, ep1, pi, dt, dx, & & kpbl2d, itimestep, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, & & kms, kme, its, ite, jts, jte, kts, kte) IMPLICIT NONE ! !------------------------------------------------------------------------------- ! !-- u3d 3d u-velocity interpolated to theta points (m/s) !-- v3d 3d v-velocity interpolated to theta points (m/s) !-- t3d temperature (k) !-- qv3d 3d water vapor mixing ratio (kg/kg) !-- p3d 3d pressure (pa) !-- p3di 3d pressure (pa) at interface level !-- pi3d 3d exner function (dimensionless) !-- rublten u tendency due to pbl parameterization (m/s/s) !-- rvblten v tendency due to pbl parameterization (m/s/s) !-- znu eta values (sigma values) !-- cp heat capacity at constant pressure for dry air (j/kg/k) !-- g acceleration due to gravity (m/s^2) !-- rd gas constant for dry air (j/kg/k) !-- z height above sea level (m) !-- rv gas constant for water vapor (j/kg/k) !-- dt time step (s) !-- dx model grid interval (m) !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless) !-- ids start index for i in domain !-- ide end index for i in domain !-- jds start index for j in domain !-- jde end index for j in domain !-- kds start index for k in domain !-- kde end index for k in domain !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !-- its start index for i in tile !-- ite end index for i in tile !-- jts start index for j in tile !-- jte end index for j in tile !-- kts start index for k in tile !-- kte end index for k in tile ! !------------------------------------------------------------------------------- INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, & & jme, kms, kme, its, ite, jts, jte, kts, kte INTEGER, INTENT(IN) :: itimestep ! REAL, INTENT(IN) :: dt, dx, cp, g, rd, rv, ep1, pi ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: qv3d, p3d, & & pi3d, t3d, z REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: qv3db, pi3db, t3db, zb REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: p3di REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: p3dib ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rublten, & & rvblten REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rubltenb REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3d, & & dtauy3d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtaux3db ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: u3d, v3d REAL, DIMENSION(ims:ime, kms:kme, jms:jme) :: u3db, v3db ! INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl2d REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcg, dvsfcg REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dusfcgb ! REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: var2d, oc12d, oa2d1, & & oa2d2, oa2d3, oa2d4, ol2d1, ol2d2, ol2d3, ol2d4 ! REAL, DIMENSION(kms:kme), OPTIONAL, INTENT(IN) :: znu, znw ! REAL, OPTIONAL, INTENT(IN) :: p_top ! !local ! REAL, DIMENSION(its:ite, kts:kte) :: delprsi, pdh REAL, DIMENSION(its:ite, kts:kte) :: delprsib, pdhb REAL, DIMENSION(its:ite, kts:kte + 1) :: pdhi REAL, DIMENSION(its:ite, kts:kte+1) :: pdhib REAL, DIMENSION(its:ite, 4) :: oa4, ol4 INTEGER :: i, j, k, kdt, kpblmax INTEGER :: branch REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: rvbltenb REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: dvsfcgb REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: dtauy3db ! DO k=kts,kte IF (znu(k) .GT. 0.6) kpblmax = k + 1 END DO ! DO j=jts,jte DO k=kts,kte+1 DO i=its,ite IF (k .LE. kte) THEN CALL PUSHREAL8(pdh(i, k)) pdh(i, k) = p3d(i, k, j) CALL PUSHCONTROL1B(0) ELSE CALL PUSHCONTROL1B(1) END IF CALL PUSHREAL8(pdhi(i, k)) pdhi(i, k) = p3di(i, k, j) END DO END DO ! DO k=kts,kte DO i=its,ite CALL PUSHREAL8(delprsi(i, k)) delprsi(i, k) = pdhi(i, k) - pdhi(i, k+1) END DO END DO DO i=its,ite CALL PUSHREAL8(oa4(i, 1)) oa4(i, 1) = oa2d1(i, j) CALL PUSHREAL8(oa4(i, 2)) oa4(i, 2) = oa2d2(i, j) CALL PUSHREAL8(oa4(i, 3)) oa4(i, 3) = oa2d3(i, j) CALL PUSHREAL8(oa4(i, 4)) oa4(i, 4) = oa2d4(i, j) CALL PUSHREAL8(ol4(i, 1)) ol4(i, 1) = ol2d1(i, j) CALL PUSHREAL8(ol4(i, 2)) ol4(i, 2) = ol2d2(i, j) CALL PUSHREAL8(ol4(i, 3)) ol4(i, 3) = ol2d3(i, j) CALL PUSHREAL8(ol4(i, 4)) ol4(i, 4) = ol2d4(i, j) END DO END DO delprsib = 0.0 pdhib = 0.0 pdhb = 0.0 DO j=jte,jts,-1 CALL GWDO2D_B(rublten(ims, kms, j), rubltenb(ims, kms, j), rvblten(& & ims, kms, j), rvbltenb(ims, kms, j), dtaux3d(ims, kms, j), & & dtaux3db(ims, kms, j), dtauy3d(ims, kms, j), dtauy3db(ims, & & kms, j), u3d(ims, kms, j), u3db(ims, kms, j), v3d(ims, kms, & & j), v3db(ims, kms, j), t3d(ims, kms, j), t3db(ims, kms, j), & & qv3d(ims, kms, j), qv3db(ims, kms, j), delprsi(its, kts), & & delprsib(its, kts), pdhi(its, kts), pdhib(its, kts), pdh(its& & , kts), pdhb(its, kts), pi3d(ims, kms, j), pi3db(ims, kms, j& & ), z(ims, kms, j), zb(ims, kms, j), rcl=1.0, kpblmax=kpblmax& & , var=var2d(ims, j), oc1=oc12d(ims, j), oa4=oa4, ol4=ol4, & & dusfc=dusfcg(ims, j), dusfcb=dusfcgb(ims, j), dvsfc=dvsfcg(& & ims, j), dvsfcb=dvsfcgb(ims, j), g=g, cp=cp, rd=rd, rv=rv, & & fv=ep1, pi=pi, dxmeter=dx, deltim=dt, kpbl=kpbl2d(ims, j), & & kdt=itimestep, lat=j, ids=ids, ide=ide, jds=jds, jde=jde, & & kds=kds, kde=kde, ims=ims, ime=ime, jms=jms, jme=jme, kms=& & kms, kme=kme, its=its, ite=ite, jts=jts, jte=jte, kts=kts, & & kte=kte) DO i=ite,its,-1 CALL POPREAL8(ol4(i, 4)) CALL POPREAL8(ol4(i, 3)) CALL POPREAL8(ol4(i, 2)) CALL POPREAL8(ol4(i, 1)) CALL POPREAL8(oa4(i, 4)) CALL POPREAL8(oa4(i, 3)) CALL POPREAL8(oa4(i, 2)) CALL POPREAL8(oa4(i, 1)) END DO DO k=kte,kts,-1 DO i=ite,its,-1 CALL POPREAL8(delprsi(i, k)) pdhib(i, k) = pdhib(i, k) + delprsib(i, k) pdhib(i, k+1) = pdhib(i, k+1) - delprsib(i, k) delprsib(i, k) = 0.0 END DO END DO DO k=kte+1,kts,-1 DO i=ite,its,-1 CALL POPREAL8(pdhi(i, k)) p3dib(i, k, j) = p3dib(i, k, j) + pdhib(i, k) pdhib(i, k) = 0.0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL8(pdh(i, k)) pdhb(i, k) = 0.0 END IF END DO END DO END DO END SUBROUTINE GWDO_B ! Differentiation of gwdo2d in reverse (adjoint) mode: ! gradient of useful results: v1 dvsfc dvdt dtauy2d prsi ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk ! with respect to varying inputs: v1 dvsfc dvdt dtauy2d prsi ! prsl dusfc del t1 q1 dudt dtaux2d u1 zl prslk !------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------- SUBROUTINE GWDO2D_B(dudt, dudtb, dvdt, dvdtb, dtaux2d, dtaux2db, dtauy2d& & , dtauy2db, u1, u1b, v1, v1b, t1, t1b, q1, q1b, del, delb, prsi, prsib& & , prsl, prslb, prslk, prslkb, zl, zlb, rcl, kpblmax, var, oc1, oa4, & & ol4, dusfc, dusfcb, dvsfc, dvsfcb, g, cp, rd, rv, fv, pi, dxmeter, & & deltim, kpbl, kdt, lat, ids, ide, jds, jde, kds, kde, ims, ime, jms, & & jme, kms, kme, its, ite, jts, jte, kts, kte) IMPLICIT NONE !------------------------------------------------------------------------------- INTEGER :: kdt, lat, latd, lond, kpblmax, ids, ide, jds, jde, kds, kde& & , ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte ! REAL :: g, rd, rv, fv, cp, pi, dxmeter, deltim, rcl REAL :: dudt(ims:ime, kms:kme), dvdt(ims:ime, kms:kme), dtaux2d(ims:& & ime, kms:kme), dtauy2d(ims:ime, kms:kme), u1(ims:ime, kms:kme), v1(ims& & :ime, kms:kme), t1(ims:ime, kms:kme), q1(ims:ime, kms:kme), zl(ims:ime& & , kms:kme), prsl(its:ite, kts:kte), prslk(ims:ime, kms:kme) REAL :: dudtb(ims:ime, kms:kme), dvdtb(ims:ime, kms:kme), dtaux2db(ims& & :ime, kms:kme), dtauy2db(ims:ime, kms:kme), u1b(ims:ime, kms:kme), v1b& & (ims:ime, kms:kme), t1b(ims:ime, kms:kme), q1b(ims:ime, kms:kme), zlb(& & ims:ime, kms:kme), prslb(its:ite, kts:kte), prslkb(ims:ime, kms:kme) REAL :: prsi(its:ite, kts:kte+1), del(its:ite, kts:kte) REAL :: prsib(its:ite, kts:kte+1), delb(its:ite, kts:kte) REAL :: oa4(its:ite, 4), ol4(its:ite, 4) ! INTEGER :: kpbl(ims:ime) REAL :: var(ims:ime), oc1(ims:ime), dusfc(ims:ime), dvsfc(ims:ime) REAL :: dusfcb(ims:ime), dvsfcb(ims:ime) ! ! critical richardson number for wave breaking : ! larger drag with larger value ! REAL, PARAMETER :: ric=0.25 ! REAL, PARAMETER :: dw2min=1. REAL, PARAMETER :: rimin=-100. REAL, PARAMETER :: bnv2min=1.0e-5 REAL, PARAMETER :: efmin=0.0 REAL, PARAMETER :: efmax=10.0 REAL, PARAMETER :: xl=4.0e4 REAL, PARAMETER :: critac=1.0e-5 REAL, PARAMETER :: gmax=1. REAL, PARAMETER :: veleps=1.0 REAL, PARAMETER :: factop=0.5 REAL, PARAMETER :: frc=1.0 REAL, PARAMETER :: ce=0.8 REAL, PARAMETER :: cg=0.5 INTEGER, PARAMETER :: kpblmin=2 ! ! local variables ! INTEGER :: i, k, lcap, lcapp1, nwd, idir, klcap, kp1, ikount, kk ! REAL :: rcs, rclcs, csg, fdir, cleff, cs, rcsks, wdir, ti, rdz, temp, & & tem2, dw2, shr2, bvf2, rdelks, wtkbj, tem, gfobnv, hd, fro, rim, temc& & , tem1, efact, temv, dtaux, dtauy REAL :: rcsksb, tib, rdzb, tem2b, dw2b, shr2b, bvf2b, rdelksb, wtkbjb& & , temb, gfobnvb, hdb, temcb, tem1b, efactb, dtauxb, dtauyb ! LOGICAL :: ldrag(its:ite), icrilv(its:ite), flag(its:ite), kloop1(its:& & ite) ! REAL :: taub(its:ite), taup(its:ite, kts:kte+1), xn(its:ite), yn(its:& & ite), ubar(its:ite), vbar(its:ite), fr(its:ite), ulow(its:ite), rulow(& & its:ite), bnv(its:ite), oa(its:ite), ol(its:ite), roll(its:ite), dtfac& & (its:ite), brvf(its:ite), xlinv(its:ite), delks(its:ite), delks1(its:& & ite), bnv2(its:ite, kts:kte), usqj(its:ite, kts:kte), taud(its:ite, & & kts:kte), ro(its:ite, kts:kte), vtk(its:ite, kts:kte), vtj(its:ite, & & kts:kte), zlowtop(its:ite), velco(its:ite, kts:kte-1), coefm(its:ite) REAL :: taubb(its:ite), taupb(its:ite, kts:kte+1), xnb(its:ite), ynb(& & its:ite), ubarb(its:ite), vbarb(its:ite), frb(its:ite), ulowb(its:ite)& & , rulowb(its:ite), bnvb(its:ite), rollb(its:ite), dtfacb(its:ite), & & brvfb(its:ite), delksb(its:ite), delks1b(its:ite), bnv2b(its:ite, kts:& & kte), usqjb(its:ite, kts:kte), taudb(its:ite, kts:kte), rob(its:ite, & & kts:kte), vtkb(its:ite, kts:kte), vtjb(its:ite, kts:kte), velcob(its:& & ite, kts:kte-1) ! INTEGER :: kbl(its:ite), klowtop(its:ite) ! LOGICAL :: iope INTEGER, PARAMETER :: mdir=8 INTEGER :: nwdir(mdir) ! ! variables for flow-blocking drag ! REAL, PARAMETER :: frmax=10. REAL, PARAMETER :: olmin=1.0e-5 REAL, PARAMETER :: odmin=0.1 REAL, PARAMETER :: odmax=10. REAL, PARAMETER :: erad=6371.315e+3 INTEGER :: komax(its:ite) INTEGER :: kblk REAL :: cd REAL :: zblk, tautem REAL :: zblkb, tautemb REAL :: pe, ke REAL :: delx, dely, dxy4(4), dxy4p(4) REAL :: dxy(its:ite), dxyp(its:ite) REAL :: ol4p(4), olp(its:ite), od(its:ite) REAL :: taufb(its:ite, kts:kte+1) REAL :: taufbb(its:ite, kts:kte+1) REAL :: tmp REAL :: tmp0 REAL :: tmp1 REAL :: tmp2 REAL :: tmp3 INTEGER :: branch INTEGER :: ad_to REAL :: temp3 REAL :: x3b REAL :: temp2 REAL :: y1b REAL :: temp1 REAL :: temp0 REAL :: tempb9 REAL :: tempb8 REAL :: tempb7 REAL :: max2b REAL :: tempb6 REAL :: tempb5 REAL :: tempb4 REAL :: tempb19 REAL :: tempb3 REAL :: tempb18 REAL :: tempb2 REAL :: tempb17 REAL :: tempb1 REAL :: tempb16 REAL :: tempb0 REAL :: tempb15 REAL :: tempb14 REAL :: tempb13 REAL :: tmpb3 REAL :: tmpb REAL :: tempb12 REAL :: tmpb2 REAL :: tmpb1 REAL :: tempb11 REAL :: tmpb0 REAL :: tempb10 REAL :: x4 REAL :: x3 REAL :: x2 INTEGER :: x1 REAL :: x2b REAL :: temp12 REAL :: temp11 REAL :: temp10 REAL :: tempb REAL :: x4b REAL :: max3 REAL :: max2 REAL :: max1 REAL :: temp9 REAL :: temp8 REAL :: tempb20 REAL :: temp7 REAL :: temp6 REAL :: temp5 REAL :: y1 REAL :: temp4 DATA nwdir /6, 7, 5, 8, 2, 3, 1, 4/ ! !---- constants ! rcs = SQRT(rcl) cs = 1./SQRT(rcl) csg = cs*g lcap = kte lcapp1 = lcap + 1 fdir = mdir/(2.0*pi) ! !--- calculate length of grid for flow-blocking drag ! delx = dxmeter dely = dxmeter dxy4(1) = delx dxy4(2) = dely dxy4(3) = SQRT(delx*delx + dely*dely) dxy4(4) = dxy4(3) dxy4p(1) = dxy4(2) dxy4p(2) = dxy4(1) dxy4p(3) = dxy4(4) dxy4p(4) = dxy4(3) ! ! !-----initialize arrays ! DO i=its,ite klowtop(i) = 0 kbl(i) = 0 END DO ! DO i=its,ite xn(i) = 0.0 yn(i) = 0.0 ubar(i) = 0.0 vbar(i) = 0.0 roll(i) = 0.0 taub(i) = 0.0 taup(i, 1) = 0.0 oa(i) = 0.0 ol(i) = 0.0 ulow(i) = 0.0 dtfac(i) = 1.0 ldrag(i) = .false. icrilv(i) = .false. END DO ! DO k=kts,kte DO i=its,ite usqj(i, k) = 0.0 bnv2(i, k) = 0.0 vtj(i, k) = 0.0 vtk(i, k) = 0.0 taup(i, k) = 0.0 taud(i, k) = 0.0 END DO END DO ! DO i=its,ite taup(i, kte+1) = 0.0 END DO ! ! initialize array for flow-blocking drag ! taufb(its:ite, kts:kte+1) = 0.0 ! DO k=kts,kte DO i=its,ite vtj(i, k) = t1(i, k)*(1.+fv*q1(i, k)) vtk(i, k) = vtj(i, k)/prslk(i, k) ! density kg/m**3 ro(i, k) = 1./rd*prsl(i, k)/vtj(i, k) END DO END DO ! ! determine reference level: maximum of 2*var and pbl heights ! DO i=its,ite zlowtop(i) = 2.*var(i) END DO ! DO i=its,ite kloop1(i) = .true. END DO ! DO k=kts+1,kte DO i=its,ite IF (kloop1(i) .AND. zl(i, k) - zl(i, 1) .GE. zlowtop(i)) THEN klowtop(i) = k + 1 kloop1(i) = .false. END IF END DO END DO ! DO i=its,ite IF (kpbl(i) .LT. klowtop(i)) THEN kbl(i) = klowtop(i) ELSE kbl(i) = kpbl(i) END IF IF (kbl(i) .GT. kpblmax) THEN x1 = kpblmax ELSE x1 = kbl(i) END IF IF (x1 .LT. kpblmin) THEN kbl(i) = kpblmin ELSE kbl(i) = x1 END IF END DO ! ! determine the level of maximum orographic height ! komax(:) = kbl(:) ! DO i=its,ite delks(i) = 1.0/(prsi(i, 1)-prsi(i, kbl(i))) delks1(i) = 1.0/(prsl(i, 1)-prsl(i, kbl(i))) END DO ! ! compute low level averages within pbl ! DO k=kts,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN rcsks = rcs*del(i, k)*delks(i) rdelks = del(i, k)*delks(i) ! pbl u mean ubar(i) = ubar(i) + rcsks*u1(i, k) ! pbl v mean vbar(i) = vbar(i) + rcsks*v1(i, k) ! ro mean roll(i) = roll(i) + rdelks*ro(i, k) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO END DO ! ! figure out low-level horizontal wind direction ! ! nwd 1 2 3 4 5 6 7 8 ! wd w s sw nw e n ne se ! DO i=its,ite wdir = ATAN2(ubar(i), vbar(i)) + pi idir = MOD(NINT(fdir*wdir), mdir) + 1 nwd = nwdir(idir) oa(i) = (1-2*INT((nwd-1)/4))*oa4(i, MOD(nwd-1, 4)+1) ol(i) = ol4(i, MOD(nwd-1, 4)+1) ! !----- compute orographic width along (ol) and perpendicular (olp) !----- the direction of wind ! ol4p(1) = ol4(i, 2) ol4p(2) = ol4(i, 1) ol4p(3) = ol4(i, 4) ol4p(4) = ol4(i, 3) olp(i) = ol4p(MOD(nwd-1, 4)+1) IF (ol(i) .LT. olmin) THEN max1 = olmin ELSE max1 = ol(i) END IF ! !----- compute orographic direction (horizontal orographic aspect ratio) ! od(i) = olp(i)/max1 IF (od(i) .GT. odmax) THEN od(i) = odmax ELSE od(i) = od(i) END IF IF (od(i) .LT. odmin) THEN od(i) = odmin ELSE od(i) = od(i) END IF ! !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions ! dxy(i) = dxy4(MOD(nwd-1, 4)+1) dxyp(i) = dxy4p(MOD(nwd-1, 4)+1) END DO ! !--- saving richardson number in usqj for migwdi ! DO k=kts,kte-1 DO i=its,ite ti = 2.0/(t1(i, k)+t1(i, k+1)) rdz = 1./(zl(i, k+1)-zl(i, k)) CALL PUSHREAL8(tem1) tem1 = u1(i, k) - u1(i, k+1) CALL PUSHREAL8(tem2) tem2 = v1(i, k) - v1(i, k+1) dw2 = rcl*(tem1*tem1+tem2*tem2) IF (dw2 .LT. dw2min) THEN CALL PUSHREAL8(max2) max2 = dw2min CALL PUSHCONTROL1B(0) ELSE CALL PUSHREAL8(max2) max2 = dw2 CALL PUSHCONTROL1B(1) END IF shr2 = max2*rdz*rdz CALL PUSHREAL8(bvf2) bvf2 = g*(g/cp+rdz*(vtj(i, k+1)-vtj(i, k)))*ti IF (bvf2/shr2 .LT. rimin) THEN usqj(i, k) = rimin CALL PUSHCONTROL1B(0) ELSE usqj(i, k) = bvf2/shr2 CALL PUSHCONTROL1B(1) END IF bnv2(i, k) = 2.0*g*rdz*(vtk(i, k+1)-vtk(i, k))/(vtk(i, k+1)+vtk(i& & , k)) IF (bnv2(i, k) .LT. bnv2min) THEN bnv2(i, k) = bnv2min CALL PUSHCONTROL1B(0) ELSE CALL PUSHCONTROL1B(1) bnv2(i, k) = bnv2(i, k) END IF END DO END DO ! !----compute the "low level" or 1/3 wind magnitude (m/s) ! DO i=its,ite x2 = SQRT(ubar(i)*ubar(i) + vbar(i)*vbar(i)) IF (x2 .LT. 1.0) THEN ulow(i) = 1.0 CALL PUSHCONTROL1B(0) ELSE ulow(i) = x2 CALL PUSHCONTROL1B(1) END IF rulow(i) = 1./ulow(i) END DO ! DO k=kts,kte-1 DO i=its,ite velco(i, k) = 0.5*rcs*((u1(i, k)+u1(i, k+1))*ubar(i)+(v1(i, k)+v1(& & i, k+1))*vbar(i)) CALL PUSHREAL8(velco(i, k)) velco(i, k) = velco(i, k)*rulow(i) IF (velco(i, k) .LT. veleps .AND. velco(i, k) .GT. 0.) THEN velco(i, k) = veleps CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO END DO ! ! no drag when critical level in the base layer ! DO i=its,ite ldrag(i) = velco(i, 1) .LE. 0. END DO ! ! no drag when velco.lt.0 ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. velco(i, k) .LE. 0. END DO END DO ! ! no drag when bnv2.lt.0 ! DO k=kts,kpblmax DO i=its,ite IF (k .LT. kbl(i)) ldrag(i) = ldrag(i) .OR. bnv2(i, k) .LT. 0. END DO END DO ! !-----the low level weighted average ri is stored in usqj(1,1; im) !-----the low level weighted average n**2 is stored in bnv2(1,1; im) !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of / ! DO i=its,ite wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i) CALL PUSHREAL8(bnv2(i, 1)) bnv2(i, 1) = wtkbj*bnv2(i, 1) CALL PUSHREAL8(usqj(i, 1)) usqj(i, 1) = wtkbj*usqj(i, 1) END DO ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i) tmp = bnv2(i, 1) + bnv2(i, k)*rdelks CALL PUSHREAL8(bnv2(i, 1)) bnv2(i, 1) = tmp tmp0 = usqj(i, 1) + usqj(i, k)*rdelks CALL PUSHREAL8(usqj(i, 1)) usqj(i, 1) = tmp0 CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO END DO ! DO i=its,ite ldrag(i) = ldrag(i) .OR. bnv2(i, 1) .LE. 0.0 ldrag(i) = ldrag(i) .OR. ulow(i) .EQ. 1.0 ldrag(i) = ldrag(i) .OR. var(i) .LE. 0.0 END DO ! ! set all ri low level values to the low level value ! DO k=kpblmin,kpblmax DO i=its,ite IF (k .LT. kbl(i)) THEN tmp1 = usqj(i, 1) CALL PUSHREAL8(usqj(i, k)) usqj(i, k) = tmp1 CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO END DO ! DO i=its,ite IF (.NOT.ldrag(i)) THEN bnv(i) = SQRT(bnv2(i, 1)) fr(i) = bnv(i)*rulow(i)*2.*var(i)*od(i) IF (fr(i) .GT. frmax) THEN fr(i) = frmax CALL PUSHCONTROL1B(0) ELSE CALL PUSHCONTROL1B(1) fr(i) = fr(i) END IF xn(i) = ubar(i)*rulow(i) yn(i) = vbar(i)*rulow(i) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO ! ! compute the base level stress and store it in taub ! calculate enhancement factor, number of mountains & aspect ! ratio const. use simplified relationship between standard ! deviation & critical hgt ! DO i=its,ite IF (.NOT.ldrag(i)) THEN CALL PUSHREAL8(efact) efact = (oa(i)+2.)**(ce*fr(i)/frc) IF (efact .LT. efmin) THEN x3 = efmin CALL PUSHCONTROL1B(0) ELSE x3 = efact CALL PUSHCONTROL1B(1) END IF IF (x3 .GT. efmax) THEN efact = efmax CALL PUSHCONTROL1B(0) ELSE efact = x3 CALL PUSHCONTROL1B(1) END IF !!!!!!! cleff (effective grid length) is highly tunable parameter !!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag cleff = SQRT(dxy(i)**2. + dxyp(i)**2.) IF (dxmeter .LT. cleff) THEN max3 = cleff ELSE max3 = dxmeter END IF cleff = 3.*max3 coefm(i) = (1.+ol(i))**(oa(i)+1.) xlinv(i) = coefm(i)/cleff tem = fr(i)*fr(i)*oc1(i) gfobnv = gmax*tem/((tem+cg)*bnv(i)) taub(i) = xlinv(i)*roll(i)*ulow(i)*ulow(i)*ulow(i)*gfobnv*efact CALL PUSHCONTROL1B(1) ELSE taub(i) = 0.0 xn(i) = 0.0 yn(i) = 0.0 CALL PUSHCONTROL1B(0) END IF END DO ! ! now compute vertical structure of the stress. ! DO k=kts,kpblmax DO i=its,ite IF (k .LE. kbl(i)) THEN taup(i, k) = taub(i) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO END DO ! ! vertical level k loop! DO k=kpblmin,kte-1 CALL PUSHINTEGER4(kp1) kp1 = k + 1 DO i=its,ite ! ! unstablelayer if ri < ric ! unstable layer if upper air vel comp along surf vel <=0 (crit lay) ! at (u-c)=0. crit layer exists and bit vector should be set (.le.) ! IF (k .GE. kbl(i)) THEN icrilv(i) = (icrilv(i) .OR. usqj(i, k) .LT. ric) .OR. velco(i, k& & ) .LE. 0.0 IF (bnv2(i, k) .LT. bnv2min) THEN CALL PUSHREAL8(brvf(i)) brvf(i) = bnv2min CALL PUSHCONTROL1B(0) ELSE CALL PUSHREAL8(brvf(i)) brvf(i) = bnv2(i, k) CALL PUSHCONTROL1B(1) END IF ! brunt-vaisala frequency CALL PUSHREAL8(brvf(i)) brvf(i) = SQRT(brvf(i)) CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF END DO ! DO i=its,ite IF (k .GE. kbl(i) .AND. (.NOT.ldrag(i))) THEN IF (.NOT.icrilv(i) .AND. taup(i, k) .GT. 0.0) THEN temv = 1.0/velco(i, k) CALL PUSHREAL8(tem1) tem1 = coefm(i)/dxy(i)*(ro(i, kp1)+ro(i, k))*brvf(i)*velco(i, & & k)*0.5 CALL PUSHREAL8(hd) hd = SQRT(taup(i, k)/tem1) fro = brvf(i)*hd*temv ! ! rim is the minimum-richardson number by shutts (1985) ! CALL PUSHREAL8(tem2) tem2 = SQRT(usqj(i, k)) tem = 1. + tem2*fro rim = usqj(i, k)*(1.-fro)/(tem*tem) ! ! check stability to employ the 'saturation hypothesis' ! of lindzen (1981) except at tropospheric downstream regions ! IF (rim .LE. ric) THEN ! saturation hypothesis! IF (oa(i) .LE. 0. .OR. kp1 .GE. kpblmin) THEN temc = 2.0 + 1.0/tem2 hd = velco(i, k)*(2.*SQRT(temc)-temc)/brvf(i) taup(i, kp1) = tem1*hd*hd CALL PUSHCONTROL3B(4) ELSE CALL PUSHCONTROL3B(3) END IF ELSE ! no wavebreaking! tmp2 = taup(i, k) taup(i, kp1) = tmp2 CALL PUSHCONTROL3B(2) END IF ELSE CALL PUSHCONTROL3B(1) END IF ELSE CALL PUSHCONTROL3B(0) END IF END DO END DO ! IF (lcap .LT. kte) THEN DO klcap=lcapp1,kte DO i=its,ite tmp3 = prsi(i, klcap)/prsi(i, lcap)*taup(i, lcap) CALL PUSHREAL8(taup(i, klcap)) taup(i, klcap) = tmp3 END DO END DO CALL PUSHCONTROL1B(1) ELSE CALL PUSHCONTROL1B(0) END IF DO i=its,ite IF (.NOT.ldrag(i)) THEN ! !------- determine the height of flow-blocking layer ! CALL PUSHINTEGER4(kblk) kblk = 0 pe = 0.0 DO k=kte,kpblmin,-1 IF (kblk .EQ. 0 .AND. k .LE. komax(i)) THEN pe = pe + bnv2(i, k)*(zl(i, komax(i))-zl(i, k))*del(i, k)/g/ro& & (i, k) ke = 0.5*((rcs*u1(i, k))**2.+(rcs*v1(i, k))**2.) ! !---------- apply flow-blocking drag when pe >= ke ! IF (pe .GE. ke) THEN CALL PUSHINTEGER4(kblk) kblk = k IF (kblk .GT. kbl(i)) THEN kblk = kbl(i) ELSE kblk = kblk END IF CALL PUSHREAL8(zblk) zblk = zl(i, kblk) - zl(i, kts) CALL PUSHCONTROL2B(2) ELSE CALL PUSHCONTROL2B(1) END IF ELSE CALL PUSHCONTROL2B(0) END IF END DO IF (kblk .NE. 0) THEN ! !--------- compute flow-blocking stress ! IF (2.0 - 1.0/od(i) .LT. 0.0) THEN CALL PUSHREAL8(cd) cd = 0.0 CALL PUSHCONTROL1B(0) ELSE CALL PUSHREAL8(cd) cd = 2.0 - 1.0/od(i) CALL PUSHCONTROL1B(1) END IF taufb(i, kts) = 0.5*roll(i)*coefm(i)/dxy(i)**2*cd*dxyp(i)*olp(i)& & *zblk*ulow(i)**2 tautem = taufb(i, kts)/FLOAT(kblk-kts) DO k=kts+1,kblk taufb(i, k) = taufb(i, k-1) - tautem END DO CALL PUSHINTEGER4(k - 1) ! !----------sum orographic GW stress and flow-blocking stress ! CALL PUSHREAL8ARRAY(taup(i, :), kte - kts + 2) taup(i, :) = taup(i, :) + taufb(i, :) CALL PUSHCONTROL2B(2) ELSE CALL PUSHCONTROL2B(1) END IF ELSE CALL PUSHCONTROL2B(0) END IF END DO ! ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy ! DO k=kts,kte DO i=its,ite taud(i, k) = 1.*(taup(i, k+1)-taup(i, k))*csg/del(i, k) END DO END DO ! ! limit de-acceleration (momentum deposition ) at top to 1/2 value ! the idea is some stuff must go out the 'top' ! DO klcap=lcap,kte DO i=its,ite taud(i, klcap) = taud(i, klcap)*factop END DO END DO ! ! if the gravity wave drag would force a critical line ! in the lower ksmm1 layers during the next deltim timestep, ! then only apply drag until that critical line is reached. ! DO k=kts,kpblmax-1 DO i=its,ite IF (k .LE. kbl(i)) THEN IF (taud(i, k) .NE. 0.) THEN x4 = velco(i, k)/(deltim*rcs*taud(i, k)) IF (x4 .GE. 0.) THEN y1 = x4 CALL PUSHCONTROL1B(0) ELSE y1 = -x4 CALL PUSHCONTROL1B(1) END IF IF (dtfac(i) .GT. y1) THEN dtfac(i) = y1 CALL PUSHCONTROL2B(2) ELSE dtfac(i) = dtfac(i) CALL PUSHCONTROL2B(3) END IF ELSE CALL PUSHCONTROL2B(1) END IF ELSE CALL PUSHCONTROL2B(0) END IF END DO END DO ! DO k=kts,kte DO i=its,ite CALL PUSHREAL8(taud(i, k)) taud(i, k) = taud(i, k)*dtfac(i) END DO END DO DO i=ite,its,-1 dvsfcb(i) = -(rcs*dvsfcb(i)/g) dusfcb(i) = -(rcs*dusfcb(i)/g) END DO taudb = 0.0 xnb = 0.0 ynb = 0.0 dtfacb = 0.0 DO k=kte,kts,-1 DO i=ite,its,-1 dtaux = taud(i, k)*xn(i) dtauy = taud(i, k)*yn(i) dtauyb = dvdtb(i, k) + dtauy2db(i, k) + del(i, k)*dvsfcb(i) delb(i, k) = delb(i, k) + dtaux*dusfcb(i) + dtauy*dvsfcb(i) dtauxb = dudtb(i, k) + dtaux2db(i, k) + del(i, k)*dusfcb(i) dtauy2db(i, k) = 0.0 dtaux2db(i, k) = 0.0 taudb(i, k) = taudb(i, k) + xn(i)*dtauxb + yn(i)*dtauyb ynb(i) = ynb(i) + taud(i, k)*dtauyb xnb(i) = xnb(i) + taud(i, k)*dtauxb CALL POPREAL8(taud(i, k)) dtfacb(i) = dtfacb(i) + taud(i, k)*taudb(i, k) taudb(i, k) = dtfac(i)*taudb(i, k) END DO END DO DO i=ite,its,-1 dvsfcb(i) = 0.0 dusfcb(i) = 0.0 END DO velcob = 0.0 DO k=kpblmax-1,kts,-1 DO i=ite,its,-1 CALL POPCONTROL2B(branch) IF (branch .GE. 2) THEN IF (branch .EQ. 2) THEN y1b = dtfacb(i) dtfacb(i) = 0.0 ELSE y1b = 0.0 END IF CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN x4b = y1b ELSE x4b = -y1b END IF temp12 = deltim*rcs*taud(i, k) velcob(i, k) = velcob(i, k) + x4b/temp12 taudb(i, k) = taudb(i, k) - velco(i, k)*deltim*rcs*x4b/temp12**2 END IF END DO END DO DO klcap=kte,lcap,-1 DO i=ite,its,-1 taudb(i, klcap) = factop*taudb(i, klcap) END DO END DO taupb = 0.0 DO k=kte,kts,-1 DO i=ite,its,-1 tempb20 = csg*taudb(i, k)/del(i, k) taupb(i, k+1) = taupb(i, k+1) + tempb20 taupb(i, k) = taupb(i, k) - tempb20 delb(i, k) = delb(i, k) - (taup(i, k+1)-taup(i, k))*tempb20/del(i& & , k) taudb(i, k) = 0.0 END DO END DO rollb = 0.0 zblkb = 0.0 taufbb = 0.0 ulowb = 0.0 DO i=ite,its,-1 CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN IF (branch .NE. 1) THEN CALL POPREAL8ARRAY(taup(i, :), kte - kts + 2) taufbb(i, :) = taufbb(i, :) + taupb(i, :) tautemb = 0.0 CALL POPINTEGER4(ad_to) DO k=ad_to,kts+1,-1 taufbb(i, k-1) = taufbb(i, k-1) + taufbb(i, k) tautemb = tautemb - taufbb(i, k) taufbb(i, k) = 0.0 END DO taufbb(i, kts) = taufbb(i, kts) + tautemb/FLOAT(kblk-kts) tempb18 = cd*0.5*olp(i)*coefm(i)*dxyp(i)*taufbb(i, kts) tempb19 = ulow(i)**2*tempb18/dxy(i)**2 rollb(i) = rollb(i) + zblk*tempb19 zblkb = zblkb + roll(i)*tempb19 ulowb(i) = ulowb(i) + roll(i)*zblk*2*ulow(i)*tempb18/dxy(i)**2 taufbb(i, kts) = 0.0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL8(cd) ELSE CALL POPREAL8(cd) END IF END IF DO k=kpblmin,kte,1 CALL POPCONTROL2B(branch) IF (branch .NE. 0) THEN IF (branch .NE. 1) THEN CALL POPREAL8(zblk) zlb(i, kblk) = zlb(i, kblk) + zblkb zlb(i, kts) = zlb(i, kts) - zblkb CALL POPINTEGER4(kblk) zblkb = 0.0 END IF END IF END DO CALL POPINTEGER4(kblk) END IF END DO CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN DO klcap=kte,lcapp1,-1 DO i=ite,its,-1 CALL POPREAL8(taup(i, klcap)) tmpb3 = taupb(i, klcap) taupb(i, klcap) = 0.0 tempb17 = tmpb3/prsi(i, lcap) prsib(i, klcap) = prsib(i, klcap) + taup(i, lcap)*tempb17 taupb(i, lcap) = taupb(i, lcap) + prsi(i, klcap)*tempb17 prsib(i, lcap) = prsib(i, lcap) - prsi(i, klcap)*taup(i, lcap)*& & tempb17/prsi(i, lcap) END DO END DO END IF brvfb = 0.0 bnv2b = 0.0 rob = 0.0 usqjb = 0.0 DO k=kte-1,kpblmin,-1 DO i=ite,its,-1 CALL POPCONTROL3B(branch) IF (branch .GE. 2) THEN IF (branch .EQ. 2) THEN tmpb2 = taupb(i, kp1) taupb(i, kp1) = 0.0 taupb(i, k) = taupb(i, k) + tmpb2 tem1b = 0.0 tem2b = 0.0 ELSE IF (branch .EQ. 3) THEN tem1b = 0.0 tem2b = 0.0 ELSE tem1b = hd**2*taupb(i, kp1) hdb = tem1*2*hd*taupb(i, kp1) taupb(i, kp1) = 0.0 temc = 2.0 + 1.0/tem2 temp11 = SQRT(temc) tempb16 = (2.*temp11-temc)*hdb/brvf(i) temp10 = velco(i, k)/brvf(i) velcob(i, k) = velcob(i, k) + tempb16 brvfb(i) = brvfb(i) - temp10*tempb16 IF (temc .EQ. 0.0) THEN temcb = -(temp10*hdb) ELSE temcb = (2.*temp10/(2.0*temp11)-temp10)*hdb END IF tem2b = -(temcb/tem2**2) END IF CALL POPREAL8(tem2) IF (.NOT.usqj(i, k) .EQ. 0.0) usqjb(i, k) = usqjb(i, k) + tem2b/& & (2.0*SQRT(usqj(i, k))) CALL POPREAL8(hd) CALL POPREAL8(tem1) temp9 = brvf(i)/dxy(i) tempb13 = coefm(i)*0.5*tem1b tempb14 = velco(i, k)*temp9*tempb13 tempb15 = (ro(i, kp1)+ro(i, k))*tempb13 rob(i, kp1) = rob(i, kp1) + tempb14 rob(i, k) = rob(i, k) + tempb14 velcob(i, k) = velcob(i, k) + temp9*tempb15 brvfb(i) = brvfb(i) + velco(i, k)*tempb15/dxy(i) END IF END DO DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN CALL POPREAL8(brvf(i)) IF (brvf(i) .EQ. 0.0) THEN brvfb(i) = 0.0 ELSE brvfb(i) = brvfb(i)/(2.0*SQRT(brvf(i))) END IF CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL8(brvf(i)) brvfb(i) = 0.0 ELSE CALL POPREAL8(brvf(i)) bnv2b(i, k) = bnv2b(i, k) + brvfb(i) brvfb(i) = 0.0 END IF END IF END DO CALL POPINTEGER4(kp1) END DO taubb = 0.0 DO k=kpblmax,kts,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN taubb(i) = taubb(i) + taupb(i, k) taupb(i, k) = 0.0 END IF END DO END DO bnvb = 0.0 frb = 0.0 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN ynb(i) = 0.0 xnb(i) = 0.0 taubb(i) = 0.0 ELSE tem = fr(i)*fr(i)*oc1(i) gfobnv = gmax*tem/((tem+cg)*bnv(i)) tempb10 = xlinv(i)*ulow(i)**3*taubb(i) ulowb(i) = ulowb(i) + roll(i)*gfobnv*efact*xlinv(i)*3*ulow(i)**2*& & taubb(i) rollb(i) = rollb(i) + gfobnv*efact*tempb10 gfobnvb = roll(i)*efact*tempb10 efactb = roll(i)*gfobnv*tempb10 taubb(i) = 0.0 temp8 = (cg+tem)*bnv(i) tempb11 = gmax*gfobnvb/temp8 tempb12 = -(tem*tempb11/temp8) temb = bnv(i)*tempb12 + tempb11 bnvb(i) = bnvb(i) + (cg+tem)*tempb12 frb(i) = frb(i) + oc1(i)*2*fr(i)*temb CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN x3b = 0.0 ELSE x3b = efactb END IF CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN efactb = 0.0 ELSE efactb = x3b END IF CALL POPREAL8(efact) IF (.NOT.oa(i) + 2. .LE. 0.0) frb(i) = frb(i) + ce*(oa(i)+2.)**(ce& & *(fr(i)/frc))*LOG(oa(i)+2.)*efactb/frc END IF END DO vbarb = 0.0 ubarb = 0.0 rulowb = 0.0 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN vbarb(i) = vbarb(i) + rulow(i)*ynb(i) rulowb(i) = rulowb(i) + ubar(i)*xnb(i) + vbar(i)*ynb(i) ynb(i) = 0.0 ubarb(i) = ubarb(i) + rulow(i)*xnb(i) xnb(i) = 0.0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) frb(i) = 0.0 tempb9 = var(i)*2.*od(i)*frb(i) bnvb(i) = bnvb(i) + rulow(i)*tempb9 rulowb(i) = rulowb(i) + bnv(i)*tempb9 frb(i) = 0.0 IF (.NOT.bnv2(i, 1) .EQ. 0.0) bnv2b(i, 1) = bnv2b(i, 1) + bnvb(i)/& & (2.0*SQRT(bnv2(i, 1))) bnvb(i) = 0.0 END IF END DO DO k=kpblmax,kpblmin,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN CALL POPREAL8(usqj(i, k)) tmpb1 = usqjb(i, k) usqjb(i, k) = 0.0 usqjb(i, 1) = usqjb(i, 1) + tmpb1 END IF END DO END DO delks1b = 0.0 DO k=kpblmax,kpblmin,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN tmpb0 = bnv2b(i, 1) rdelks = (prsl(i, k)-prsl(i, k+1))*delks1(i) CALL POPREAL8(usqj(i, 1)) tmpb = usqjb(i, 1) usqjb(i, 1) = tmpb usqjb(i, k) = usqjb(i, k) + rdelks*tmpb CALL POPREAL8(bnv2(i, 1)) rdelksb = bnv2(i, k)*tmpb0 + usqj(i, k)*tmpb bnv2b(i, 1) = tmpb0 bnv2b(i, k) = bnv2b(i, k) + rdelks*tmpb0 prslb(i, k) = prslb(i, k) + delks1(i)*rdelksb prslb(i, k+1) = prslb(i, k+1) - delks1(i)*rdelksb delks1b(i) = delks1b(i) + (prsl(i, k)-prsl(i, k+1))*rdelksb END IF END DO END DO DO i=ite,its,-1 wtkbj = (prsl(i, 1)-prsl(i, 2))*delks1(i) CALL POPREAL8(usqj(i, 1)) CALL POPREAL8(bnv2(i, 1)) wtkbjb = bnv2(i, 1)*bnv2b(i, 1) + usqj(i, 1)*usqjb(i, 1) usqjb(i, 1) = wtkbj*usqjb(i, 1) bnv2b(i, 1) = wtkbj*bnv2b(i, 1) prslb(i, 1) = prslb(i, 1) + delks1(i)*wtkbjb prslb(i, 2) = prslb(i, 2) - delks1(i)*wtkbjb delks1b(i) = delks1b(i) + (prsl(i, 1)-prsl(i, 2))*wtkbjb END DO DO k=kte-1,kts,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) velcob(i, k) = 0.0 CALL POPREAL8(velco(i, k)) rulowb(i) = rulowb(i) + velco(i, k)*velcob(i, k) velcob(i, k) = rulow(i)*velcob(i, k) tempb8 = rcs*0.5*velcob(i, k) u1b(i, k) = u1b(i, k) + ubar(i)*tempb8 u1b(i, k+1) = u1b(i, k+1) + ubar(i)*tempb8 ubarb(i) = ubarb(i) + (u1(i, k)+u1(i, k+1))*tempb8 v1b(i, k) = v1b(i, k) + vbar(i)*tempb8 v1b(i, k+1) = v1b(i, k+1) + vbar(i)*tempb8 vbarb(i) = vbarb(i) + (v1(i, k)+v1(i, k+1))*tempb8 velcob(i, k) = 0.0 END DO END DO DO i=ite,its,-1 ulowb(i) = ulowb(i) - rulowb(i)/ulow(i)**2 rulowb(i) = 0.0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN ulowb(i) = 0.0 x2b = 0.0 ELSE x2b = ulowb(i) ulowb(i) = 0.0 END IF IF (ubar(i)**2 + vbar(i)**2 .EQ. 0.0) THEN tempb7 = 0.0 ELSE tempb7 = x2b/(2.0*SQRT(ubar(i)**2+vbar(i)**2)) END IF ubarb(i) = ubarb(i) + 2*ubar(i)*tempb7 vbarb(i) = vbarb(i) + 2*vbar(i)*tempb7 END DO vtjb = 0.0 vtkb = 0.0 DO k=kte-1,kts,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) bnv2b(i, k) = 0.0 rdz = 1./(zl(i, k+1)-zl(i, k)) temp6 = vtk(i, k+1) + vtk(i, k) temp7 = vtk(i, k+1) - vtk(i, k) tempb5 = g*2.0*bnv2b(i, k)/temp6 tempb6 = -(rdz*temp7*tempb5/temp6) rdzb = temp7*tempb5 vtkb(i, k+1) = vtkb(i, k+1) + tempb6 + rdz*tempb5 vtkb(i, k) = vtkb(i, k) + tempb6 - rdz*tempb5 bnv2b(i, k) = 0.0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN usqjb(i, k) = 0.0 bvf2b = 0.0 shr2b = 0.0 ELSE shr2 = max2*rdz*rdz bvf2b = usqjb(i, k)/shr2 shr2b = -(bvf2*usqjb(i, k)/shr2**2) usqjb(i, k) = 0.0 END IF ti = 2.0/(t1(i, k)+t1(i, k+1)) CALL POPREAL8(bvf2) temp5 = vtj(i, k+1) - vtj(i, k) tempb4 = g*ti*bvf2b rdzb = rdzb + max2*2*rdz*shr2b + temp5*tempb4 vtjb(i, k+1) = vtjb(i, k+1) + rdz*tempb4 vtjb(i, k) = vtjb(i, k) - rdz*tempb4 tib = g*(g/cp+rdz*temp5)*bvf2b max2b = rdz**2*shr2b CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN CALL POPREAL8(max2) dw2b = 0.0 ELSE CALL POPREAL8(max2) dw2b = max2b END IF tem1b = rcl*2*tem1*dw2b tem2b = rcl*2*tem2*dw2b CALL POPREAL8(tem2) v1b(i, k) = v1b(i, k) + tem2b v1b(i, k+1) = v1b(i, k+1) - tem2b CALL POPREAL8(tem1) u1b(i, k) = u1b(i, k) + tem1b u1b(i, k+1) = u1b(i, k+1) - tem1b temp4 = zl(i, k+1) - zl(i, k) tempb2 = -(rdzb/temp4**2) zlb(i, k+1) = zlb(i, k+1) + tempb2 zlb(i, k) = zlb(i, k) - tempb2 temp3 = t1(i, k) + t1(i, k+1) tempb3 = -(2.0*tib/temp3**2) t1b(i, k) = t1b(i, k) + tempb3 t1b(i, k+1) = t1b(i, k+1) + tempb3 END DO END DO delksb = 0.0 DO k=kpblmax,kts,-1 DO i=ite,its,-1 CALL POPCONTROL1B(branch) IF (branch .NE. 0) THEN rdelks = del(i, k)*delks(i) rdelksb = ro(i, k)*rollb(i) rob(i, k) = rob(i, k) + rdelks*rollb(i) rcsks = rcs*del(i, k)*delks(i) rcsksb = u1(i, k)*ubarb(i) + v1(i, k)*vbarb(i) v1b(i, k) = v1b(i, k) + rcsks*vbarb(i) u1b(i, k) = u1b(i, k) + rcsks*ubarb(i) delb(i, k) = delb(i, k) + rcs*delks(i)*rcsksb + delks(i)*rdelksb delksb(i) = delksb(i) + rcs*del(i, k)*rcsksb + del(i, k)*rdelksb END IF END DO END DO DO i=ite,its,-1 temp2 = prsl(i, 1) - prsl(i, kbl(i)) tempb0 = -(delks1b(i)/temp2**2) prslb(i, 1) = prslb(i, 1) + tempb0 prslb(i, kbl(i)) = prslb(i, kbl(i)) - tempb0 delks1b(i) = 0.0 temp1 = prsi(i, 1) - prsi(i, kbl(i)) tempb1 = -(delksb(i)/temp1**2) prsib(i, 1) = prsib(i, 1) + tempb1 prsib(i, kbl(i)) = prsib(i, kbl(i)) - tempb1 delksb(i) = 0.0 END DO DO k=kte,kts,-1 DO i=ite,its,-1 tempb = vtkb(i, k)/prslk(i, k) temp0 = rd*vtj(i, k) prslb(i, k) = prslb(i, k) + rob(i, k)/temp0 vtjb(i, k) = vtjb(i, k) + tempb - prsl(i, k)*rd*rob(i, k)/temp0**2 rob(i, k) = 0.0 prslkb(i, k) = prslkb(i, k) - vtj(i, k)*tempb/prslk(i, k) vtkb(i, k) = 0.0 t1b(i, k) = t1b(i, k) + (fv*q1(i, k)+1.)*vtjb(i, k) q1b(i, k) = q1b(i, k) + t1(i, k)*fv*vtjb(i, k) vtjb(i, k) = 0.0 END DO END DO DO k=kte,kts,-1 DO i=ite,its,-1 dtauy2db(i, k) = 0.0 dtaux2db(i, k) = 0.0 END DO END DO END SUBROUTINE GWDO2D_B end module a_module_bl_gwdo