subroutine convtrans(igrid,dtinp,ncol,nrow,nlay,deltat,idfin, & depth,press,temp,cigfrc,cigtim,cigent,cigdet, & cigmas) use filunit c c----CAMx v7Beta6 190902 c c CONVTRANS integrates cloud-in-grid mass transport matrices based on input c K-F entrainment/detrainment fluxes. Assuming a steady-state sub-grid c cloud environment, these CiG matrices are calculated for each convective c column over an entire grid for the duration of the column-specific c convective time scale and then used in CIGDRIVE for all timesteps between c the meteorological update times (usually 1 hour). c c Copyright 1996 - 2018 c Ramboll c c Modifications: c None c c Input arguments: c igrid grid index c dtinp met input frequency (min) c ncol number of columns c nrow number of rows c nlay number of layers c deltat time step (s) c idfin map of nested grids in this grid c depth layer depth field (m) c press pressure field (mb) c temp temperature field (K) c cigfrc CiG cloud fraction field c cigtim CiG cloud timescale (s) c cigent CiG sub-grid cloud entrainment rate (kg/m2/s) c cigdet CiG sub-grid cloud detrainment rate (kg/m2/s) c c Output arguments: c cigmas CiG transport matrix (unitless) c c Routines called: c None c c Called by: c TSTEP_INIT c implicit none include "camx.prm" integer igrid,ncol,nrow,nlay integer idfin(ncol,nrow) real deltat,dtinp real depth(ncol,nrow,nlay) real press(ncol,nrow,nlay) real temp(ncol,nrow,nlay) real cigfrc(ncol,nrow) real cigtim(ncol,nrow) real cigent(ncol,nrow,nlay) real cigdet(ncol,nrow,nlay) real cigmas(2,nlay,nlay,ncol,nrow) integer i,j,k,kk,l,n,ntimes,ks,ke real timescale,dt,dcflx,daflx real fm,fp,crat real tmassi,tmassf,terr real dz(MXLAYER) real rho(MXLAYER) real entc(MXLAYER),detc(MXLAYER) real enta(MXLAYER),deta(MXLAYER) real cflx(MXLAYER),aflx(MXLAYER) real aatrac(MXLAYER,MXLAYER),actrac(MXLAYER,MXLAYER) real cctrac(MXLAYER,MXLAYER),catrac(MXLAYER,MXLAYER) real aanew(MXLAYER,MXLAYER),acnew(MXLAYER,MXLAYER) real ccnew(MXLAYER,MXLAYER),canew(MXLAYER,MXLAYER) c c-----Entry point c do j = 1,nrow do i = 1,ncol do l = 1,nlay do k = 1,nlay do n = 1,2 cigmas(n,k,l,i,j) = 0. enddo enddo enddo enddo enddo c c-----Loop over vertical grid columns, skipping any covered by a nest or c without active sub-grid convection c do 20 j = 2,nrow-1 do 10 i = 2,ncol-1 if (idfin(i,j).gt.igrid) goto 10 if (cigfrc(i,j).eq.0. .or. cigtim(i,j).eq.0.) goto 10 c c-----Set hydrostatic density profiles (kg/m3) and entrainment/detrainment c rates (kg/m2/s) c timescale = min(cigtim(i,j),60.*dtinp) crat = cigfrc(i,j)/(1. - cigfrc(i,j)) do k = 1,nlay dz(k) = depth(i,j,k) rho(k) = 100.*press(i,j,k)/(287.*temp(i,j,k)) entc(k) = cigent(i,j,k) detc(k) = cigdet(i,j,k) enta(k) = entc(k)*crat deta(k) = detc(k)*crat enddo c c-----Calculate in-cloud and ambient net vertical fluxes (kg/m2/s) c cflx(1) = entc(1) - detc(1) aflx(1) = deta(1) - enta(1) do k = 2,nlay cflx(k) = cflx(k-1) + entc(k) - detc(k) aflx(k) = aflx(k-1) - enta(k) + deta(k) enddo c c-----Initialize tracer matrix c do l = 1,nlay do k = 1,nlay cctrac(k,l) = 0. aatrac(k,l) = 0. catrac(k,l) = 0. actrac(k,l) = 0. enddo cctrac(l,l) = 1. aatrac(l,l) = 1. enddo c c-----Check time step, determine integration step for first-order upstream c transport technique c dt = timescale 100 continue do k = 1,nlay dcflx = cflx(k) daflx = aflx(k) if (k.gt.1) then dcflx = max(abs(cflx(k-1)),abs(cflx(k))) daflx = max(abs(aflx(k-1)),abs(aflx(k))) endif if ((dcflx + detc(k))*dt/dz(k).gt.rho(k)/2. .or. & (daflx + enta(k))*dt/dz(k).gt.rho(k)/2.) then dt = dt/2. goto 100 endif enddo c c-----Loop for the duration of the input timestep c ntimes = nint(timescale/dt) dt = timescale/float(ntimes) do n = 1,ntimes c c-----Cloud tracer change c do l = 1,nlay do k = 1,nlay if (k.eq.1) then dcflx = -cflx(k)*cctrac(k,l) if (cflx(k).lt.0.) dcflx = -cflx(k)*cctrac(k+1,l) else fm = cflx(k-1)*cctrac(k-1,l) if (cflx(k-1).lt.0.) fm = cflx(k-1)*cctrac(k,l) fp = cflx(k)*cctrac(k,l) if (cflx(k).lt.0.) then fp = 0. if (k.lt.nlay) fp = cflx(k)*cctrac(k+1,l) endif dcflx = fm - fp endif ccnew(k,l) = cctrac(k,l) + (dt/dz(k))*(1./rho(k))* & (dcflx + entc(k)*catrac(k,l)*crat - detc(k)*cctrac(k,l)) if (k.eq.1) then daflx = -aflx(k)*catrac(k,l) if (aflx(k).lt.0.) daflx = -aflx(k)*catrac(k+1,l) else fm = aflx(k-1)*catrac(k-1,l) if (aflx(k-1).lt.0.) fm = aflx(k-1)*catrac(k,l) fp = aflx(k)*catrac(k,l) if (aflx(k).lt.0.) then fp = 0. if (k.lt.nlay) fp = aflx(k)*catrac(k+1,l) endif daflx = fm - fp endif canew(k,l) = catrac(k,l) + (dt/dz(k))*(1./rho(k))* & (daflx - enta(k)*catrac(k,l) + deta(k)*cctrac(k,l)/crat) enddo do k = 1,nlay cctrac(k,l) = ccnew(k,l) catrac(k,l) = canew(k,l) enddo enddo c c-----Ambient tracer change c do l = 1,nlay do k = 1,nlay if (k.eq.1) then daflx = -aflx(k)*aatrac(k,l) if (aflx(k).lt.0.) daflx = -aflx(k)*aatrac(k+1,l) else fm = aflx(k-1)*aatrac(k-1,l) if (aflx(k-1).lt.0.) fm = aflx(k-1)*aatrac(k,l) fp = aflx(k)*aatrac(k,l) if (aflx(k).lt.0.) then fp = 0. if (k.lt.nlay) fp = aflx(k)*aatrac(k+1,l) endif daflx = fm - fp endif aanew(k,l) = aatrac(k,l) + (dt/dz(k))*(1./rho(k))* & (daflx - enta(k)*aatrac(k,l) + deta(k)*actrac(k,l)/crat) if (k.eq.1) then dcflx = -cflx(k)*actrac(k,l) if (cflx(k).lt.0.) dcflx = -cflx(k)*actrac(k+1,l) else fm = cflx(k-1)*actrac(k-1,l) if (cflx(k-1).lt.0.) fm = cflx(k-1)*actrac(k,l) fp = cflx(k)*actrac(k,l) if (cflx(k).lt.0.) then fp = 0. if (k.lt.nlay) fp = cflx(k)*actrac(k+1,l) endif dcflx = fm - fp endif acnew(k,l) = actrac(k,l) + (dt/dz(k))*(1./rho(k))* & (dcflx + entc(k)*aatrac(k,l)*crat - detc(k)*actrac(k,l)) enddo do k = 1,nlay aatrac(k,l) = aanew(k,l) actrac(k,l) = acnew(k,l) enddo enddo enddo ! End time loop c c-----Load CiG transport matrix array c do ks = 1,nlay !Starting layer do ke = 1,nlay !Ending layer cigmas(1,ke,ks,i,j) = cctrac(ke,ks) + actrac(ke,ks)/crat cigmas(2,ke,ks,i,j) = aatrac(ke,ks) + catrac(ke,ks)*crat enddo enddo c c-----Check that transport matrices conserve mass c do ks = 1,nlay tmassi = rho(ks)*dz(ks) tmassf = 0. do ke = 1,nlay if (cigmas(1,ke,ks,i,j).lt.0. .or. & cigmas(2,ke,ks,i,j).lt.0.) then write(iout,'(//,a)') 'ERROR in CONVTRANS:' write(iout,'(2a)') 'Negative tracer after', & ' convective transport' write(iout,*) 'igrid, i, j, ks, ke = ',igrid,i,j,ks,ke write(iout,*) 'Spec Name In-cloud Ambient' do kk = 1,nlay write(iout,'(2i3,2e10.3)') & ks,kk,cigmas(1,kk,ks,i,j),cigmas(2,kk,ks,i,j) enddo call camxerr() endif tmassf = tmassf + rho(ke)*dz(ke)* & (cigmas(1,ke,ks,i,j)*cigfrc(i,j) + & cigmas(2,ke,ks,i,j)*(1.-cigfrc(i,j))) enddo terr = abs(tmassf - tmassi)/tmassi if (terr.gt.0.001) then write(iout,'(//,a)')'ERROR in CONVTRANS:' write(iout,*)'CiG transport matrix not conserving mass' write(iout,*)'Grid,I,J,K : ',igrid,i,j,ks write(iout,*)'Starting/ending mass: ',tmassi,tmassf call camxerr() endif enddo c 10 continue 20 continue return end