LCOV - code coverage report
Current view: | top level - physics/pumas - pumas_stochastic_collect_tau.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 285 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module pumas_stochastic_collect_tau
           2             : ! From Morrison (Lebo, originally TAU bin code)
           3             : ! Gettelman and Chen 2018
           4             : !the subroutines take in air density, air temperature, and the bin mass boundaries, and
           5             : !output the mass and number mixing ratio tendencies in each bin directly.
           6             : !this is then wrapped for CAM.
           7             :
           8             : ! note, this is now coded locally. Want the CAM interface to be over i,k I think.
           9             :
          10             : #ifndef HAVE_GAMMA_INTRINSICS
          11             : use shr_spfn_mod, only: gamma => shr_spfn_gamma
          12             : #endif
          13             :
          14             : use shr_kind_mod,      only: r8=>shr_kind_r8
          15             : use cam_history,       only: addfld
          16             : use micro_pumas_utils, only: pi, rhow, qsmall, VLENS
          17             : use cam_logfile,       only: iulog
          18             :
          19             : implicit none
          20             : private
          21             : save
          22             :
          23             : ! Subroutines
          24             : public :: pumas_stochastic_kernel_init, pumas_stochastic_collect_tau_tend
          25             :
          26             : !In the module top, declare the following so that these can be used throughout the module:
          27             :
          28             : integer, parameter, public  :: ncd = 35
          29             : integer, parameter, public  :: ncdp = ncd + 1
          30             : integer, parameter, public  :: ncdl = ncd
          31             : integer, parameter, public  :: ncdpl = ncdl+1
          32             :
          33             : ! for Zach's collision-coalescence code
          34             :
          35             : real(r8), private :: knn(ncd,ncd)
          36             :
          37             : real(r8), private :: mmean(ncd), diammean(ncd)       ! kg & m at bin mid-points
          38             : real(r8), private :: medge(ncdp), diamedge(ncdp)     ! kg & m at bin edges
          39             : integer, private  :: cutoff_id                       ! cutoff between cloud water and rain drop, D = 40 microns
          40             :
          41             : ! Assume 6 microns for each...
          42             : real(r8), parameter :: m1 = 4._r8/3._r8*pi*rhow*(6.e-6_r8)**3
          43             :
          44             : !$acc declare create(knn,cutoff_id,mmean,diammean,medge,diamedge)
          45             :
          46             : !===============================================================================
          47             : contains
          48             : !===============================================================================
          49             :
          50           0 : subroutine calc_bins
          51             :
          52             :   real(r8) :: DIAM(ncdp)
          53             :   real(r8) :: X(ncdp)
          54             :   real(r8) :: radsl(ncdp)
          55             :   real(r8) :: radl(ncd)
          56             :   integer  :: L, lcl
          57             :   real(r8) :: kkfac
          58             :
          59             : !Then before doing any calculations you'll need to calculate the bin mass grid
          60             : ! (note this code could be cleaned up, I'm just taking it as it's used in our bin scheme).
          61             : ! This only needs to be done once, since we'll use the same bin mass grid for all calculations.
          62             :
          63             : ! use mass doubling bins from Graham Feingold (note cgs units)
          64             :
          65           0 :   DIAM(1)=1.5625*2.E-04_r8                ! cm
          66           0 :   X(1)=PI/6._r8*DIAM(1)**3*rhow/1000._r8  ! rhow kg/m3 --> g/cm3
          67           0 :   radsl(1) = X(1)                         ! grams
          68             :
          69           0 :   DO l=2,ncdp
          70           0 :      X(l)=2._r8*X(l-1)
          71           0 :      DIAM(l)=(6._r8/pi*X(l)*1000._r8/rhow)**(1._r8/3._r8)  ! cm
          72           0 :      radsl(l)=X(l)
          73             :   ENDDO
          74             :
          75             : ! now get bin mid-points
          76             :
          77           0 :   do l=1,ncd
          78           0 :      radl(l)=(radsl(l)+radsl(l+1))/2._r8         ! grams
          79           0 :      diammean(l) = (6._r8/pi*radl(l)*1000._r8/rhow)**(1._r8/3._r8) ! cm
          80             :   end do
          81             :
          82             : ! set bin grid for method of moments
          83             :
          84             :   ! for method of moments
          85             :
          86           0 :   do lcl = 1,ncd+1
          87           0 :      medge(lcl) = radsl(lcl)               ! grams
          88           0 :      diamedge(lcl) = DIAM(lcl)             ! cm
          89             :   enddo
          90             :
          91           0 :   do lcl = 1,ncd
          92           0 :      mmean(lcl) = radl(lcl)
          93           0 :      diammean(lcl) = diammean(lcl)
          94             :   enddo
          95             :
          96           0 :   do lcl = ncdp,1,-1
          97           0 :      if( diamedge(lcl).ge.40.e-4_r8 ) cutoff_id = lcl
          98             :   end do
          99             :
         100           0 : end subroutine calc_bins
         101             :
         102           0 : subroutine pumas_stochastic_kernel_init(kernel_filename)
         103             :
         104             :   use cam_history_support, only: add_hist_coord
         105             :
         106             :   character(len=*), intent(in) :: kernel_filename  ! Full pathname to kernel file
         107             :
         108             :   integer :: iunit ! unit number of opened file for collection kernel code from a lookup table.
         109             :
         110             :   integer :: idd, jdd
         111             :   real(r8) :: kkfac
         112             :
         113           0 :   call calc_bins
         114             :
         115             : ! Read in the collection kernel code from a lookup table. Again, this only needs to be done once.
         116             : ! use kernel from Zach (who got it from Jerry)
         117             :
         118           0 :   KNN(:,:)=0._r8 ! initialize values
         119           0 :   kkfac=1.5_r8   ! from Zach
         120             :
         121           0 :   open(newunit=iunit,file=kernel_filename,status='old')
         122             :
         123           0 :   do idd=1,ncd
         124           0 :      do jdd=1,idd
         125           0 :         READ(iunit,941) KNN(IDD,JDD)
         126             : 941     FORMAT(2X,E12.5)
         127             :
         128           0 :         KNN(IDD,JDD)=(mmean(IDD)*kkfac+mmean(JDD)*kkfac)*KNN(IDD,JDD)
         129           0 :         if (knn(idd,jdd) < 0._r8) knn(idd,jdd)=0._r8
         130             :      end do
         131             :   end do
         132             :
         133             :   !$acc update device(knn,cutoff_id,mmean,diammean,medge,diamedge)
         134             :
         135           0 : end subroutine pumas_stochastic_kernel_init
         136             :
         137             : !main driver routine
         138             : !needs to pull in i,k fields (so might need dimensions here too)
         139             :
         140           0 : subroutine pumas_stochastic_collect_tau_tend(deltatin,t,rho,qc,qr,qcin,     &
         141           0 :                         ncin,qrin,nrin,lcldm,precip_frac,mu_c,lambda_c,     &
         142           0 :                         n0r,lambda_r,qcin_new,ncin_new,qrin_new,nrin_new,   &
         143           0 :                         qctend,nctend,qrtend,nrtend,qctend_TAU,nctend_TAU,  &
         144           0 :                         qrtend_TAU,nrtend_TAU,scale_qc,scale_nc,scale_qr,   &
         145           0 :                         scale_nr,amk_c,ank_c,amk_r,ank_r,amk,ank,amk_out,   &
         146           0 :                         ank_out,gmnnn_lmnnn_TAU,mgncol,nlev)
         147             :
         148             :   !inputs: mgncol,nlev,t,rho,qcin,ncin,qrin,nrin
         149             :   !outputs: qctend,nctend,qrtend,nrtend
         150             :   !not sure if we want to output bins (extra dimension). Good for testing?
         151             :
         152             :   integer, intent(in) :: mgncol,nlev
         153             :   real(r8), intent(in) :: deltatin
         154             :   real(r8), intent(in) :: t(mgncol,nlev)
         155             :   real(r8), intent(in) :: rho(mgncol,nlev)
         156             :   real(r8), intent(in) :: qc(mgncol,nlev)
         157             :   real(r8), intent(in) :: qr(mgncol,nlev)
         158             :   real(r8), intent(in) :: qcin(mgncol,nlev)
         159             :   real(r8), intent(in) :: ncin(mgncol,nlev)
         160             :   real(r8), intent(in) :: qrin(mgncol,nlev)
         161             :   real(r8), intent(in) :: nrin(mgncol,nlev)
         162             :   real(r8), intent(in) :: lcldm(mgncol,nlev)
         163             :   real(r8), intent(in) :: precip_frac(mgncol,nlev)
         164             :   real(r8), intent(inout) :: qctend(mgncol,nlev)
         165             :   real(r8), intent(inout) :: nctend(mgncol,nlev)
         166             :   real(r8), intent(inout) :: qrtend(mgncol,nlev)
         167             :   real(r8), intent(inout) :: nrtend(mgncol,nlev)
         168             :   real(r8), intent(out) :: qctend_TAU(mgncol,nlev)
         169             :   real(r8), intent(out) :: nctend_TAU(mgncol,nlev)
         170             :   real(r8), intent(out) :: qrtend_TAU(mgncol,nlev)
         171             :   real(r8), intent(out) :: nrtend_TAU(mgncol,nlev)
         172             :
         173             :   real(r8), intent(out) :: scale_qc(mgncol,nlev)
         174             :   real(r8), intent(out) :: scale_nc(mgncol,nlev)
         175             :   real(r8), intent(out) :: scale_qr(mgncol,nlev)
         176             :   real(r8), intent(out) :: scale_nr(mgncol,nlev)
         177             :
         178             :   real(r8), intent(out) :: amk_c(mgncol,nlev,ncd)
         179             :   real(r8), intent(out) :: ank_c(mgncol,nlev,ncd)
         180             :   real(r8), intent(out) :: amk_r(mgncol,nlev,ncd)
         181             :   real(r8), intent(out) :: ank_r(mgncol,nlev,ncd)
         182             :   real(r8), intent(out) :: amk(mgncol,nlev,ncd)
         183             :   real(r8), intent(out) :: ank(mgncol,nlev,ncd)
         184             :   real(r8), intent(out) :: amk_out(mgncol,nlev,ncd)
         185             :   real(r8), intent(out) :: ank_out(mgncol,nlev,ncd)
         186             :
         187             :   real(r8), intent(out) :: mu_c(mgncol,nlev)
         188             :   real(r8), intent(out) :: lambda_c(mgncol,nlev)
         189             :   real(r8), intent(out) :: lambda_r(mgncol,nlev)
         190             :   real(r8), intent(out) :: n0r(mgncol,nlev)
         191             :
         192             :   real(r8), intent(out) :: qcin_new(mgncol,nlev)
         193             :   real(r8), intent(out) :: ncin_new(mgncol,nlev)
         194             :   real(r8), intent(out) :: qrin_new(mgncol,nlev)
         195             :   real(r8), intent(out) :: nrin_new(mgncol,nlev)
         196             :   real(r8), intent(out) :: gmnnn_lmnnn_TAU(mgncol,nlev)
         197             :
         198             :   ! Local variables
         199             :
         200             :   integer :: i,k,n,lcl
         201           0 :   integer :: cutoff_amk(mgncol,nlev),cutoff(mgncol,nlev)
         202             :
         203             :   real(r8) :: all_gmnnn,all_lmnnn
         204             :   real(r8) :: qscl
         205             :
         206           0 :   real(r8) :: qcin_old(mgncol,nlev)
         207           0 :   real(r8) :: ncin_old(mgncol,nlev)
         208           0 :   real(r8) :: qrin_old(mgncol,nlev)
         209           0 :   real(r8) :: nrin_old(mgncol,nlev)
         210             :
         211           0 :   real(r8) :: amk0(mgncol,nlev,ncd)
         212           0 :   real(r8) :: ank0(mgncol,nlev,ncd)
         213           0 :   real(r8) :: gnnnn(mgncol,nlev,ncd)
         214           0 :   real(r8) :: gmnnn(mgncol,nlev,ncd)
         215           0 :   real(r8) :: lnnnn(mgncol,nlev,ncd)
         216           0 :   real(r8) :: lmnnn(mgncol,nlev,ncd)
         217           0 :   real(r8) :: gnnnn0(mgncol,nlev,ncd)
         218           0 :   real(r8) :: gmnnn0(mgncol,nlev,ncd)
         219           0 :   real(r8) :: lnnnn0(mgncol,nlev,ncd)
         220           0 :   real(r8) :: lmnnn0(mgncol,nlev,ncd)
         221             :
         222             :   integer, parameter :: sub_step = 60
         223             :
         224             :   !$acc data create (cutoff_amk,cutoff,qcin_old,ncin_old,qrin_old, &
         225             :   !$acc              nrin_old,amk0,ank0,gnnnn,gmnnn,lnnnn,lmnnn,   &
         226             :   !$acc              gnnnn0,gmnnn0,lnnnn0,lmnnn0)
         227             :
         228             :   !$acc parallel vector_length(VLENS) default(present)
         229             :   !$acc loop gang vector collapse(2)
         230           0 :   do k=1,nlev
         231           0 :      do i=1,mgncol
         232           0 :         cutoff(i,k) = cutoff_id - 1
         233             :      end do
         234             :   end do
         235             :   !$acc end parallel
         236             :
         237             :   ! First make bins from cam size distribution (bins are diagnostic)
         238             :
         239             :   call cam_bin_distribute(qc,qr,qcin,ncin,qrin,nrin,mu_c,lambda_c, &
         240             :                           lambda_r,n0r,lcldm,precip_frac,scale_qc, &
         241             :                           scale_nc,scale_qr,scale_nr,amk_c,ank_c,  &
         242           0 :                           amk_r,ank_r,amk,ank,cutoff_amk,mgncol,nlev)
         243             :
         244             :   !$acc parallel vector_length(VLENS) default(present)
         245             :   !$acc loop gang vector collapse(2)
         246           0 :   do k=1,nlev
         247           0 :      do i=1,mgncol
         248           0 :         if ( cutoff_amk(i,k) > 0 ) then
         249           0 :            cutoff(i,k) = cutoff_amk(i,k)
         250             :         end if
         251             :      end do
         252             :   end do
         253             :   !$acc end parallel
         254             :
         255             : !Then call the subroutines that actually do the calculations. The inputs/ouputs are described in comments below.
         256             :
         257             : !This part of the code needs to be called each time for each process rate calculation
         258             : ! (i.e., for each sampled cloud/rain gamma distribution):
         259             :
         260             : ! note: variables passed to compute_column_params are all inputs,
         261             : ! outputs from this subroutine are stored as global variables
         262             :
         263             : ! inputs: t --> input air temperature (K)
         264             : !         rho --> input air density (kg/m^3)
         265             : !         medge --> bin mass boundary (g)
         266             : !         amk --> array of bin mass mixing ratio, i.e., the input drop mass distribution (kg/kg)
         267             : !         ank --> array of bin number mixing ratio, i.e., the input drop number distribution (kg^-1)
         268             :
         269             : ! inputs: medge --> bin mass boundary (g), same as above
         270             :
         271             : ! outputs: gnnnn --> bin number mass mixing tendency gain, array in bins (#/cm^3/s)
         272             : !          lnnnn --> bin number mass mixing tendency loss, array in bins (#/cm^3/s)
         273             : !          gmnnn --> bin mass mixing ratio tendency gain, array in bins (g/cm^3/s)
         274             : !          lmnnn --> bin mass mixing ratio tendency loss, array in bins (g/cm^3/s)
         275             :
         276             :
         277             : ! Call Kernel
         278             :
         279             :   !$acc parallel vector_length(VLENS) default(present)
         280             :   !$acc loop gang vector collapse(2)
         281           0 :   do k=1,nlev
         282           0 :      do i=1,mgncol
         283           0 :         qcin_new(i,k) = 0._r8
         284           0 :         ncin_new(i,k) = 0._r8
         285           0 :         qrin_new(i,k) = 0._r8
         286           0 :         nrin_new(i,k) = 0._r8
         287             :
         288           0 :         qcin_old(i,k) = 0._r8
         289           0 :         ncin_old(i,k) = 0._r8
         290           0 :         qrin_old(i,k) = 0._r8
         291           0 :         nrin_old(i,k) = 0._r8
         292             :
         293           0 :         qctend_TAU(i,k) = 0._r8
         294           0 :         nctend_TAU(i,k) = 0._r8
         295           0 :         qrtend_TAU(i,k) = 0._r8
         296           0 :         nrtend_TAU(i,k) = 0._r8
         297             :      end do
         298             :   end do
         299             :   !$acc end parallel
         300             :
         301             :   !$acc parallel vector_length(VLENS) default(present)
         302             :   !$acc loop gang vector collapse(3)
         303           0 :   do lcl=1,ncd
         304           0 :      do k=1,nlev
         305           0 :         do i=1,mgncol
         306           0 :            gnnnn(i,k,lcl) = 0._r8
         307           0 :            gmnnn(i,k,lcl) = 0._r8
         308           0 :            lnnnn(i,k,lcl) = 0._r8
         309           0 :            lmnnn(i,k,lcl) = 0._r8
         310             :         end do
         311             :      end do
         312             :   end do
         313             :   !$acc end parallel
         314             :
         315             : ! update qc, nc, qr, nr
         316             :
         317             :   !$acc parallel vector_length(VLENS) default(present)
         318             :   !$acc loop gang vector collapse(2)
         319           0 :   do k=1,nlev
         320           0 :      do i=1,mgncol
         321             :         !$acc loop seq
         322           0 :         do lcl=1,ncd
         323           0 :            amk0(i,k,lcl) = amk(i,k,lcl)
         324           0 :            ank0(i,k,lcl) = ank(i,k,lcl)
         325             :         end do
         326             :         ! substep bin code
         327             :         !$acc loop seq
         328           0 :         do n=1,sub_step
         329           0 :            call compute_coll_params(rho(i,k),medge,amk0(i,k,1:ncd),ank0(i,k,1:ncd),gnnnn0(i,k,1:ncd),gmnnn0(i,k,1:ncd),lnnnn0(i,k,1:ncd),lmnnn0(i,k,1:ncd))
         330             :
         331           0 :            all_gmnnn=0._r8
         332           0 :            all_lmnnn=0._r8
         333             :            !scaling gmnnn, lmnnn
         334             :            !$acc loop seq
         335           0 :            do lcl=1,ncd
         336           0 :               all_gmnnn = all_gmnnn+gmnnn0(i,k,lcl)
         337           0 :               all_lmnnn = all_lmnnn+lmnnn0(i,k,lcl)
         338             :            end do
         339             :
         340           0 :            if ( (all_gmnnn == 0._r8) .or. (all_lmnnn == 0._r8) ) then
         341             :               !$acc loop seq
         342           0 :               do lcl=1,ncd
         343           0 :                  gmnnn0(i,k,lcl) = 0._r8
         344           0 :                  lmnnn0(i,k,lcl) = 0._r8
         345             :               end do
         346             :            else
         347             :               !$acc loop seq
         348           0 :               do lcl=1,ncd
         349           0 :                  lmnnn0(i,k,lcl) = lmnnn0(i,k,lcl)*(all_gmnnn/all_lmnnn)
         350             :               end do
         351             :            end if
         352             :
         353             :            !$acc loop seq
         354           0 :            do lcl=1,ncd
         355           0 :               amk0(i,k,lcl) = amk0(i,k,lcl)+(gmnnn0(i,k,lcl)-lmnnn0(i,k,lcl))*1.e3_r8/ &
         356           0 :                           rho(i,k)*deltatin/dble(sub_step)
         357             :               ank0(i,k,lcl) = ank0(i,k,lcl)+(gnnnn0(i,k,lcl)-lnnnn0(i,k,lcl))*1.e6_r8/ &
         358           0 :                           rho(i,k)*deltatin/dble(sub_step)
         359           0 :               gmnnn(i,k,lcl) = gmnnn(i,k,lcl)+gmnnn0(i,k,lcl)/sub_step
         360           0 :               gnnnn(i,k,lcl) = gnnnn(i,k,lcl)+gnnnn0(i,k,lcl)/sub_step
         361           0 :               lmnnn(i,k,lcl) = lmnnn(i,k,lcl)+lmnnn0(i,k,lcl)/sub_step
         362           0 :               lnnnn(i,k,lcl) = lnnnn(i,k,lcl)+lnnnn0(i,k,lcl)/sub_step
         363             :            end do
         364             :         end do ! end of loop "sub_step"
         365             :
         366             :         ! cloud water
         367             :         !$acc loop seq
         368           0 :         do lcl=1,cutoff(i,k)
         369           0 :            qcin_old(i,k) = qcin_old(i,k)+amk(i,k,lcl)
         370           0 :            ncin_old(i,k) = ncin_old(i,k)+ank(i,k,lcl)
         371           0 :            qcin_new(i,k) = qcin_new(i,k)+(gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
         372           0 :            ncin_new(i,k) = ncin_new(i,k)+(gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
         373           0 :            qctend_TAU(i,k) = qctend_TAU(i,k)+(amk0(i,k,lcl)-amk(i,k,lcl))/deltatin
         374           0 :            nctend_TAU(i,k) = nctend_TAU(i,k)+(ank0(i,k,lcl)-ank(i,k,lcl))/deltatin
         375           0 :            gmnnn_lmnnn_TAU(i,k) = gmnnn_lmnnn_TAU(i,k)+gmnnn(i,k,lcl)-lmnnn(i,k,lcl)
         376             :         end do
         377             :
         378             :         ! rain
         379             :         !$acc loop seq
         380           0 :         do lcl=cutoff(i,k)+1,ncd
         381           0 :            qrin_old(i,k) = qrin_old(i,k)+amk(i,k,lcl)
         382           0 :            nrin_old(i,k) = nrin_old(i,k)+ank(i,k,lcl)
         383           0 :            qrin_new(i,k) = qrin_new(i,k)+(gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
         384           0 :            nrin_new(i,k) = nrin_new(i,k)+(gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
         385           0 :            qrtend_TAU(i,k) = qrtend_TAU(i,k)+(amk0(i,k,lcl)-amk(i,k,lcl))/deltatin
         386           0 :            nrtend_TAU(i,k) = nrtend_TAU(i,k)+(ank0(i,k,lcl)-ank(i,k,lcl))/deltatin
         387           0 :            gmnnn_lmnnn_TAU(i,k) = gmnnn_lmnnn_TAU(i,k)+gmnnn(i,k,lcl)-lmnnn(i,k,lcl)
         388             :         end do
         389             :
         390             :         !$acc loop seq
         391           0 :         do lcl=1,ncd
         392           0 :            amk_out(i,k,lcl) = amk(i,k,lcl) + (gmnnn(i,k,lcl)-lmnnn(i,k,lcl))*1.e3_r8/rho(i,k)*deltatin
         393           0 :            ank_out(i,k,lcl) = ank(i,k,lcl) + (gnnnn(i,k,lcl)-lnnnn(i,k,lcl))*1.e6_r8/rho(i,k)*deltatin
         394             :         end do
         395             :
         396           0 :         qcin_new(i,k) = qcin_new(i,k)+qcin_old(i,k)
         397           0 :         ncin_new(i,k) = ncin_new(i,k)+ncin_old(i,k)
         398           0 :         qrin_new(i,k) = qrin_new(i,k)+qrin_old(i,k)
         399           0 :         nrin_new(i,k) = nrin_new(i,k)+nrin_old(i,k)
         400             :      end do
         401             :   end do
         402             :   !$acc end parallel
         403             :
         404             : ! Conservation checks
         405             : ! AG: Added May 2023
         406             :
         407             :   !$acc parallel vector_length(VLENS) default(present)
         408             :   !$acc loop gang vector collapse(2)
         409           0 :   do k=1,nlev
         410           0 :      do i=1,mgncol
         411             :
         412             :         ! First make sure all not negative
         413           0 :         qcin_new(i,k)=max(qcin_new(i,k),0._r8)
         414           0 :         ncin_new(i,k)=max(ncin_new(i,k),0._r8)
         415           0 :         qrin_new(i,k)=max(qrin_new(i,k),0._r8)
         416           0 :         nrin_new(i,k)=max(nrin_new(i,k),0._r8)
         417             :
         418             : ! Now adjust so that sign is correct. qc_new,nc_new <= input, qr_new >= input
         419             : ! NOte that due to self collection nr can be larger or smaller than input....
         420             : ! Makes above check redundant I think.
         421             :
         422           0 :         qcin_new(i,k)=min(qcin_new(i,k),qcin(i,k))
         423           0 :         ncin_new(i,k)=min(ncin_new(i,k),ncin(i,k))
         424           0 :         qrin_new(i,k)=max(qrin_new(i,k),qrin(i,k))
         425             :
         426             : ! Next scale mass...so output qc+qr is the same as input
         427             :
         428           0 :         if ( (qcin_new(i,k)+qrin_new(i,k)) > 0._r8 ) then
         429           0 :            qscl = (qcin(i,k)+qrin(i,k))/(qcin_new(i,k)+qrin_new(i,k))
         430             :         else
         431             :            qscl = 0._r8
         432             :         end if
         433           0 :         qcin_new(i,k) = qcin_new(i,k) * qscl
         434           0 :         qrin_new(i,k) = qrin_new(i,k) * qscl
         435             :
         436             : ! Now zero nr,nc if either small or no mass?
         437             :
         438           0 :         if ( qcin_new(i,k) < qsmall ) then
         439           0 :            ncin_new(i,k) = 0._r8
         440             :         end if
         441             :
         442           0 :         if ( qrin_new(i,k) < qsmall ) then
         443           0 :            nrin_new(i,k) = 0._r8
         444             :         end if
         445             :
         446             : !Finally add number if mass but no (or super small) number
         447             :
         448           0 :         if ( qcin_new(i,k) > qsmall .and. ncin_new(i,k) < qsmall ) then
         449           0 :            ncin_new(i,k) = qcin_new(i,k)/m1
         450             :         end if
         451             :
         452           0 :         if ( qrin_new(i,k) > qsmall .and. nrin_new(i,k) < qsmall) then
         453           0 :            nrin_new(i,k) = qrin_new(i,k)/m1
         454             :         end if
         455             :
         456             : ! Then recalculate tendencies based on difference
         457             : ! Clip tendencies for cloud (qc,nc) to be <= 0.
         458             : ! Qrtend is not used in pumas (-qctend is used) but clip that too).
         459             : ! Nr tend can be muliply signed.
         460             :
         461           0 :         qctend_TAU(i,k)= min((qcin_new(i,k) - qcin(i,k)) / deltatin,0._r8)
         462           0 :         nctend_TAU(i,k)= min((ncin_new(i,k) - ncin(i,k)) / deltatin,0._r8)
         463           0 :         qrtend_TAU(i,k)= max((qrin_new(i,k) - qrin(i,k)) / deltatin,0._r8)
         464           0 :         nrtend_TAU(i,k)= (nrin_new(i,k) - nrin(i,k)) / deltatin
         465             :
         466             :      end do
         467             :   end do
         468             :   !$acc end parallel
         469             :
         470             :   !$acc end data
         471             :
         472           0 : end subroutine pumas_stochastic_collect_tau_tend
         473             :
         474           0 : subroutine cam_bin_distribute(qc_all,qr_all,qc,nc,qr,nr,mu_c,lambda_c, &
         475           0 :                               lambda_r,n0r,lcldm,precip_frac,scale_qc, &
         476           0 :                               scale_nc,scale_qr,scale_nr,amk_c,ank_c,  &
         477           0 :                               amk_r,ank_r,amk,ank,cutoff_amk,mgncol,nlev)
         478             :
         479             :   implicit none
         480             :
         481             :   integer, intent(in) :: mgncol,nlev
         482             :   real(r8), dimension(mgncol,nlev), intent(in) :: qc_all,qr_all,qc,nc,qr,nr,mu_c, &
         483             :                                                   lambda_c,lambda_r,n0r,lcldm,    &
         484             :                                                   precip_frac
         485             :   real(r8), dimension(mgncol,nlev,ncd), intent(out) :: amk_c,ank_c,amk_r,ank_r,amk,ank
         486             :   real(r8), dimension(mgncol,nlev), intent(out) :: scale_nc,scale_qc,scale_nr,scale_qr
         487             :   integer, dimension(mgncol,nlev), intent(out) :: cutoff_amk
         488             :
         489             :   ! Local variables
         490             :
         491             :   integer  :: i,j,k
         492             :   real(r8) :: phi
         493             :   integer  :: id_max_qc, id_max_qr
         494             :   real(r8) :: max_qc, max_qr, min_amk
         495             :
         496             :   !$acc parallel vector_length(VLENS) default(present)
         497             :   !$acc loop gang vector collapse(3)
         498           0 :   do j=1,ncd
         499           0 :      do k=1,nlev
         500           0 :         do i=1,mgncol
         501           0 :            ank_c(i,k,j) = 0._r8
         502           0 :            amk_c(i,k,j) = 0._r8
         503           0 :            ank_r(i,k,j) = 0._r8
         504           0 :            amk_r(i,k,j) = 0._r8
         505           0 :            ank(i,k,j) = 0._r8
         506           0 :            amk(i,k,j) = 0._r8
         507             :         end do
         508             :      end do
         509             :   end do
         510             :   !$acc end parallel
         511             :
         512             :   !$acc parallel vector_length(VLENS) default(present)
         513             :   !$acc loop gang vector collapse(2)
         514           0 :   do k=1,nlev
         515           0 :      do i=1,mgncol
         516           0 :         scale_nc(i,k) = 0._r8
         517           0 :         scale_qc(i,k) = 0._r8
         518           0 :         scale_nr(i,k) = 0._r8
         519           0 :         scale_qr(i,k) = 0._r8
         520           0 :         cutoff_amk(i,k) = 0
         521             :
         522           0 :         id_max_qc = 0
         523           0 :         id_max_qr = 0
         524           0 :         max_qc = 0._r8
         525           0 :         max_qr = 0._r8
         526             :
         527             :         ! cloud water, nc in #/m3 --> #/cm3
         528           0 :         if ( (qc_all(i,k) > qsmall) .and. (qc(i,k) > qsmall) ) then
         529             :            !$acc loop seq
         530           0 :            do j=1,ncd
         531             :               phi = nc(i,k)*lambda_c(i,k)**(mu_c(i,k)+1._r8)/ &
         532           0 :                     gamma(mu_c(i,k)+1._r8)*(diammean(j)*1.e-2_r8)**mu_c(i,k)* &
         533           0 :                     exp(-lambda_c(i,k)*diammean(j)*1.e-2_r8)  ! D cm --> m
         534           0 :               ank_c(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8   ! D cm --> m
         535           0 :               amk_c(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8*mmean(j)*1.e-3_r8  ! mass in bin g --> kg
         536           0 :               scale_nc(i,k) = scale_nc(i,k)+ank_c(i,k,j)
         537           0 :               scale_qc(i,k) = scale_qc(i,k)+amk_c(i,k,j)
         538             :            end do
         539           0 :            scale_nc(i,k) = scale_nc(i,k)/nc(i,k)
         540           0 :            scale_qc(i,k) = scale_qc(i,k)/qc(i,k)
         541             :
         542             :            !$acc loop seq
         543           0 :            do j=1,ncd
         544           0 :               ank_c(i,k,j) = ank_c(i,k,j)/scale_nc(i,k)*lcldm(i,k)
         545           0 :               amk_c(i,k,j) = amk_c(i,k,j)/scale_qc(i,k)*lcldm(i,k)
         546           0 :               if ( amk_c(i,k,j) > max_qc ) then
         547           0 :                  id_max_qc = j
         548           0 :                  max_qc = amk_c(i,k,j)
         549             :               end if
         550             :            end do
         551             :         end if
         552             :
         553             :         ! rain drop
         554           0 :         if ( (qr_all(i,k) > qsmall) .and. (qr(i,k) > qsmall) ) then
         555             :            !$acc loop seq
         556           0 :            do j=1,ncd
         557           0 :               phi = n0r(i,k)*exp(-lambda_r(i,k)*diammean(j)*1.e-2_r8)   ! D cm --> m
         558           0 :               ank_r(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8   ! D cm --> m
         559           0 :               amk_r(i,k,j) = phi*(diamedge(j+1)-diamedge(j))*1.e-2_r8*mmean(j)*1.e-3_r8
         560           0 :               scale_nr(i,k) = scale_nr(i,k) + ank_r(i,k,j)
         561           0 :               scale_qr(i,k) = scale_qr(i,k) + amk_r(i,k,j)
         562             :            end do
         563           0 :            scale_nr(i,k) = scale_nr(i,k)/nr(i,k)
         564           0 :            scale_qr(i,k) = scale_qr(i,k)/qr(i,k)
         565             :
         566             :            !$acc loop seq
         567           0 :            do j=1,ncd
         568           0 :               ank_r(i,k,j) = ank_r(i,k,j)/scale_nr(i,k)*precip_frac(i,k)
         569           0 :               amk_r(i,k,j) = amk_r(i,k,j)/scale_qr(i,k)*precip_frac(i,k)
         570           0 :               if ( amk_r(i,k,j) > max_qr ) then
         571           0 :                  id_max_qr = j
         572           0 :                  max_qr = amk_r(i,k,j)
         573             :               end if
         574             :            end do
         575             :         end if
         576             :
         577             :         !$acc loop seq
         578           0 :         do j=1,ncd
         579           0 :            amk(i,k,j) = amk_c(i,k,j) + amk_r(i,k,j)
         580           0 :            ank(i,k,j) = ank_c(i,k,j) + ank_r(i,k,j)
         581             :         end do
         582             :
         583           0 :         if ( (id_max_qc > 0) .and. (id_max_qr > 0) ) then
         584           0 :            if ( (max_qc/max_qr < 10._r8) .or. (max_qc/max_qr > 0.1_r8) ) then
         585           0 :               min_amk = amk(i,k,id_max_qc)
         586             :               !$acc loop seq
         587           0 :               do j=id_max_qc,id_max_qr
         588           0 :                  if ( amk(i,k,j) <= min_amk ) then
         589           0 :                     cutoff_amk(i,k) = j
         590           0 :                     min_amk = amk(i,k,j)
         591             :                  end if
         592             :               end do
         593             :            end if
         594             :         end if
         595             :      end do  ! end of loop "mgncol"
         596             :   end do     ! end of loop "nlev"
         597             :   !$acc end parallel
         598             :
         599             : !input: qc,nc,qr,nr, medge (bin edges). May also need # bins?
         600             : !output: amk, ank (mixing ratio and number in each bin)
         601             :
         602             : !this part will take a bit of thinking about.
         603             : !use size distribution parameters (mu, lambda) to generate the values at discrete size points
         604             : !need to also ensure mass conservation
         605             :
         606           0 : end subroutine cam_bin_distribute
         607             :
         608             : ! here are the subroutines called above that actually do the collision-coalescence calculations:
         609             :
         610             : ! The Kernel is from Jerry from many moons ago (included)
         611             :
         612             : ! I read in the file data and multiply by the summed mass of the individual bins
         613             : ! (with a factor of 1.5 so that the values represent the middle of the bin
         614             :
         615             : ! 941 FORMAT(2X,E12.5)
         616             : !     READ(iunit,941) KNN(IDD,JDD)
         617             : !     KNN(IDD,JDD)=(XK_GR(IDD)*kkfac+XK_GR(JDD)*kkfac)*KNN(IDD,JDD)
         618             :
         619             : !where idd and jdd are the indexes for the bins and xk_gr is the mass of drops in a bin in grams
         620             : !
         621             :
         622             : !************************************************************************************
         623             : ! Setup variables needed for collection
         624             : ! Either pass in or define globally the following variables
         625             : ! tbase(height) - temperature in K as a function of height
         626             : ! rhon(height) - air density as a function of height in kg/m^3
         627             : ! xk_gr(bins) - mass of single drop in each bin in grams
         628             : ! lsmall - small number
         629             : ! QC - mass mixing ratio in kg/kg
         630             : ! QN - number mixing ratio in #/kg
         631             : ! All parameters are defined to be global in my version so that they are readily available throughout the code:
         632             : ! SMN0,SNN0,SMCN,APN,AMN2,AMN3,PSIN,FN,FPSIN,XPSIN,HPSIN,FN2,XXPSIN (all arrays of drop bins)
         633             : !************************************************************************************
         634             :
         635             : !AG: Global arrays need to be passed around I think? Right now at the module level. Is that okay?
         636             :
         637           0 : SUBROUTINE COMPUTE_COLL_PARAMS(rhon,xk_gr,qc,qn,gnnnn,gmnnn,lnnnn,lmnnn)
         638             :
         639             :   !$acc routine seq
         640             :
         641             :   IMPLICIT NONE
         642             :
         643             : ! variable declarations (added by hm, 020118)
         644             : ! note: vertical array looping is stripped out, this subroutine operates
         645             : ! only on LOCAL values
         646             :
         647             :   real(r8), dimension(ncd) :: qc,qn
         648             :   real(r8), dimension(ncdp) :: xk_gr
         649             :   real(r8) :: tbase,rhon
         650             :   integer :: lk
         651             :   integer :: l
         652             :   real(r8), parameter :: lsmall = 1.e-12_r8
         653             :   real(r8), dimension(ncd) :: smn0,snn0,smcn,amn2,amn3,psin,fn,fpsin, &
         654             :                                xpsin,hpsin,fn2,xxpsin
         655             :   real(r8) :: apn
         656             :
         657             :   real(r8), dimension(ncd) :: gnnnn,gmnnn,lnnnn,lmnnn
         658             :   integer :: lm1,ll
         659             :
         660           0 :   lk=ncd
         661             :
         662           0 :   DO L=1,LK
         663           0 :      SMN0(L)=QC(L)*RHON/1.E3_r8
         664           0 :      SNN0(L)=QN(L)*RHON/1.E6_r8
         665             :
         666           0 :      IF(SMN0(L).LT.lsmall.OR.SNN0(L).LT.lsmall)THEN
         667           0 :         SMN0(L)=0.0_r8
         668           0 :         SNN0(L)=0.0_r8
         669             :      ENDIF
         670             :   ENDDO
         671             :
         672           0 :   DO L=1,LK
         673           0 :      IF(SMN0(L) .gt. 0._r8.AND.SNN0(L) .gt. 0._r8)THEN
         674           0 :         SMCN(L)=SMN0(L)/SNN0(L)
         675           0 :         IF((SMCN(L) .GT. 2._r8*XK_GR(L)))THEN
         676           0 :            SMCN(L) = (2._r8*XK_GR(L))
         677             :         ENDIF
         678           0 :         IF((SMCN(L) .LT. XK_GR(L)))THEN
         679           0 :            SMCN(L) = XK_GR(L)
         680             :         ENDIF
         681             :      ELSE
         682           0 :         SMCN(L)=0._r8
         683             :      ENDIF
         684           0 :      IF (SMCN(L).LT.XK_GR(L).OR.SMCN(L).GT.(2._r8*XK_GR(L)).OR.SMCN(L).EQ.0.0_r8)THEN
         685             :         APN=1.0_r8
         686             :      ELSE
         687           0 :         APN=0.5_r8*(1._r8+3._r8*(XK_GR(L)/SMCN(L))-2._r8*((XK_GR(L)/SMCN(L))**2._r8))
         688             :      ENDIF
         689             :
         690           0 :      IF(SNN0(L) .GT. LSMALL)THEN
         691           0 :         AMN2(L)=APN*SMN0(L)*SMN0(L)/SNN0(L)
         692           0 :         AMN3(L)=APN*APN*APN*SMN0(L)*SMN0(L)*SMN0(L)/(SNN0(L)*SNN0(L))
         693             :      ELSE
         694           0 :         AMN2(L)=0._r8
         695           0 :         AMN3(L)=0._r8
         696             :      ENDIF
         697             :
         698           0 :      IF(SMCN(L).LT.XK_GR(L))THEN
         699           0 :         PSIN(L)=0.0_r8
         700           0 :         FN(L)=2._r8*SNN0(L)/XK_GR(L)
         701             :      ELSE
         702           0 :         IF(SMCN(L).GT.(2._r8*XK_GR(L)))THEN
         703           0 :            FN(L)=0.0_r8
         704           0 :            PSIN(L)=2._r8*SNN0(L)/XK_GR(L)
         705             :         ELSE
         706           0 :            PSIN(L)=2._r8/XK_GR(L)*(SMN0(L)/XK_GR(L)-SNN0(L))
         707           0 :            FN(L)=2._r8/XK_GR(L)*(2._r8*SNN0(L)-SMN0(L)/XK_GR(L))
         708             :         ENDIF
         709             :      ENDIF
         710             :
         711           0 :      IF(SNN0(L).LT.LSMALL.OR.SMN0(L).LT.LSMALL)THEN
         712           0 :         PSIN(L)=0.0_r8
         713           0 :         FN(L)=0.0_r8
         714             :      ENDIF
         715             :
         716           0 :      FPSIN(L)=0.5_r8/XK_GR(L)*(PSIN(L)-FN(L))
         717           0 :      XPSIN(L)=2._r8*XK_GR(L)*PSIN(L)
         718           0 :      HPSIN(L)=PSIN(L)-0.5_r8*FN(L)
         719           0 :      FN2(L)=FN(L)/2._r8
         720             :
         721           0 :      IF(L.GT.1)THEN
         722           0 :         XXPSIN(L)=XK_GR(L)*PSIN(L-1)
         723             :      ENDIF
         724             :   ENDDO
         725             :
         726             : !************************************************************************************
         727             : ! Compute collision coalescence
         728             : ! Either pass in or define globally the following variables
         729             : ! Gain terms begin with G, loss terms begin with L
         730             : ! Second letter defines mass (M) or number (N)
         731             : ! Third and fourth letters define the types of particles colling, i.e., NN means drops colliding with drops
         732             : ! Last letter defines the category the new particles go into, in this case just N for liquid drops
         733             : ! The resulting rates are in units of #/cm^3/s and g/cm^3/s
         734             : ! Relies on predefined kernel array KNN(bins,bins) - see top of this file
         735             : !************************************************************************************
         736             :
         737           0 :    GMNNN = 0._r8
         738           0 :    GNNNN = 0._r8
         739           0 :    LMNNN = 0._r8
         740           0 :    LNNNN = 0._r8
         741             : ! remove verical array index, calculate gain/loss terms locally
         742             :
         743           0 :   DO L=3,LK-1
         744           0 :      LM1=L-1
         745           0 :      DO LL=1,L-2
         746           0 :         GNNNN(L)=GNNNN(L)+(PSIN(LM1)*SMN0(LL)-FPSIN(LM1)*AMN2(LL))*KNN(LM1,LL)
         747           0 :         GMNNN(L)=GMNNN(L)+(XK_GR(L)*PSIN(LM1)*SMN0(LL)+FN2(LM1)*AMN2(LL)-FPSIN(LM1)*AMN3(LL))*KNN(LM1,LL)
         748             :      ENDDO
         749             :   ENDDO
         750             :
         751           0 :   DO L=2,LK-1
         752           0 :      LM1=L-1
         753           0 :      GNNNN(L)=GNNNN(L)+0.5_r8*SNN0(LM1)*SNN0(LM1)*KNN(LM1,LM1)
         754           0 :      GMNNN(L)=GMNNN(L)+0.5_r8*(SNN0(LM1)*SMN0(LM1)+SMN0(LM1)*SNN0(LM1))*KNN(LM1,LM1)
         755           0 :      DO LL=1,L-1
         756           0 :         LNNNN(L)=LNNNN(L)+(PSIN(L)*SMN0(LL)-FPSIN(L)*AMN2(LL))*KNN(L,LL)
         757           0 :         GMNNN(L)=GMNNN(L)+(SMN0(LL)*SNN0(L)-PSIN(L)*AMN2(LL)+FPSIN(L)*AMN3(LL))*KNN(L,LL)
         758           0 :         LMNNN(L)=LMNNN(L)+(XPSIN(L)*SMN0(LL)-HPSIN(L)*AMN2(LL))*KNN(L,LL)
         759             :      ENDDO
         760             :   ENDDO
         761             :
         762           0 :   DO L=1,LK-1
         763           0 :      DO LL=L,LK-1
         764           0 :         LNNNN(L)=LNNNN(L)+(SNN0(LL)*SNN0(L))*KNN(LL,L)
         765           0 :         LMNNN(L)=LMNNN(L)+(SNN0(LL)*SMN0(L))*KNN(LL,L)
         766             :      ENDDO
         767             :   ENDDO
         768             :
         769           0 : END SUBROUTINE COMPUTE_COLL_PARAMS
         770             :
         771             :
         772             : end module pumas_stochastic_collect_tau
         773             :
         774             :

Generated by: LCOV version 1.14