MODULE simple_pdf_mod !MNS NOTE: THIS IS EXPERIMENTAL... use fms_mod, only: write_version_number use beta_dist_mod, only: incomplete_beta, beta_dist_init, & beta_dist_end use strat_cloud_utilities_mod, only: strat_cloud_utilities_init, & diag_id_type, diag_pt_type implicit none private !----------------------------------------------------------------------- !----interfaces-------------------------------------------------------- public simple_pdf, simple_pdf_init, simple_pdf_end !---------------------------------------------------------------------- !----version number---------------------------------------------------- Character(len=128) :: Version = '$Id: simple_pdf.F90,v 20.0 2013/12/13 23:22:07 fms Exp $' Character(len=128) :: Tagname = '$Name: tikal $' logical :: module_is_initialized = .false. CONTAINS !######################################################################### SUBROUTINE simple_pdf_init !----------------------------------------------------------------------- if (module_is_initialized) return !----------------------------------------------------------------------- ! write version number to output file. !----------------------------------------------------------------------- call write_version_number(version, tagname) !----------------------------------------------------------------------- ! make sure needed modules have been initialized. !----------------------------------------------------------------------- call strat_cloud_utilities_init call beta_dist_init !----------------------------------------------------------------------- ! mark this module initialized. !----------------------------------------------------------------------- module_is_initialized = .true. !------------------------------------------------------------------------ END SUBROUTINE simple_pdf_init !######################################################################### SUBROUTINE simple_pdf (j, idim, jdim, kdim, qmin, qa, qtot, qs, gamma, & qthalfwidth, betaP, inv_dtcloud, SA_0, n_diag_4d, & diag_4d, diag_id, diag_pt, SA, qa_upd) !------------------------------------------------------------------------- INTEGER, INTENT(IN) :: j, idim, jdim, kdim REAL, INTENT(IN) :: qmin, inv_dtcloud REAL, dimension(idim,kdim), INTENT(IN) :: qtot, qa, qs, & gamma, SA_0 INTEGER, INTENT(IN) :: betaP REAL, INTENT(IN) :: qthalfwidth INTEGER, INTENT(IN) :: n_diag_4d REAL, dimension(idim, jdim, kdim, 0:n_diag_4d), & INTENT(INOUT) :: diag_4d TYPE(diag_id_type), intent(in) :: diag_id TYPE(diag_pt_type), intent(inout) :: diag_pt REAL, dimension(idim, kdim), INTENT(OUT) :: SA, qa_upd !------------------------------------------------------------------------ !local variables REAL, dimension(idim, kdim) :: qta3, qtqsa3, qcg, qag, qa1, qa0, & qabar real, dimension(idim) :: qtbar, deltaQ, qtmin, qs_norm, & qagtmp, qcgtmp, qvgtmp real :: icbp, icbp1, pnorm INTEGER :: k, id !------------------------------------------------------------------------ ! compute pdf cloud fraction and condensate. Note that the SYMMETRIC ! beta distribution is used here. Initialize grid-box mean values of ! cloud fraction (qag) and cloud condensate(qcg). !------------------------------------------------------------------------ qta3(:,:) = max(qmin, qtot) qtqsa3(:,:) = qta3(:,:) - qs qcg(:,:) = 0. qag(:,:) = 0. DO k=1,kdim !------------------------------------------------------------------------ ! calculate normalized vertical level. 0. = top of gridbox, 1. = bottom ! of gridbox !------------------------------------------------------------------------ pnorm = 0.5 !------------------------------------------------------------------------ ! First step is to calculating the minimum (qtmin) of the total water ! distribution and the width of the qt distribution (deltaQ). For ! diagnostic variance this is set to (1.-qthalfwidth)*qtbar and ! 2*qthalfwidth*qtbar, respectively, where qtbar is the mean total water ! in the grid box. !------------------------------------------------------------------------ qtbar = qta3(:,k) qtbar = max(qmin ,qtbar) deltaQ = 2.*qthalfwidth*qtbar qtmin = (1. - qthalfwidth)*qtbar !------------------------------------------------------------------------ ! From this the variable normalized saturation specific humidity qs_norm ! is calculated. ! qs_norm = (qs(Tl) - qtmin)/(qtmax-qtmin) ! ! = 0.5 - (qtbar - qs(Tl))/deltaQ ! ! Note that if qs_norm > 1., the grid box is fully clear. If qs_norm ! < 0., the grid box is fully cloudy. !------------------------------------------------------------------------ qs_norm = qtqsa3(:,k) qs_norm = 0.5 - (qs_norm/deltaQ) !------------------------------------------------------------------------ ! Calculation of cloud fraction (qagtmp), cloud condensate (qcgtmp), and ! water vapor in clear air part of the grid box (qvgtmp). ! ! Formulas (from Tompkins, and personal derivations): ! ! Define icbp = incomplete_beta(qs_norm,p,q) ! icbp1 = incomplete_beta(qs_norm,p+1,q) ! ! qagtmp = 1. - icbp ! ! qcgtmp = aThermo * { (qtbar-qtmin)*(1.-icbp1) - ! qs_norm*deltaQ*(1.-icbp ) } ! ! ! qvgtmp = qtmin + (p/(p+q))*(icbp1/icbp)*deltaQ ! ! ! where aThermo = 1./(1.+(L/cp)*dqsdT) ! ! note that in the qvg formula below the factor of 0.5 is equal ! to (p/(p+q)). ! !------------------------------------------------------------------------ do id = 1,idim if (qs_norm(id).le.1.) then icbp = incomplete_beta (max(0., qs_norm(id)), & p = betaP, q = betaP) icbp1 = incomplete_beta (max(0., qs_norm(id)), & p = betaP + 1, q = betaP) qagtmp(id) = 1. - icbp qcgtmp(id) = (qtbar(id) - qtmin(id))*(1. - icbp1) - & qs_norm(id)*deltaQ(id)*(1. - icbp) qcgtmp(id) = qcgtmp(id)/(1. + gamma(id,k)) qvgtmp(id) = qtmin(id) + 0.5*(icbp1/max(icbp, qmin))*deltaQ(id) !------------------------------------------------------------------------ ! bound very very small cloud fractions which may cause negative cloud ! condensates due to roundoff errors or similar errors in the beta table ! lookup. !------------------------------------------------------------------------ if ((qagtmp(id) .lt. 0.).or.(qcgtmp(id) .le. 0.)) then qagtmp(id) = 0. qcgtmp(id) = 0. qvgtmp(id) = qtbar(id) end if else qagtmp(id) = 0. qcgtmp(id) = 0. qvgtmp(id) = qtbar(id) end if end do !------------------------------------------------------------------------ ! sum vertically. note special averaging of clear-sky water vapor ! this is weighting clear-sky relative humidity by the clear-sky fraction !------------------------------------------------------------------------ qag(:,k) = qag(:,k) + qagtmp qcg(:,k) = qcg(:,k) + qcgtmp enddo ! ( k loop) !------------------------------------------------------------------------ ! compute grid-box average cloud fraction, cloud condensate and ! water vapor. do adjustment of cloud fraction. set total tendency term ! (SA) and update cloud fraction (qa_upd). !------------------------------------------------------------------------ qa0 = qa qa1 = qag SA = SA_0 + qa1 - qa0 qa_upd = qa1 !-------------------------------------------------------------------------- ! fill desired diagnostics. !-------------------------------------------------------------------------- if (diag_id%qadt_lsform + diag_id%qa_lsform_col > 0) then diag_4d(:,j,:,diag_pt%qadt_lsform ) = & max(qa1 - qa0, 0.)*inv_dtcloud end if if (diag_id%qadt_lsdiss + diag_id%qa_lsdiss_col > 0) then diag_4d(:,j,:,diag_pt%qadt_lsdiss ) = & max(qa0 - qa1, 0.)*inv_dtcloud end if !------------------------------------------------------------------------ END SUBROUTINE simple_pdf !######################################################################### SUBROUTINE simple_pdf_end call beta_dist_end module_is_initialized = .false. END SUBROUTINE simple_pdf_end !######################################################################### END MODULE simple_pdf_mod