!----------------------------------------------------------------------- ! Copyright 2011 NOAA Geophysical Fluid Dynamics Lab, Princeton, NJ ! This program is distributed under the terms of the GNU General Public ! License. See the file COPYING contained in this directory !----------------------------------------------------------------------- module moisture_convert_mod !----------------------------------------------------------------------- use plev_constants_mod, only: RDGAS, RVGAS, TFREEZE implicit none private public sphum_to_mixrat, mixrat_to_sphum, & e_to_mixrat, e_to_sphum, & mixrat_to_rh, sphum_to_rh, & rh_to_mixrat, rh_to_sphum !----------------------------------------------------------------------- interface sphum_to_mixrat module procedure sphum_to_mixrat_0d, sphum_to_mixrat_3d end interface interface mixrat_to_sphum module procedure mixrat_to_sphum_0d, mixrat_to_sphum_3d end interface interface e_to_mixrat module procedure e_to_mixrat_0d, e_to_mixrat_3d end interface interface e_to_sphum module procedure e_to_sphum_0d, e_to_sphum_3d end interface interface mixrat_to_rh module procedure mixrat_to_rh_0d, mixrat_to_rh_3d end interface interface sphum_to_rh module procedure sphum_to_rh_0d, sphum_to_rh_3d end interface interface rh_to_mixrat module procedure rh_to_mixrat_0d, rh_to_mixrat_3d end interface interface rh_to_sphum module procedure rh_to_sphum_0d, rh_to_sphum_3d end interface interface lookup_es module procedure lookup_es_0d, lookup_es_3d end interface !----------------------------------------------------------------------- real, parameter :: eps = RDGAS/RVGAS, one_eps = 1.-eps !----------------------------------------------------------------------- ! parameters for es table size and resolution integer, parameter :: tcmin = -160 ! minimum temperature (degC) in lookup table integer, parameter :: tcmax = 100 ! maximum temperature (degC) in lookup table integer, parameter :: esres = 10 ! table resolution (increments per degree) integer, parameter :: nsize = (tcmax-tcmin)*esres+1 ! lookup table size integer, parameter :: nlim = nsize-1 real :: tmin, tmax ! lookup table limits in degK real :: dtres, dtinv, teps real :: TABLE(nsize) ! sat vapor pres (es) real :: DTABLE(nsize) ! first derivative of es real :: D2TABLE(nsize) ! second derivative of es logical :: table_is_initialized = .FALSE. contains !####################################################################### function sphum_to_mixrat_0d (sphum) result (mixrat) real, intent(in) :: sphum real :: mixrat mixrat = sphum / (1.0 - sphum) end function sphum_to_mixrat_0d !----------------------------------------------------------------------- function sphum_to_mixrat_3d (sphum) result (mixrat) real, intent(in) :: sphum(:,:,:) real, dimension(size(sphum,1),size(sphum,2),size(sphum,3)) :: mixrat mixrat = sphum / (1.0 - sphum) end function sphum_to_mixrat_3d !####################################################################### function mixrat_to_sphum_0d (mixrat) result (sphum) real, intent(in) :: mixrat real :: sphum sphum = mixrat / (1.0 + mixrat) end function mixrat_to_sphum_0d !----------------------------------------------------------------------- function mixrat_to_sphum_3d (mixrat) result (sphum) real, intent(in) :: mixrat(:,:,:) real, dimension(size(mixrat,1),size(mixrat,2),size(mixrat,3)) :: & sphum sphum = mixrat / (1.0 + mixrat) end function mixrat_to_sphum_3d !####################################################################### function e_to_mixrat_0d (e, pres) result (mixrat) real, intent(in) :: e, pres real :: mixrat mixrat = eps * e / (pres - e ) end function e_to_mixrat_0d !----------------------------------------------------------------------- function e_to_mixrat_3d (e, pres) result (mixrat) real, intent(in), dimension(:,:,:) :: e, pres real, dimension(size(e,1),size(e,2),size(e,3)) :: mixrat mixrat = eps * e / (pres - e ) end function e_to_mixrat_3d !####################################################################### function e_to_sphum_0d (e, pres) result (sphum) real, intent(in) :: e, pres real :: sphum sphum = eps * e / (pres - one_eps * e ) end function e_to_sphum_0d !----------------------------------------------------------------------- function e_to_sphum_3d (e, pres) result (sphum) real, intent(in), dimension(:,:,:) :: e, pres real, dimension(size(e,1),size(e,2),size(e,3)) :: sphum sphum = eps * e / (pres - one_eps * e ) end function e_to_sphum_3d !####################################################################### function mixrat_to_rh_0d (mixrat, temp, pres) result (rh) real, intent(in) :: mixrat, temp, pres real :: rh, esat, mixrat_sat call lookup_es (temp, esat) mixrat_sat = e_to_mixrat (esat, pres) rh = mixrat / mixrat_sat end function mixrat_to_rh_0d !----------------------------------------------------------------------- function mixrat_to_rh_3d (mixrat, temp, pres) result (rh) real, intent(in), dimension(:,:,:) :: mixrat, temp, pres real, dimension(size(temp,1),size(temp,2),size(temp,3)) :: & rh, esat, mixrat_sat call lookup_es (temp, esat) mixrat_sat = e_to_mixrat (esat, pres) rh = mixrat / mixrat_sat end function mixrat_to_rh_3d !####################################################################### function sphum_to_rh_0d (sphum, temp, pres) result (rh) real, intent(in) :: sphum, temp, pres real :: rh, mixrat mixrat = sphum_to_mixrat (sphum) rh = mixrat_to_rh (mixrat, temp, pres) end function sphum_to_rh_0d !----------------------------------------------------------------------- function sphum_to_rh_3d (sphum, temp, pres) result (rh) real, intent(in), dimension(:,:,:) :: sphum, temp, pres real, dimension(size(temp,1),size(temp,2),size(temp,3)) :: & rh, mixrat mixrat = sphum_to_mixrat (sphum) rh = mixrat_to_rh (mixrat, temp, pres) end function sphum_to_rh_3d !####################################################################### function rh_to_mixrat_0d (rh, temp, pres) result (mixrat) real, intent(in) :: rh, temp, pres real :: mixrat, esat, mixrat_sat call lookup_es (temp, esat) mixrat_sat = e_to_mixrat (esat, pres) mixrat = mixrat_sat * rh end function rh_to_mixrat_0d !----------------------------------------------------------------------- function rh_to_mixrat_3d (rh, temp, pres) result (mixrat) real, intent(in), dimension(:,:,:) :: rh, temp, pres real, dimension(size(temp,1),size(temp,2),size(temp,3)) :: & mixrat, esat, mixrat_sat call lookup_es (temp, esat) mixrat_sat = e_to_mixrat (esat, pres) mixrat = mixrat_sat * rh end function rh_to_mixrat_3d !####################################################################### function rh_to_sphum_0d (rh, temp, pres) result (sphum) real, intent(in) :: rh, temp, pres real :: sphum, mixrat mixrat = rh_to_mixrat (rh, temp, pres) sphum = mixrat_to_sphum (mixrat) end function rh_to_sphum_0d !----------------------------------------------------------------------- function rh_to_sphum_3d (rh, temp, pres) result (sphum) real, intent(in), dimension(:,:,:) :: rh, temp, pres real, dimension(size(temp,1),size(temp,2),size(temp,3)) :: & sphum, mixrat mixrat = rh_to_mixrat (rh, temp, pres) sphum = mixrat_to_sphum (mixrat) end function rh_to_sphum_3d !####################################################################### ! saturation vapor pressure lookup table !####################################################################### subroutine lookup_es_0d (temp, esat) real, intent(in) :: temp real, intent(out) :: esat real :: tmp, del integer :: ind call sat_vapor_pres_init tmp = temp-tmin ind = int(dtinv*(tmp+teps)) del = tmp-dtres*real(ind) esat = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) if (ind < 0 .or. ind > nlim) then print *, 'ERROR: saturation vapor pressure error' stop 111 endif end subroutine lookup_es_0d !----------------------------------------------------------------------- subroutine lookup_es_3d (temp, esat) real, intent(in) :: temp(:,:,:) real, intent(out) :: esat(:,:,:) real :: tmp, del integer :: ind integer :: i, j, k call sat_vapor_pres_init do k = 1, size(temp,3) do j = 1, size(temp,2) do i = 1, size(temp,1) tmp = temp(i,j,k)-tmin ind = int(dtinv*(tmp+teps)) del = tmp-dtres*real(ind) esat(i,j,k) = TABLE(ind+1) + del*(DTABLE(ind+1) + del*D2TABLE(ind+1)) if (ind < 0 .or. ind > nlim) then print *, 'ERROR: saturation vapor pressure error' stop 111 endif enddo enddo enddo end subroutine lookup_es_3d !----------------------------------------------------------------------- ! Computation of the es values ! ! Saturation vapor pressure (es) values are computed from ! equations in the Smithsonian meteorological tables page 350. ! For temperatures < 0C, sat vapor pres is computed over ice. ! For temperatures > -20C, sat vapor pres is computed over water. ! Between -20C and 0C the returned value is blended (over water ! and over ice). All sat vapor pres values are returned in pascals. ! ! Reference: Smithsonian meteorological tables, page 350. function compute_es (tem) result (es) real, intent(in) :: tem(:) real :: es(size(tem,1)) real, parameter :: TBASW = TFREEZE+100. real, parameter :: TBASI = TFREEZE real, parameter :: ESBASW = 101324.60 real, parameter :: ESBASI = 610.71 real :: x, esice, esh2o integer :: i do i = 1, size(tem(:)) ! compute es over ice if (tem(i) < TBASI) then x = -9.09718*(TBASI/tem(i)-1.0) - 3.56654*log10(TBASI/tem(i)) & +0.876793*(1.0-tem(i)/TBASI) + log10(ESBASI) esice =10.**(x) else esice = 0. endif ! compute es over water greater than -20 c. ! values over 100 c may not be valid ! see smithsonian meteorological tables page 350. if (tem(i) > -20.+TBASI) then x = -7.90298*(TBASW/tem(i)-1) + 5.02808*log10(TBASW/tem(i)) & -1.3816e-07*(10**((1-tem(i)/TBASW)*11.344)-1) & +8.1328e-03*(10**((TBASW/tem(i)-1)*(-3.49149))-1) & +log10(ESBASW) esh2o = 10.**(x) else esh2o = 0. endif ! derive blended es over ice and supercooled water between -20c and 0c if (tem(i) <= -20.+TBASI) then es(i) = esice else if (tem(i) >= TBASI) then es(i) = esh2o else es(i) = 0.05*((TBASI-tem(i))*esice + (tem(i)-TBASI+20.)*esh2o) endif enddo end function compute_es !----------------------------------------------------------------------- subroutine sat_vapor_pres_init real :: tem(3), es(3), hdtinv integer :: i ! increment used to generate derivative table real, parameter :: tinrc = .01 real, parameter :: tfact = 1./(2.*tinrc) ! return silently if this routine has already been called if (table_is_initialized) return ! global variables tmin = real(tcmin)+TFREEZE ! minimum valid temp in table tmax = real(tcmax)+TFREEZE ! maximum valid temp in table dtinv = real(esres) dtres = 1./dtinv teps = 1./real(2*esres) ! local variables hdtinv = dtinv*0.5! global variables tmin = real(tcmin)+TFREEZE ! minimum valid temp in table tmax = real(tcmax)+TFREEZE ! maximum valid temp in table dtinv = real(esres) dtres = 1./dtinv teps = 1./real(2*esres) ! local variables hdtinv = dtinv*0.5 ! compute es tables from tcmin to tcmax ! estimate es derivative with small +/- difference do i = 1, nsize tem(1) = tmin + dtres*real(i-1) tem(2) = tem(1)-tinrc tem(3) = tem(1)+tinrc es = compute_es (tem) TABLE(i) = es(1) DTABLE(i) = (es(3)-es(2))*tfact enddo ! compute one-half second derivative using centered differences ! differencing des values in the table do i = 2, nsize-1 D2TABLE(i) = 0.25*dtinv*(DTABLE(i+1)-DTABLE(i-1)) enddo ! one-sided derivatives at boundaries D2TABLE(1) = 0.50*dtinv*(DTABLE(2) -DTABLE(1)) D2TABLE(nsize) = 0.50*dtinv*(DTABLE(nsize)-DTABLE(nsize-1)) table_is_initialized = .true. end subroutine sat_vapor_pres_init end module moisture_convert_mod