MODULE CloudModule USE CloudTab IMPLICIT none PRIVATE PUBLIC :: initCloud, setCloudOptProfile, setCloudOptLayer, destroyCloud INTEGER,PARAMETER,PUBLIC :: iCLW=1,iIce=2,iRain=3,iSnow=4,iGraupel=5 INTEGER,PARAMETER,PUBLIC :: nHydNonPcp=2 INTEGER,PARAMETER :: iKabs=1,iKscat=2,iAsym=3 INTEGER,PARAMETER :: iCTPice=0,iThkIce=1,iTempIce=2,iAmtIce=3,iDeffIce=4 INTEGER,PARAMETER :: iCTPliq=0,iThkLiq=1,iAmtLiq=2,iDeffLiq=3 INTEGER,PARAMETER :: mxLev=120 REAL,DIMENSION(:),ALLOCATABLE :: kabs,kscat,asymi,itscat INTEGER,DIMENSION(:),ALLOCATABLE :: hydIndx,sizIndx,amtIndx,topIndx,thkIndx INTEGER,DIMENSION(:,:),ALLOCATABLE :: indxArr INTEGER :: nSpcPts,nHydLoc REAL,DIMENSION(:),ALLOCATABLE :: intAmt CONTAINS SUBROUTINE initCloud(ctab_file,dataType_in,spectPtsGrd,iCldLiq,iCldIce, & cldTyp,nHyd) USE constants, ONLY: MISSING_CHAR CHARACTER(LEN=*),INTENT(IN) :: ctab_file CHARACTER(LEN=*),DIMENSION(:),INTENT(IN) :: dataType_in CHARACTER(LEN=1),PARAMETER :: spc='_' CHARACTER(LEN=20),DIMENSION(:),ALLOCATABLE :: dataType INTEGER, INTENT(IN),OPTIONAL :: iCldLiq,iCldIce INTEGER :: dSize INTEGER :: i,j,dim1,ict INTEGER,PARAMETER :: nHydDflt=5 INTEGER,PARAMETER :: nOptical=3 INTEGER,DIMENSION(:),ALLOCATABLE :: indxArr1d CHARACTER(LEN=5),DIMENSION(nOptical),PARAMETER :: optIDs=(/'kabs ','kscat','gasym'/) REAL,DIMENSION(:),POINTER :: spectPtsGrd CHARACTER(LEN=1) :: cdum CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: cldTyp INTEGER,INTENT(IN),OPTIONAL :: nHyd if(present(nHyd))then if(nHyd.gt.nHydDflt)then print*,'err[CloudModule::initCloud]:' print*,'Input nHyd exceeds the default (maximum) nHyd value: ' print*,nHyd,nHydDflt call errorHalt(1) endif nHydLoc=nHyd else nHydLoc=nHydDflt endif if(present(cldTyp))then if(trim(adjustl(cldTyp)).eq.'layer')nHydLoc=nHydNonPcp endif ict=0 do i=1,size(dataType_in) cdum = dataType_in(i) if(cdum(1:1).ne.MISSING_CHAR)ict=ict+1 enddo if(ict.lt.nHydLoc)then print*,'err[CloudModule::initCloud]:' print*,'nHydLoc parameter greater than number of requested ' print*,' hydrometeors in dataType character vector.' print*,'nHydLoc = ',nHydLoc print*,'# of requested hydrometeors: ',ict call errorHalt(1) endif allocate(dataType(nHydLoc*nOptical)) dSize = size(dataType) allocate(indxArr1d(dSize)) indxArr1d = (/(i,i=1,dSize)/) dim1 = dsize/nOptical allocate(indxArr(dim1,nOptical)) indxArr = reshape(indxArr1d,(/dim1,nOptical/),order=(/2,1/)) if(present(cldTyp))then if(trim(adjustl(cldTyp)).eq.'layer')then allocate(hydIndx(nHydLoc),sizIndx(nHydLoc),amtIndx(nHydLoc)) allocate(thkIndx(nHydLoc),topIndx(nHydLoc)) hydIndx = (/iCldLiq,iCldIce/) topIndx = (/iCTPliq,iCTPice/) thkIndx = (/iThkLiq,iThkIce/) amtIndx = (/iAmtLiq,iAmtIce/) sizIndx = (/iDeffLiq,iDeffIce/) else if(trim(adjustl(cldTyp)).ne.'profile')then print*,'err[CldModule::initCloud]: Invalid cldTyp parameter.' call errorHalt(1) endif endif endif ! indxArr logic assumes all three optical properties are requested ! for each hydrometeor type; thus, 'nOptical' evenly divides 'dSize' !--Prepare dataType by combining hydrometeor and optical property ID strings: ict=1 do i=1,nHydLoc do j=1,nOptical dataType(ict) = trim(dataType_in(i)) // spc // trim(optIDs(j)) ict=ict+1 enddo enddo CALL cloudTabInit(ctab_file,dataType,outBndsExit_in=.false.,nSpcPts=nSpcPts, & spectPtsGrd=spectPtsGrd) allocate(kabs(nSpcPts),kscat(nSpcPts),asymi(nSpcPts),itscat(nSpcPts)) allocate(intAmt(nHydLoc)) deallocate(indxArr1d,dataType) RETURN END SUBROUTINE initCloud ! =================================== SUBROUTINE setCloudOptProfile(xG,tavl,hydr,sizP,nSurf,pref,psfc,ih2o,tabs,tscat,asym) USE constants, ONLY : rd,g0,ratioMwt !---Input variables: INTEGER, INTENT(IN) :: nSurf,ih2o REAL, DIMENSION(:), INTENT(IN) :: xG,tavl,pref REAL, DIMENSION(:,:),INTENT(IN) :: hydr,sizP REAL, INTENT(IN) :: psfc !---Output variables: REAL, DIMENSION(:,:),INTENT(INOUT) :: tabs,tscat,asym !---Local variables: REAL, DIMENSION(mxLev,nHydLoc) :: uAmt,mcon,rdum1,rdum2 REAL, DIMENSION(nHydLoc) :: mconSfc INTEGER :: iLyr,iLev,iHyd,iSpc,n2 REAL :: hydrSum,scal REAL,PARAMETER :: amtScl=1e5 tabs = 0. tscat = 0. asym = 0. n2 = nSurf-1 !---Compute mass concentration of each hydrometeor class on the fixed ! pressure levels: do iLev=1,nSurf hydrSum = sum(hydr(iLev,:)) do iHyd=1,nHydLoc mcon(iLev,iHyd) = hydr(iLev,iHyd) / (1 + xG(ih2o+(iLev)-1) + hydrSum) enddo enddo !---Compute hydrometeor layer amount: scal = amtScl/g0 ! convert p(mb) to g/m^2 ; (1mb = 100 kg/(m s^2)) do iLyr=1,n2-1 do iHyd=1,nHydLoc if(max(mcon(iLyr,iHyd),mcon(iLyr+1,iHyd)).gt.epsilon(1.))then call lpsum(pref(iLyr),pref(iLyr+1),mcon(iLyr,iHyd), & mcon(iLyr+1,iHyd),scal, & uAmt(iLyr,iHyd),rdum1(iLyr,iHyd),rdum2(iLyr,iHyd)) else uAmt(iLyr,iHyd)=0. endif enddo enddo !---Surface layer: do iHyd=1,nHydLoc if(max(mcon(nSurf,iHyd),mcon(nSurf-1,iHyd)).gt.epsilon(1.))then call lint(mcon(:,iHyd),pref,nSurf,psfc,mconSfc(iHyd)) call lpsum(pref(n2),psfc,mcon(n2,iHyd),mconSfc(iHyd),scal,uAmt(n2,iHyd), & rdum1(n2,iHyd),rdum2(n2,iHyd)) else uAmt(n2,iHyd)=0. endif enddo itscat = 0. do iLyr=1,nSurf-1 !---Get tabulated optical properties, convert to optical depths, and sum ! over hydrometeor classes: do iHyd=1,nHydLoc if(uAmt(iLyr,iHyd).lt.epsilon(1.)) cycle call getCloudTab(indxArr(iHyd,iKabs),sizP(iLyr,iHyd),0.,tavl(iLyr),kabs(:)) call getCloudTab(indxArr(iHyd,iKscat),sizP(iLyr,iHyd),0.,tavl(iLyr),kscat(:)) call getCloudTab(indxArr(iHyd,iAsym),sizP(iLyr,iHyd),0.,tavl(iLyr),asymi(:)) itscat = kscat(:) * uAmt(iLyr,iHyd) tabs(iLyr,:) = tabs(iLyr,:) + kabs(:) *uAmt(iLyr,iHyd) tscat(iLyr,:) = tscat(iLyr,:) + itscat(:) asym(iLyr,:) = asym(iLyr,:) + itscat(:)*asymi(:) enddo !---The following loop could be replaced with a WHERE construct, although the ! current version of the SGI compiler does not have the proper libraries ! to support WHERE in DEBUG mode. do iSpc=1,nSpcPts if(tscat(iLyr,iSpc).gt.epsilon(1.))then asym(iLyr,iSpc) = asym(iLyr,iSpc) / tscat(iLyr,iSpc) else asym(iLyr,iSpc) = 0. endif enddo enddo RETURN END SUBROUTINE setCloudOptProfile ! =================================== SUBROUTINE setCloudOptLayer(xG,tavl,nSurf,pref,psfc,iCldLiq,iCldIce,tabs,tscat,asym) !---Input variables: INTEGER, INTENT(IN) :: nSurf,iCldLiq,iCldIce REAL, DIMENSION(:), INTENT(IN) :: xG,tavl,pref REAL, INTENT(IN) :: psfc !---Output variables: REAL, DIMENSION(:,:),INTENT(INOUT) :: tabs,tscat,asym !---Local variables: REAL, DIMENSION(mxLev,nHydLoc) :: cFrac INTEGER,DIMENSION(mxLev) :: ltyp INTEGER :: iLyr,iLev,iHyd,iSpc,n2 INTEGER :: topID,thkID,amtID INTEGER :: j,j1,j2 REAL :: hydrSum,scal,sizPs,cldTop,cldBot,cldThk REAL,PARAMETER :: amtScl=1e5,clwScl=1000. REAL,PARAMETER,DIMENSION(nHydNonPcp) :: sclPath=(/1000.,1./) real :: ttemp tabs = 0. tscat = 0. asym = 0. n2 = nSurf-1 ltyp(:) = 0 cFrac(:,:) = 0. intAmt(:) = 0. !-- Input: several-parameter cloud model ! - Determine fraction of cloud falling within each layer (for each hydro type). do iHyd=1,nHydLoc topID = hydIndx(iHyd)+topIndx(iHyd) thkID = hydIndx(iHyd)+thkIndx(iHyd) amtID = hydIndx(iHyd)+amtIndx(iHyd) intAmt(iHyd) = xG(hydIndx(iHyd)+amtIndx(iHyd)) * sclPath(iHyd) if (xG(amtID).gt.epsilon(1.))then cldTop = xG(topID) cldBot = xG(topID)+xG(thkID) cldThk = xG(thkID) do j=1,n2 if(cldTop.lt.pref(j)) exit enddo j1=MIN(j,n2) !-----Find index of first standard level at or below cloud base. do j=j1-1,n2 if(cldBot.le.pref(j)) exit enddo j2=MIN(j,n2) !-----Calculate cloud amt in each layer cFrac(j1-1,iHyd)=(pref(j1)-cldTop)/cldThk ltyp(j1-1)=iHyd do j=j1,j2-2 cFrac(j,iHyd)=(pref(j+1)-pref(j))/cldthk ltyp(j)=iHyd enddo cFrac(j2-1,iHyd)=(cldbot-amax1(pref(j2-1),cldtop))/cldthk !bottom layer ltyp(j2-1)=iHyd endif enddo ! - Given constant Deff, and T from profile, get opt. depth profiles do iLyr=1,nSurf-1 do iHyd=1,nHydLoc if (cFrac(iLyr,iHyd).lt.epsilon(1.)) cycle sizPs=xG(hydIndx(iHyd)+sizIndx(iHyd)) call getCloudTab(indxArr(iHyd,iKabs),sizPs,0.,tavl(iLyr),kabs(:)) call getCloudTab(indxArr(iHyd,iKscat),sizPs,0.,tavl(iLyr),kscat(:)) call getCloudTab(indxArr(iHyd,iAsym),sizPs,0.,tavl(iLyr),asymi(:)) itscat = kscat(:) * cFrac(iLyr,iHyd) * intAmt(iHyd) tabs(iLyr,:) = tabs(iLyr,:) + kabs(:) * cFrac(iLyr,iHyd) * intAmt(iHyd) tscat(iLyr,:) = tscat(iLyr,:) + itscat(:) asym(iLyr,:) = asym(iLyr,:) + itscat(:)*asymi(:) enddo !---The following loop could be replaced with a WHERE construct, although the ! current version of the SGI compiler does not have the proper libraries ! to support WHERE in DEBUG mode. do iSpc=1,nSpcPts if(tscat(iLyr,iSpc).gt.epsilon(1.))then asym(iLyr,iSpc) = asym(iLyr,iSpc) / tscat(iLyr,iSpc) else asym(iLyr,iSpc) = 0. endif enddo enddo ! ltyp may be passed to ossdrv and input to scattRT, ! to avoid repeatedly determining ltyp there. RETURN END SUBROUTINE setCloudOptLayer ! =================================== SUBROUTINE destroyCloud(nindx) integer, intent(in) :: nindx call cloudTabDestroy(nindx) if (allocated(kabs)) deallocate(kabs) if (allocated(kscat)) deallocate(kscat) if (allocated(asymi)) deallocate(asymi) if (allocated(itscat)) deallocate(itscat) if (allocated(indxArr))deallocate(indxArr) if (allocated(intAmt)) deallocate(intAmt) if (allocated(hydIndx))deallocate(hydIndx) if (allocated(sizIndx))deallocate(sizIndx) if (allocated(amtIndx))deallocate(amtIndx) if (allocated(thkIndx))deallocate(thkIndx) if (allocated(topIndx))deallocate(topIndx) RETURN END SUBROUTINE destroyCloud ! The following interpolation and layer averaging routines were copied here from ! oss_mw as a short-term fix. They should either become public subroutines of ! oss_mw, or be placed in a more general module. !---------------------------------------------------------------------------- ! PURPOSE: Computes average layer quantities (or integrated amount) ! using a linear dependence on P. !---------------------------------------------------------------------------- SUBROUTINE lpsum(pu,pl,xu,xl,scal,xint,dxu,dxl) !---Input variables REAL, INTENT(IN) :: pu,pl,xu,xl,scal !---Output variables REAL, INTENT(OUT) :: xint,dxu,dxl xint = 0.5*(pl-pu)*scal dxu = xint dxl = xint xint = xint*(xu+xl) RETURN END SUBROUTINE lpsum !---------------------------------------------------------------------------- ! PURPOSE: Computes average layer quantities (or integrated amount) ! using a log-x dependence on log-p. !---------------------------------------------------------------------------- SUBROUTINE lpsum_log(pu,pl,xu,xl,scal,xint,dxu,dxl) REAL, PARAMETER :: epsiln=1.e-05 !---Input variables REAL, INTENT(IN) :: pu,pl,xu,xl,scal !---Output variables REAL, INTENT(OUT) :: xint,dxu,dxl !---Local variables REAL :: hp,x0,zeta,alza,alpha hp = alog(pl/pu) x0 = pl*xl*hp zeta = pu*xu/(pl*xl) IF(ABS(zeta-1.).GT.epsiln)THEN alza = alog(zeta) xint = x0*(zeta-1.)/alza alpha = zeta/(zeta-1.)-1./alza ELSE xint = x0*2./(3.-zeta) alpha = zeta/(3.-zeta) END IF xint = xint*scal dxu = xint*alpha/xu dxl = xint*(1.-alpha)/xl RETURN END SUBROUTINE lpsum_log SUBROUTINE lint(x,p,N0,p0,x0) INTEGER :: n0 REAL :: p0 REAL, DIMENSION(:), INTENT(IN) :: x,p !---Output variable REAL :: x0 !---Local variable REAL :: xx xx = (x(N0)-x(N0-1))/(p(N0)-p(N0-1)) x0 = x(N0-1)+(p0-p(N0-1))*xx RETURN END SUBROUTINE lint SUBROUTINE lint_log(x,p,N0,p0,x0) INTEGER :: n0 REAL :: p0 REAL, DIMENSION(:), INTENT(IN) :: x,p !---Output variable REAL :: x0 !---Local variable REAL :: xx xx = alog(x(N0)/x(N0-1))/alog(p(N0)/p(N0-1)) x0 = x(N0-1)*(p0/p(N0-1))**xx RETURN END SUBROUTINE lint_log END MODULE CloudModule