MODULE polysvp_mod use fms_mod, only : FATAL, error_mesg, write_version_number use sat_vapor_pres_mod, only : compute_qs, sat_vapor_pres_init use constants_mod, only : rdgas, rvgas, hlv, hls, tfreeze, cp_air use mg_const_mod, only : mg_const_init, mg_pr use strat_cloud_utilities_mod, & only : strat_cloud_utilities_init, & atmos_state_type, cloud_state_type,& strat_nml_type IMPLICIT NONE PRIVATE !------------------------------------------------------------------------ !---interfaces----------------------------------------------------------- interface polysvp_l module procedure polysvp_l1d, polysvp_l3d end interface polysvp_l interface polysvp_i module procedure polysvp_i1d, polysvp_i3d end interface polysvp_i PUBLIC polysvp_l, polysvp_i, polysvp_init, polysvp_end, & compute_qs_a, compute_qs_x1 private compute_qs_x2 !------------------------------------------------------------------------- !----version number------------------------------------------------------- Character(len=128) :: Version = '$Id: polysvp.F90,v 20.0 2013/12/13 23:22:03 fms Exp $' Character(len=128) :: Tagname = '$Name: tikal $' !--------------------------------------------------------------------- real, parameter :: d622 = rdgas / rvgas real, parameter :: d378 = 1. - d622 logical :: module_is_initialized = .false. CONTAINS !######################################################################### SUBROUTINE polysvp_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 sat_vapor_pres_init call strat_cloud_utilities_init call mg_const_init !------------------------------------------------------------------------ ! mark this module as initialized. !------------------------------------------------------------------------ module_is_initialized = .TRUE. END SUBROUTINE polysvp_init !######################################################################### SUBROUTINE compute_qs_a (idim, jdim, kdim, Nml, Atmos_state, Cloud_state) !------------------------------------------------------------------------- INTEGER, INTENT(IN ) :: idim, jdim, kdim type(strat_nml_type), intent(in) :: Nml type(atmos_state_type), intent(inout) :: Atmos_state type(cloud_state_type), intent(inout) :: Cloud_state !----------------------------------------------------------------------- !----local variables INTEGER :: i, j, k !------------------------------------------------------------------------- ! compute qs and associated parameters, using different algorithms as ! appropriate. !------------------------------------------------------------------------- if (Nml%do_pdf_clouds .and. Nml%pdf_org) then if (Nml%super_ice_opt .LT. 1) then !------------------------------------------------------------------------- ! if using pdf clouds and not allowing supersaturation: !------------------------------------------------------------------------- call compute_qs (Atmos_state%T_in - ((hlv*Cloud_state%ql_in + & hls*Cloud_state%qi_in)/cp_air), & Atmos_state%pfull, Atmos_state%qs, & dqsdT=Atmos_state%dqsdT, esat=Atmos_state%esat0) Atmos_state%gamma = Atmos_state%dqsdT*(min(1., & max(0., 0.05*(Atmos_state%T_in - ((hlv*Cloud_state%ql_in + & hls*Cloud_state%qi_in)/cp_air) - tfreeze + 20.)))*hlv + & min(1., max(0., 0.05*(tfreeze - Atmos_state%T_in + & ((hlv*Cloud_state%ql_in + & hls*Cloud_state%qi_in)/cp_air))))*hls)/cp_air else !------------------------------------------------------------------------- ! if using pdf clouds and allowing supersaturation: !------------------------------------------------------------------------- call error_mesg ( 'strat_cloud_mod', & 'super_ice_opt > 1 and pdf_opt currently not supported', FATAL) ! this code is pattern for what might be used when this option supported: do k=1,kdim do j=1,jdim do i=1,idim !call compute_qs_x1 (T(i,j,k), pfull(i,j,k), qs(i,j,k), & ! qsl(i,j,k), qsi(i,j,k), dqsdT(i,j,k), & ! gamma(i,j,k)) end do end do end do Atmos_state%gamma = Atmos_state%dqsdT*(min(1., & max(0., 0.05*(Atmos_state%T_in - tfreeze + 20.)))*hlv + & min(1., max(0., 0.05*(tfreeze - & Atmos_state%T_in)))*hls)/cp_air end if else !------------------------------------------------------------------------- ! if not using pdf clouds and not allowing supersaturation: !------------------------------------------------------------------------- if (Nml%super_ice_opt .LT. 1) then !------------------------------------------------------------------------ ! Calculate saturation specific humidity and its temperature ! derivative, and thermal conductivity plus vapor diffusivity factor. ! ! FOR ORIGINAL SCHEME These are calculated according to the formulas: ! ! (1) qs = d622*esat/ [pfull - (1.-d622)*esat] ! ! (2) dqsdT = d622*pfull*(desat/dT)/[pfull-(1.-d622)*esat]**2. ! ! (3) gamma = (L/cp) * dqsdT ! ! where d622 = rdgas/rvgas; esat = saturation vapor pressure; ! and desat/dT is the temperature derivative of esat. ! Note that in the calculation of gamma, ! ! { hlv for T > tfreeze } ! L = { 0.05*(T-tfreeze+20.)*hlv + 0.05*(tfreeze-T)*hls } ! { for tfreeze-20.< T < tfreeze} ! { hls for T < tfreeze-20. } ! ! This linear form is chosen because at tfreeze-20. es = esi, and ! at tfreeze, es = esl, with linear interpolation in between. ! ! ! (5) chi = 2.21 E-05 (m*m)/s * (1.E+05)/pfull ! ! where p is the pressure in Pascals. ! ! ! Note that qs, dqsdT, and gamma do not have their proper values ! until all of the following code has been executed. That ! is qs and dqsdT are used to store intermediary results ! in forming the full solution. !------------------------------------------------------------------------ !------------------------------------------------------------------------ ! calculate water saturated vapor pressure from table and store ! temporarily in the variable esat0. !------------------------------------------------------------------------ call compute_qs & (Atmos_state%T_in, Atmos_state%pfull, Atmos_state%qs,& dqsdT=Atmos_state%dqsdT, esat=Atmos_state%esat0) Atmos_state%gamma = Atmos_state%dqsdT*(min(1., & max(0., 0.05*(Atmos_state%T_in - tfreeze + 20.)))*hlv + & min(1., max(0., 0.05*(tfreeze - Atmos_state%T_in)))*hls)/cp_air else !------------------------------------------------------------------------- ! if not using pdf clouds and allowing supersaturation: !------------------------------------------------------------------------- call compute_qs_x1 & (idim, jdim, kdim, Atmos_state%T_in, Atmos_state%pfull, & Atmos_state%qs, Atmos_state%qsl, Atmos_state%qsi, & Atmos_state%dqsdT, Atmos_state%gamma) end if end if !-------------------------------------------------------------------------- ! Calculate relative humidity outside of convective cloud portion of ! grid box (U_ca) under the assumption that the temp is uniform across ! the gridbox. Define the maximum portion of the remaining grid-box ! area which can be cloudy while maintaining the grid-box-mean relative ! humidity (U01). Here qrat is the ratio of the gridbox RH to that in the ! environment of the convective cloud. !-------------------------------------------------------------------------- IF (Nml%super_ice_opt .LT. 1) THEN !-------------------------------------------------------------------------- ! when supersaturation not permitted, upper limit of 100% placed on ! external relative humidity. !-------------------------------------------------------------------------- where (Atmos_state%qrat .gt. 0.) Atmos_state%U_ca = min(max(0., & (Atmos_state%qv_in/(Atmos_state%qrat*Atmos_state%qs))), 1.) elsewhere Atmos_state%U_ca = 0. end where Atmos_state%U01 = Atmos_state%U_ca ELSE !------------------------------------------------------------------------ ! when supersaturation is permitted, no upper limit placed on external ! relative humidity. Max cloud fraction depends on whether clouds are ! either liquid or ice. !-------------------------------------------------------------------------- do k=1,kdim do j=1,jdim do i=1,idim if (Atmos_state%qrat(i,j,k) .gt. 0. ) then Atmos_state%U_ca(i,j,k) = max(0., & (Atmos_state%qv_in(i,j,k)/(Atmos_state%qrat(i,j,k)* & Atmos_state%qs(i,j,k)))) if (Atmos_State%T_in(i,j,k) .LT. tfreeze) then Atmos_state%U01(i,j,k) = max(0., & (Atmos_state%qv_in(i,j,k)/(Atmos_state%qrat(i,j,k)* & Atmos_state%qsi(i,j,k)))) else Atmos_state%U01(i,j,k) = max(0., & (Atmos_state%qv_in(i,j,k)/(Atmos_state%qrat(i,j,k)* & Atmos_state%qsl(i,j,k)))) endif else Atmos_state%U_ca(i,j,k) = 0. endif end do end do end do END IF !----------------------------------------------------------------------- END SUBROUTINE compute_qs_a !######################################################################### SUBROUTINE polysvp_end module_is_initialized = .FALSE. END SUBROUTINE polysvp_end !######################################################################### function polysvp_l1d (T) !------------------------------------------------------------------------ ! Compute saturation vapor pressure with respect to liquid by using ! function from Goff and Gatch (1946). polysvp returned in units of pa. ! T is input in units of K. !------------------------------------------------------------------------ REAL(kind=mg_pr) :: T,polysvp_l1d REAL(kind=mg_pr) :: dum !------------------------------------------------------------------------- ! Goff Gatch equation, uncertain below -70 C !------------------------------------------------------------------------- polysvp_l1d = 10._mg_pr**(-7.90298_mg_pr*(373.16_mg_pr/t-1._mg_pr) + & 5.02808_mg_pr*log10(373.16_mg_pr/t) - & 1.3816e-7_mg_pr*(10._mg_pr**(11.344_mg_pr*(1._mg_pr - & t/373.16_mg_pr)) - 1._mg_pr) + & 8.1328e-3_mg_pr*(10._mg_pr**(-3.49149_mg_pr* & (373.16_mg_pr/t - 1._mg_pr)) - 1._mg_pr) + & log10(1013.246_mg_pr))*100._mg_pr !------------------------------------------------------------------------- end function polysvp_l1d !########################################################################## function polysvp_i1d (T) !------------------------------------------------------------------------ ! Compute saturation vapor pressure with respect to ice by using ! function from Goff and Gatch (1946). Polysvp returned in units of pa. ! T is input in units of K. !------------------------------------------------------------------------ REAL(kind=mg_pr) :: T,polysvp_i1d REAL(kind=mg_pr) :: dum !----------------------------------------------------------------------- ! Goff Gatch equation for ice (good down to -100 C) !----------------------------------------------------------------------- polysvp_i1d = 10._mg_pr**(-9.09718_mg_pr*(273.16_mg_pr/t - & 1._mg_pr) - 3.56654_mg_pr* & log10(273.16_mg_pr/t) + 0.876793_mg_pr* & (1._mg_pr - t/273.16_mg_pr) + & log10(6.1071_mg_pr))*100._mg_pr !------------------------------------------------------------------------- end function polysvp_i1d !######################################################################## function polysvp_l3d (T, idim, jdim, kdim) !------------------------------------------------------------------------ ! Compute saturation vapor pressure with respect to liquid by using ! function from Goff and Gatch (1946). polysvp returned in units of pa. ! T is input in units of K. !------------------------------------------------------------------------ REAL(kind=mg_pr), dimension(idim,jdim,kdim) :: T,polysvp_l3d REAL(kind=mg_pr) :: dum integer :: idim,jdim,kdim !------------------------------------------------------------------------- ! Goff Gatch equation, uncertain below -70 C !------------------------------------------------------------------------- polysvp_l3d = 10._mg_pr**(-7.90298_mg_pr*(373.16_mg_pr/t-1._mg_pr) + & 5.02808_mg_pr*log10(373.16_mg_pr/t) - & 1.3816e-7_mg_pr*(10._mg_pr**(11.344_mg_pr*(1._mg_pr - & t/373.16_mg_pr)) - 1._mg_pr) + & 8.1328e-3_mg_pr*(10._mg_pr**(-3.49149_mg_pr* & (373.16_mg_pr/t - 1._mg_pr)) - 1._mg_pr) + & log10(1013.246_mg_pr))*100._mg_pr !------------------------------------------------------------------------- end function polysvp_l3d !######################################################################### function polysvp_i3d (T, idim, jdim, kdim) !------------------------------------------------------------------------ ! Compute saturation vapor pressure with respect to ice by using ! function from Goff and Gatch (1946). Polysvp returned in units of pa. ! T is input in units of K. !------------------------------------------------------------------------ REAL(kind=mg_pr), dimension(idim,jdim,kdim) :: T, polysvp_i3d REAL(kind=mg_pr) :: dum integer :: idim, jdim, kdim !----------------------------------------------------------------------- ! Goff Gatch equation for ice (good down to -100 C) !----------------------------------------------------------------------- polysvp_i3d = 10._mg_pr**(-9.09718_mg_pr*(273.16_mg_pr/t - & 1._mg_pr) - 3.56654_mg_pr* & log10(273.16_mg_pr/t) + 0.876793_mg_pr* & (1._mg_pr - t/273.16_mg_pr) + & log10(6.1071_mg_pr))*100._mg_pr !------------------------------------------------------------------------- end function polysvp_i3d !######################################################################### SUBROUTINE compute_qs_x1 (idim, jdim, kdim, ttmp, pfull, qs, qs_l, qs_i,& dqsdT, gamma ) !-------------------------------------------------------------------------- integer, intent(in) :: idim, jdim, kdim REAL, dimension(idim,jdim,kdim), INTENT(IN ) :: ttmp, pfull REAL, dimension(idim,jdim,kdim), INTENT(INOUT ) :: qs, dqsdT, gamma, & qs_l, qs_i !-------------------------------------------------------------------------- !---local variables-------------------------------------------------------- REAL, dimension (idim,jdim,kdim) :: eslt, esit integer :: i,j,k !------------------------------------------------------------------------- ! compute es with respect to liquid and ice. !------------------------------------------------------------------------- eslt = polysvp_l3d (ttmp, idim, jdim, kdim) esit = polysvp_i3d (ttmp, idim, jdim, kdim) !------------------------------------------------------------------------- ! make sure ice saturation doesn't exceed water saturation at temps ! near freezing. !------------------------------------------------------------------------- where (esit .gt. eslt) esit = eslt !------------------------------------------------------------------------- ! compute saturation specific humidity over liquid. calculate denominator ! in qsat formula. limit denominator to esat, and thus qs to d622. this ! is done to avoid blow up in the upper stratosphere where pfull ~ esat. !------------------------------------------------------------------------- qs_l = pfull - d378*eslt qs_l = max(qs_l, eslt) qs_l = d622*eslt/qs_l !------------------------------------------------------------------------- ! compute saturation specific humidity over ice. calculate denominator ! in qsat formula. limit denominator to esat, and thus qs to d622. this ! is done to avoid blow up in the upper stratosphere where pfull ~ esat. ! compute saturation specific humidity over liquid. !------------------------------------------------------------------------- qs_i = pfull - d378*esit qs_i = max(qs_i, esit) qs_i = d622*esit/qs_i !-------------------------------------------------------------------------- ! define the appropriate qs, dqsdT and gamma to use, dependent on ambient ! temperature. !-------------------------------------------------------------------------- do k=1,kdim do j=1,jdim do i=1,idim IF (ttmp(i,j,k) .GE. tfreeze - 23.16) THEN qs(i,j,k) = qs_l(i,j,k) ELSE qs(i,j,k) = qs_i(i,j,k) END IF !REV#1 !888 ! IF (ttmp(i,j,k) .GE. 233.15) THEN ! IF (ttmp(i,j,k) .GE. 233.16) THEN IF (ttmp(i,j,k) .GE. tfreeze - 40.) THEN !END REV#1 dqsdT(i,j,k) = hlv*qs_l(i,j,k)/(rvgas*ttmp(i,j,k)**2) gamma(i,j,k) = dqsdT(i,j,k)*hlv/cp_air ELSE dqsdT(i,j,k) = hls*qs_i(i,j,k)/(rvgas*ttmp(i,j,k)**2) gamma(i,j,k) = dqsdT(i,j,k)*hls/cp_air END IF end do end do end do !------------------------------------------------------------------------ END SUBROUTINE compute_qs_x1 !######################################################################## SUBROUTINE compute_qs_x2 (ttmp, pfull, qs, qs_l, qs_i, dqsdT, gamma, & ql, qi, qmin ) REAL, INTENT(IN ) :: ttmp, pfull, qmin, ql, qi REAL, INTENT(INOUT ) :: qs, dqsdT, gamma, qs_l, qs_i REAL :: eslt, esit, ifrac !------------------------------------------------------------------------- ! determine ice fraction of total condensate. !------------------------------------------------------------------------- if (qi .gt. qmin) then ifrac = qi/(ql + qi) else ifrac = 0. end if !------------------------------------------------------------------------- ! compute es with respect to liquid and ice. !------------------------------------------------------------------------- eslt = polysvp_l1d(ttmp) esit = polysvp_i1d(ttmp) !------------------------------------------------------------------------- ! make sure ice saturation doesn't exceed water saturation at temps ! near freezing. !------------------------------------------------------------------------- IF (esit .GT. eslt) esit = eslt !------------------------------------------------------------------------- ! compute saturation specific humidity over liquid. calculate denominator ! in qsat formula. limit denominator to esat, and thus qs to d622. this ! is done to avoid blow up in the upper stratosphere where pfull ~ esat. !------------------------------------------------------------------------- qs_l = pfull - d378*eslt qs_l = max(qs_l, eslt) qs_l = d622*eslt/qs_l !------------------------------------------------------------------------- ! compute saturation specific humidity over ice. calculate denominator ! in qsat formula. limit denominator to esat, and thus qs to d622. this ! is done to avoid blow up in the upper stratosphere where pfull ~ esat. ! compute saturation specific humidity over liquid. !------------------------------------------------------------------------- qs_i = pfull - d378*esit qs_i = max(qs_i, esit) qs_i = d622*esit/qs_i !-------------------------------------------------------------------------- ! define the appropriate qs, dqsdT and gamma to use, dependent on ambient ! temperature and ice fraction. !-------------------------------------------------------------------------- IF (ttmp .GE. tfreeze - 23.16 .and. ifrac .LT. 0.9 ) THEN qs = qs_l ELSE qs = qs_i END IF !REV#2 !888 ! IF (ttmp .GE. 233.15 .and. ifrac .LT. 0.9 ) THEN ! IF (ttmp .GE. 233.16 .and. ifrac .LT. 0.9 ) THEN IF (ttmp .GE. tfreeze - 40. .and. ifrac .LT. 0.9 ) THEN !END REV#2 dqsdT = hlv*qs_l/(rvgas*ttmp**2) gamma = dqsdT * hlv /cp_air ELSE dqsdT = hls*qs_i/(rvgas*ttmp**2) gamma = dqsdT * hls/ cp_air END IF !------------------------------------------------------------------------- END SUBROUTINE compute_qs_x2 !######################################################################### END MODULE polysvp_mod