subroutine VERTDF (akm, akh, el, fxh, fxe, fxmx, fxmy, & vfu, vfv, vft, vfr, cdm, cdm2,xxfh) !----------------------------------------------------------------------- ! ! VERTDF computes vertical diffusion using a mrf non-local mixing ! !----------------------------------------------------------------------- use allocate_mod include 'RESOLUTION.h' include 'PARAMETERS.h' include 'AREAS.h' include 'CLCLE.h' include 'CONMLEV.h' include 'CONSLEV.h' include 'ESTAB.h' include 'FILEC.h' include 'FILEPC.h' include 'FLAGS.h' include 'GDINFO.h' include 'TKE2P5.h' include 'LIMIT.h' include 'QLOGS.h' include 'RCON.h' include 'SET.h' include 'STORJ.h' include 'TIME.h' include 'VDIFS.h' include 'XMINI.h' c common /shali/ pblh(imx),heats(imx),evaps(imx) parameter (b1=16.6) !----------------------------------------------------------------------- ! user interface variables !----------------------------------------------------------------------- real, intent (in), dimension (i1:i2,kmax) :: akm real, intent (in), dimension (i1:i2,kmax) :: akh real, intent (in), dimension (ibeg:iend) :: fxh real, intent (in), dimension (ibeg:iend) :: fxe real, intent (in), dimension (ibeg:iend) :: fxmx real, intent (in), dimension (ibeg:iend) :: fxmy real, intent (out), dimension (kmax,ibeg:iend) :: vfu real, intent (out), dimension (kmax,ibeg:iend) :: vfv real, intent (out), dimension (kmax,ibeg:iend) :: vft real, intent (out), dimension (kmax,ibeg:iend) :: vfr !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! internal variables !----------------------------------------------------------------------- real, dimension(i1:i2,1) :: el0 logical, dimension(i1:i2,1) :: flagsea real, dimension(kmax,i1:i2) :: tv real, dimension(i1:i2,kmax) :: theta real, dimension(i1:i2,kmax) :: zfull real, dimension(i1:i2,kmax+1) :: zhalf real, dimension(i1:i2,kmax+1) :: el !----------------------------------------------------------------------- real karman real, dimension(imx) :: cdm real, dimension(imx) :: cdm2 real, dimension(imx) :: cht real, dimension(imx) :: xxfh real, dimension(imx) :: fluxw real, dimension(imx) :: shfx real, dimension(i1:i2) :: tab1 real, dimension(i1:i2) :: tab2 real, dimension(i1:i2) :: t1 real, dimension(i1:i2) :: t2 real, dimension(i1:i2) :: estso real, dimension(i1:i2) :: rstso integer, dimension(i1:i2) :: it ! znt used but not needed........ real, dimension(i1:i2,1) :: znt ! zol not needed..... real, dimension(i1:i2,1) :: zol ! regime not needed.... real, dimension(i1:i2,1) :: regime real, dimension(i1:i2,1) :: hol real, dimension(i1:i2,1) :: pbl real, dimension(i1:i2,1) :: xland ! tsk used but not needed........ real, dimension(i1:i2,1) :: gz1oz0 real, dimension(i1:i2,1) :: ust real, dimension(i1:i2,1) :: hfx real, dimension(i1:i2,1) :: qfx real, dimension(i1:i2,1) :: psih real, dimension(kmax ) :: u1d real, dimension(kmax ) :: v1d real, dimension(kmax ) :: t1d real, dimension(kmax ) :: q1d real, dimension(kmax ) :: ublten real, dimension(kmax ) :: vblten real, dimension(kmax ) :: tblten real, dimension(kmax ) :: qblten real, dimension(kmax ) :: sl real, dimension(kmax ) :: slk real, dimension(kmax ) :: del real, dimension(kmax+1 ) :: si real, dimension(kmax,imx) :: dish real, dimension(kmax,imx) :: bbbb real, dimension(i1:i2,kmax,1) :: z_ikj real,parameter :: atet=3.8e+2, btet=17.27, ctet=273., dtet=36.0 real,parameter :: ES0 = 611.2 ! !----------------------------------------------------------------------- ! internal variables with storage association !----------------------------------------------------------------------- integer, dimension (4,jmx) :: iflxmx integer, dimension (4,jmx) :: ishfmx real, dimension(4,jmx) :: flxmx real, dimension(4,jmx) :: shfmx equivalence (iflxmx(1,1),flxmx(1,1)) equivalence (ishfmx(1,1),shfmx(1,1)) equivalence (iflxmx,fmxj), (ishfmx,shmxj) ! ! ! determine whether to add in dissipative heating ! idish = 1 ! !----------------------------------------------------------------------- ! code !----------------------------------------------------------------------- ! r gas constant for dry air (j/kg/k) rmks = 287. ! rv gas constant for water vapor (j/kg/k) rvmks = 461.6 ! cp heat capacity at constant pressure for dry air(j/kg/k) cpmks = 7.*rmks/2.0 rovcp = rmks/cpmks ! g acceleration due to gravity (m/s^2 gmks = 9.81 rovg = rmks/gmks ! rvovrd rv divided by r (dimensionless) rvovrd = rvmks/rmks ! eomeg angular velocity of earth's rotation (rad/s) eomeg = 7.2921e-5 ! eomeg not needed... ! p_qi species index for cloud ice (dimensionless) p_qi = 0 p_index = 1 svp1 = es0*1.0e-03 svp2 = btet svp3 = dtet svpt0 = ctet ! ep1 constant for virtual temperature ep1 = rvmks/rmks - 1.0 ! ep2 constant for specific humidity ep2 = rmks/rvmks ! karman von karman constant karman = 0.4 ! xlv latent heat of vaporization (j/kg) xlvmks = 2.5e6 dw = 1.0 ! dtmin time step in minutes/real,pa dtmin = 60.*deltx ! a and c = constants used in evaluating universal function c = .76 j = IFIX(tjloc(i1)) do i = i1,i2 shfl(i,j) = 0.0 flxw(i,j) = 0.0 stpc(i) = stpc(i) + deltat enddo !----------------------------------------------------------------------- !----------------------------------------------------------------------- do i = i1,i2 do k = 1,kmax tv(k,i) = tpc(k,i ) *(1.+0.61*rpc(k,i,1)) enddo enddo call height (tv, pspc, zc, q, qmh, kmax, i1,i2,imx,zhalf, zfull) call thet (tv, pspc, zc, q, qmh, kmax, i1,i2,imx, theta) ! adjust units.. geopotential (cgs) to heights(m) do i = i1,i2 psfc =pspc(i) * 1.0e-4 znt(i,1)=abs(zoc(i))/100. tsk =tstrc(i) tab1 (i) = tstrc(i) - 153.16 it (i) = IFIX(tab1(i)) tab2 (i) = tab1(i) - FLOAT(it(i)) t1 (i) = tab(it(i) + 1) t2 (i) = table(it(i) + 1) estso(i) = t1(i) + tab2(i)*t2(i) psps1 = pspc(i) - estso(i) if(psps1 .EQ. 0.0)then psps1 = .1 endif rstso(i) = 0.622*estso(i)/psps1 qss = rstso(i)/(1.+rstso(i)) vrts = 1. + 0.61*wetc(i)*rstso(i) dthv = theta(i,kmax) - tstrc(i)*vrts wspd = SQRT(uc(kmax,i)*uc(kmax,i) + vc(kmax,i)*vc(kmax,i)) wspd = amax1(wspd ,100.) ! br =0.0001*zfull(i,kmax)*dthv/ ! & (gmks*theta(i,kmax)*wspd *wspd ) thetkm = tpc(kmax,i)*rqc9 vrtkx = 1.0 + 0.61*rpc (kmax,i,1) zkmax = rgas*tpc(kmax,i)*qqlog(kmax)*og tempa1 = thetkm *vrtkx - tstrc(i)*vrts tempa2 = tempa1 *(thetkm - tstrc(i)) if (tempa2 .LT. 0.) tempa1 = 1.0e-4 tab1(i) = ABS(tempa1 ) if (tab1(i) .LT. 1.0e-4) tempa1 = 1.0e-4 !------------------------------------------------------------------------ ! compute bulk richardson number "br" at each point. if "rb" ! exceeds 95% of critical richardson number "brcr" then "rb = brcr" !------------------------------------------------------------------------ br = ggg*zkmax *tempa1 / & (tpc(kmax,i)*vrtkx *wspd *wspd ) brcr = 0.95/(c*(1. - ABS(zoc(i))/zkmax )) if (br .GT. brcr ) br = brcr !------------------------------------------------------------------------ wspd = wspd/100. pi = 4.*atan(1.0) pi180 = 180 / pi tlat = 10.*rlatpc(i) * pi180 tlon = 3600. - 10.*rlonpc(i) * pi180 ilat = ifix(tlat) ilon = ifix(tlon) if (zoc(i) .GT. 0.0) then xland(i,1) = 2.0 else xland(i,1) = 1.0 endif gz1oz0(i,1)=alog(0.0001*zfull(i,kmax)/(gmks*znt(i,1))) ust (i,1)= 0.01*sqrt(cdm(i)* & (uc(kmax,i)*uc(kmax,i) + vc(kmax,i)*vc(kmax,i))) qfx (i,1)=10.*fxe(i) hfx (i,1)=0.001*cp*fxh(i) fm = karman/sqrt(cdm(i)) ! fm = karman/sqrt(cdm2(i)) fh = karman*xxfh(i) ! cht(i) = karman*karman/(fm * fh) cht(i) = karman*sqrt(cdm(i))/fh ch = cht(i) cd = cdm(i) wet = wetc(i) ! initial surface tendencies... DUSFC = 0.0 DVSFC = 0.0 DTSFC = 0.0 DQSFC = 0.0 si (1) = qmh(kmax+1) do k = 1,kmax u1d(kmax+1-k) = upc(k,i)/100. v1d(kmax+1-k) = vpc(k,i)/100. t1d (kmax+1-k) = tpc (k,i) q1d(kmax+1-k) = rpc (k,i,1)/(1. + rpc (k,i,1)) sl (kmax+1-k) = q (k) slk(kmax+1-k) = q (k) ** rovcp del(kmax+1-k) = dq (k) si (kmax+2-k) = qmh(k) ! initialize boundary layer tendencies ublten(kmax+1-k) = 0.0 vblten(kmax+1-k) = 0.0 tblten(kmax+1-k) = 0.0 qblten(kmax+1-k) = 0.0 enddo rcl=1.0 alpha = .60 call moninp(kmax,1,vblten,ublten,tblten,qblten, & u1d,v1d,t1d,q1d, & psfc,br ,CD,CH,FM,FH,tsk,QSS,wet,wspd,KPBL, & SI,DEL,SL,SLK,RCL,deltat,THOUR, & DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,HGAMQ, & heat, evap,alpha) !----------------------------------------------------------------------- ! convert momentum back to cgs and specific hum to mixing ratio !----------------------------------------------------------------------- do k = 1,kmax vfu(kmax+1-k,i) = 100.*pspc(i)*ublten(k) vfv(kmax+1-k,i) = 100.*pspc(i)*vblten(k) vft(kmax+1-k,i) = pspc(i)*tblten(k) vfr(kmax+1-k,i) = pspc(i)*qblten(k)* & (1.+rpc(kmax+1-k,i,1))/(1.-q1d(k)) enddo c c store the boundary layer height and surface heat for shallow convection c pblh(i) = HPBL heats(i)= heat evaps(i) = evap c enddo !----------------------------------------------------------------------- ! add in dissipative heating !----------------------------------------------------------------------- ! if(idish .EQ. 1) then !----------------------------------------------------------------------- ! add in dissipative heating ! ! determine if you should use the entire dissapative heating tendency ! ccc factd = .55 ! c c high res test c factd = .8 !----------------------------------------------------------------------- ! do i=i1,i2 do k = 1,kmax dish(k,i) = upc(k,i)*vfu(k,i)+ vpc(k,i)*vfv(k,i) enddo enddo do i=i1,i2 do k = 1,kmax cpinv = cp*(1.+0.8*rpc(k,i,1)) dish(k,i) = -dish(k,i)/cpinv enddo enddo c do i=i1,i2 do k = 1,kmax vft(k,i) = vft(k,i) + factd * dish(k,i) enddo enddo endif ! !----------------------------------------------------------------------- ! add to the integrals of accumulated heat flux and evaporation !----------------------------------------------------------------------- do i = i1,i2 fluxw(i) = - dw*fxe (i)*deltat shfx (i) = - cpp*fxh(i)*deltat enddo if (inc .LE. 1) go to 510 do i = i1,i2,inc do ii = 2,inc iii = ii + i - 1 shfx (iii) = -1.e10 fluxw(iii) = -1.e10 enddo enddo 510 continue do i = i1,i2,inc ccc acshc(i) = acshc(i) + shfx (i) cccc acevc(i) = acevc(i) + fluxw(i) <<<<<< store large scale precip here now enddo if (inwrb .EQ. 0) go to 560 !----------------------------------------------------------------------- ! if it is time to compute the integrals, the area averaged heat ! flux and evaporation will now be obtained !----------------------------------------------------------------------- do i = i1,i2,inc if (j .LE. js .OR. j .GE. jn) go to 540 if (i .LE. iw .OR. i .GE. ie) go to 540 fluxw(i) = -1.e10 shfx(i) = -1.e10 go to 550 540 continue areams(i,j,4,ngd) = areams(i,j,4,ngd) + & dlsc(i)*fluxw(i)*deltat/dstp areams(i,j,5,ngd) = areams(i,j,5,ngd) + & dlsc(i)*shfx (i)*deltat/dstp shfl(i,j) = shfx (i)/dstp flxw(i,j) = fluxw(i)/dstp 550 continue enddo 560 continue return end subroutine VERTDF