68 subroutine gwdc_run (im,km,lat,u1,v1,t1,q1,deltim, &
69 & pmid1,pint1,dpmid1,qmax,ktop,kbot,kcnv,cldf, &
70 & grav,cp,rd,fv,pi,dlength,lprnt,ipr,fhour, &
71 & utgwc,vtgwc,tauctx,taucty,errmsg,errflg)
86 USE machine ,
ONLY : kind_phys
112 integer,
intent(in) :: im, km, lat, ipr
113 integer,
intent(in) :: ktop(:),kbot(:),kcnv(:)
114 real(kind=kind_phys),
intent(in) :: grav,cp,rd,fv,fhour,deltim,pi
115 real(kind=kind_phys),
dimension(:),
intent(in) :: qmax
116 real(kind=kind_phys),
dimension(:),
intent(out) :: tauctx,taucty
117 real(kind=kind_phys),
dimension(:),
intent(in) :: cldf,dlength
118 real(kind=kind_phys),
dimension(:,:),
intent(in) :: u1,v1,t1, &
120 real(kind=kind_phys),
dimension(:,:),
intent(out) :: utgwc,vtgwc
121 real(kind=kind_phys),
dimension(:,:),
intent(in) :: pint1
123 logical,
intent(in) :: lprnt
125 character(len=*),
intent(out) :: errmsg
126 integer,
intent(out) :: errflg
181 integer i,ii,k,k1,kk,kb,ilev,npt,kcb,kcldm,npr
182 integer,
dimension(im) :: ipt
184 real(kind=kind_phys) tem, tem1, tem2, qtem, wtgwc, tauct, &
185 & windcltop, shear, nonlinct, nonlin, nonlins,&
186 & n2, dtdp, crit1, crit2, p1, p2, &
191 integer,
allocatable :: kcldtop(:),kcldbot(:)
192 logical,
allocatable :: do_gwc(:)
193 real(kind=kind_phys),
allocatable :: tauctxl(:), tauctyl(:),
194 & gwdcloc(:), break(:),
197 & cosphi(:), sinphi(:),
198 & xstress(:), ystress(:),
199 & ucltop(:), vcltop(:),
201 & dlen(:), gqmcldlen(:)
206 real(kind=kind_phys),
allocatable :: plnint(:,:), velco(:,:),
207 & taugwci(:,:), bruni(:,:),
208 & rhoi(:,:), basicui(:,:),
209 & ti(:,:), riloc(:,:),
210 & rimin(:,:), pint(:,:)
212 real(kind=kind_phys),
allocatable ::
214 & plnmid(:,:), taugw(:,:),
215 & utgwcl(:,:), vtgwcl(:,:),
216 & basicum(:,:), u(:,:),v(:,:),
218 & pmid(:,:), dpmid(:,:),
220 & brunm(:,:), rhom(:,:)
246 real(kind=kind_phys),
parameter ::
247 & c1=1.41, c2=-0.38, ricrit=0.25
248 &, n2min=1.e-32, zero=0.0, one=1.0
249 &, taumin=1.0e-20, tauctmax=-20.
251 &, qmin=1.0e-10, shmin=1.0e-20
252 &, rimax=1.0e+20, rimaxm=0.99e+20
253 &, rimaxp=1.01e+20, rilarge=0.9e+20
254 &, riminx=-1.0e+20, riminm=-1.01e+20
255 &, riminp=-0.99e+20, rismall=-0.9e+20
266 if (kcnv(i) /= 0 .and. qmax(i) > zero)
then
336 allocate (kcldtop(npt), kcldbot(npt), do_gwc(npt))
337 allocate (tauctxl(npt), tauctyl(npt), dtfac(npt),
338 & gwdcloc(npt), break(npt), cosphi(npt),
340 & sinphi(npt), xstress(npt), ystress(npt), wrk(npt),
341 & ucltop(npt), vcltop(npt),dlen(npt), gqmcldlen(npt))
346 allocate (plnint(npt,2:km+1),
347 & taugwci(npt,km+1), bruni(npt,km+1),
348 & rhoi(npt,km+1), basicui(npt,km+1),
349 & ti(npt,km+1), riloc(npt,km+1),
350 & rimin(npt,km+1), pint(npt,km+1))
355 & (plnmid(npt,km), velco(npt,km),
356 & utgwcl(npt,km), vtgwcl(npt,km),
357 & basicum(npt,km), u(npt,km), v(npt,km),
358 & t(npt,km), spfh(npt,km), pmid(npt,km),
359 & dpmid(npt,km), taugw(npt,km),
361 & brunm(npt,km), rhom(npt,km))
373 if (ipr == ipt(i))
then
387 spfh(i,k) = max(q1(ii,k1),qmin)
388 pmid(i,k) = pmid1(ii,k1)
389 dpmid(i,k) = dpmid1(ii,k1) * onebg
392 rhom(i,k) = pmid(i,k) / (rd*t(i,k)*(1.0+fv*spfh(i,k)))
393 plnmid(i,k) = log(pmid(i,k))
407 pint(i,k) = pint1(ii,k1)
419 plnint(i,k) = log(pint(i,k))
425 kcldtop(i) = km - ktop(ii) + 1
426 kcldbot(i) = km - kbot(ii) + 1
427 dlen(i) = dlength(ii)
429 gqmcldlen(i) = grav*qmax(ii)*cldf(ii)*dlen(i)
502 rhoi(i,1) = pint(i,1)/(rd*ti(i,1))
503 bruni(i,1) = sqrt( gsqr / (cp*ti(i,1)) )
510 rhoi(i,km+1) = pint(i,km+1)/(rd*ti(i,km+1)*(1.0+fv*spfh(i,km)))
511 bruni(i,km+1) = sqrt( gsqr / (cp*ti(i,km+1)) )
524 tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
526 ti(i,k) = t(i,k-1) * tem1 + t(i,k) * tem2
527 qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2
528 rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) )
529 dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1))
530 n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp )
531 bruni(i,k) = sqrt(max(n2min, n2))
544 dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k))
545 n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp )
546 brunm(i,k) = sqrt(max(n2min, n2))
603 kcldm = max(kcldm,kk)
616 windcltop = 1.0 / sqrt( ucltop(i)*ucltop(i)
617 & + vcltop(i)*vcltop(i) )
618 cosphi(i) = ucltop(i)*windcltop
619 sinphi(i) = vcltop(i)*windcltop
637 basicum(i,k) = u(i,k)*cosphi(i) + v(i,k)*sinphi(i)
652 basicui(i,1) = basicum(i,1)
653 basicui(i,km+1) = basicum(i,km)
657 tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
659 basicui(i,k) = basicum(i,k)*tem2 + basicum(i,k-1)*tem1
693 shear = grav*rhoi(i,k) * (basicum(i,k) - basicum(i,k-1))
694 & / (pmid(i,k) - pmid(i,k-1))
695 if ( abs(shear) < shmin )
then
698 tem = bruni(i,k) / shear
699 riloc(i,k) = tem * tem
700 if (riloc(i,k) >= rimax ) riloc(i,k) = rilarge
706 riloc(i,1) = riloc(i,2)
707 riloc(i,km+1) = riloc(i,km)
824 if ( abs(basicui(i,kk)) > zero .and. riloc(i,kk) > ricrit)
then
828 nonlinct = gqmcldlen(i) / (bruni(i,kk)*t(i,kk)*tem1)
831 tauct = - rhom(i,kk) * tem * tem1 * c1 * tem2 * tem2
832 & / (bruni(i,kk)*dlen(i))
836 tauct = max(tauctmax, tauct)
837 tauctxl(i) = tauct * cosphi(i)
838 tauctyl(i) = tauct * sinphi(i)
839 taugwci(i,kk) = tauct
900 tem1 = (u(i,k)+u(i,k-1))*0.5
901 tem2 = (v(i,k)+v(i,k-1))*0.5
902 crit1 = ucltop(i)*tem1
903 crit2 = vcltop(i)*tem2
904 velco(i,k) = tem1 * cosphi(i) + tem2 * sinphi(i)
906 crit1 = ucltop(i)*u(i,1)
907 crit2 = vcltop(i)*v(i,1)
908 velco(i,1) = u(i,1) * cosphi(i) + v(i,1) * sinphi(i)
913 if ( abs(basicui(i,k)) > zero .and. crit1 > zero
914 & .and. crit2 > zero )
then
915 tem = basicui(i,k) * basicui(i,k)
916 nonlin = gqmcldlen(i) / (bruni(i,k)*ti(i,k)*tem)
918 if ( riloc(i,k) < rimaxm )
then
919 tem1 = 1 + tem*sqrt(riloc(i,k))
920 rimin(i,k) = riloc(i,k) * (1-tem) / (tem1*tem1)
921 else if((riloc(i,k) > rimaxm) .and.
922 & (riloc(i,k) < rimaxp))
then
923 rimin(i,k) = ( 1 - tem) / (tem*tem)
925 if ( rimin(i,k) <= riminx )
then
1007 if (k < kk .and. k > 1)
then
1008 if ( abs(taugwci(i,k+1)) > taumin )
then
1009 if ( riloc(i,k) > ricrit )
then
1010 if ( rimin(i,k) > ricrit )
then
1011 taugwci(i,k) = taugwci(i,k+1)
1012 elseif (rimin(i,k) > riminp)
then
1013 tem = 2.0 + 1.0 / sqrt(riloc(i,k))
1014 nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem)
1016 tem2 = c2*nonlins*tem1
1017 taugwci(i,k) = - rhoi(i,k) * c1 * tem1 * tem2 * tem2
1018 & / (bruni(i,k)*dlen(i))
1019 elseif (rimin(i,k) > riminm)
then
1035 if ( (basicum(i,k+1)*basicum(i,k) ) < 0. )
then
1036 taugwci(i,k+1) = zero
1040 if (abs(taugwci(i,k)) > abs(taugwci(i,k+1)))
then
1041 taugwci(i,k) = taugwci(i,k+1)
1044 elseif (k == 1)
then
1049 taugwci(i,1) = taugwci(i,2)
1072 taugw(i,k) = (taugwci(i,k+1) - taugwci(i,k)) / dpmid(i,k)
1073 if (taugw(i,k) /= 0.0)
then
1074 tem = deltim * taugw(i,k)
1075 dtfac(i) = min(dtfac(i), abs(velco(i,k)/tem))
1097 wtgwc = taugw(i,k) * dtfac(i)
1098 utgwcl(i,k) = wtgwc * cosphi(i)
1099 vtgwcl(i,k) = wtgwc * sinphi(i)
1128 xstress(i) = xstress(i) + utgwcl(i,k)*dpmid(i,k)
1129 ystress(i) = ystress(i) + vtgwcl(i,k)*dpmid(i,k)
1149 wrk(i) = 0.5 * pi / (pint(i,kcldbot(i)+1)-pint(i,kcldtop(i)))
1156 if (k >= kk .and. k <= kcldbot(i))
then
1157 p1 = sin(wrk(i) * (pint(i,k) -pint(i,kk)))
1158 p2 = sin(wrk(i) * (pint(i,k+1)-pint(i,kk)))
1159 tem = - (p2-p1) / dpmid(i,k)
1160 utgwcl(i,k) = tem*xstress(i)
1161 vtgwcl(i,k) = tem*ystress(i)
1311 utgwc(ii,k1) = utgwcl(i,k)
1313 vtgwc(ii,k1) = vtgwcl(i,k)
1326 tauctx(ii) = tauctxl(i)
1327 taucty(ii) = tauctyl(i)
1345 deallocate (kcldtop,kcldbot,do_gwc)
1346 deallocate (tauctxl, tauctyl, dtfac,
1348 & gwdcloc, break, cosphi,
1349 & sinphi, xstress, ystress,
1350 & dlen, ucltop, vcltop, gqmcldlen, wrk)
1352 deallocate (plnint, taugwci, velco,
1353 & bruni, rhoi, basicui,
1354 & ti, riloc, rimin, pint)
1356 deallocate (plnmid, utgwcl, vtgwcl, basicum, u, v, t,
1357 & pmid, dpmid, brunm, rhom, taugw)