Changing cold trap detection

CMakeLists.txt

@@ -7,6 +7,7 @@ set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules")
 
 # includes
 include(cmake/CPM.cmake)
+include(cmake/fypp.cmake)
 
 # options
 option(SKBUILD "Should be ON of being build by skbuild, 

cmake/fypp.cmake

@@ -0,0 +1,35 @@
+#     srcext [in]: File extension of the source files
+#     trgext [in]: File extension of the target files
+#     srcfiles [in]: List of the source files
+#     trgfiles [out]: Contains the list of the preprocessed files on exit
+#
+function(preprocess preproc preprocopts srcext trgext srcfiles trgfiles)
+
+  set(_trgfiles)
+  foreach(srcfile IN LISTS srcfiles)
+    string(REGEX REPLACE "\\.${srcext}$" ".${trgext}" trgfile ${srcfile})
+    add_custom_command(
+      OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${trgfile}
+      COMMAND ${preproc} ${preprocopts} ${CMAKE_CURRENT_SOURCE_DIR}/${srcfile} ${CMAKE_CURRENT_BINARY_DIR}/${trgfile}
+      MAIN_DEPENDENCY ${CMAKE_CURRENT_SOURCE_DIR}/${srcfile})
+    list(APPEND _trgfiles ${CMAKE_CURRENT_BINARY_DIR}/${trgfile})
+  endforeach()
+  set(${trgfiles} ${_trgfiles} PARENT_SCOPE)
+
+endfunction()
+
+# Preprocesses fortran files with fypp.
+#
+# It assumes that source files have the ".fypp" extension. Target files will be
+# created with the extension ".f90". The FYPP variable must contain the path to
+# the fypp-preprocessor.
+#
+# Args:
+#     fyppopts [in]: Options to pass to fypp.
+#     fyppfiles [in]: Files to be processed by fypp
+#     f90files [out]: List of created f90 files on exit
+#
+function (fypp_f90 fyppopts fyppfiles f90files)
+  preprocess("${FYPP}" "${fyppopts}" "fypp" "f90" "${fyppfiles}" _f90files)
+  set(${f90files} ${_f90files} PARENT_SCOPE)
+endfunction()

src/CMakeLists.txt

@@ -5,12 +5,17 @@ find_package(LAPACK REQUIRED)
 
 configure_file(clima_version.f90.in ${CMAKE_CURRENT_BINARY_DIR}/clima_version.f90)
 
+set(fppFiles 
+  clima_saturationdata.f90
+)
+
+fypp_f90("" "${fppFiles}" outFiles)
+
 set(BASE_SOURCES
   ${CMAKE_CURRENT_BINARY_DIR}/clima_version.f90
   clima_const.f90
   clima_useful.f90
   clima_eqns.f90
-  clima_saturationdata.f90
   clima_eqns_water.f90
   clima_types.f90
   clima_types_create.f90
@@ -42,6 +47,7 @@ add_library(clima
   ${RADIATE_SOURCES}
   ${CLIMATE_SOURCES}
   ${ADIABAT_SOURCES}
+  ${outFiles}
   clima.f90
 )
 target_link_libraries(clima 
@@ -51,6 +57,7 @@ target_link_libraries(clima
   finterp
   dop853
   minpack
+  forwarddiff
   ${LAPACK_LIBRARIES}
 )
 

src/adiabat/clima_adiabat_rc.f90

@@ -468,7 +468,7 @@ subroutine root_fcn(d, P, gout)
     real(dp), intent(in) :: P
     real(dp), intent(out) :: gout(:) 
 
-    real(dp) :: T, f_dry, P_sat, dTdP
+    real(dp) :: T, f_dry, P_sat, dTdlog10P
     logical :: super_saturated
     integer :: i
 
@@ -477,7 +477,7 @@ subroutine root_fcn(d, P, gout)
     d%P_i_cur = d%f_i_cur*P
 
     if (any(d%sp_type == CondensingSpeciesType)) then
-      dTdP = dT_dP(d, P)
+      call d%T%evaluate_derivative(log10(P), dTdlog10P)
     endif
 
     do i = 1,d%sp%ng
@@ -487,10 +487,28 @@ subroutine root_fcn(d, P, gout)
         P_sat = d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)
       endif
 
-      if (d%sp_type(i) == CondensingSpeciesType) then
-        ! Search for cold trap (place where Temperature decreases with altitude
-        gout(i) = dTdP - 1.0e-8_dp
-      elseif (d%sp_type(i) == DrySpeciesType) then
+      if (d%sp_type(i) == CondensingSpeciesType) then; block
+        real(dp) :: dPi_dT, dTdP, dPi_dP, dfi_dP, dlog10fi_dP
+
+        ! Derivative of SVP curve for species i
+        dPi_dT = d%RH(i)*d%sp%g(i)%sat%sat_pressure_derivative(T)
+
+        ! Derivative of temperature profile
+        dTdP = dTdlog10P*(1.0_dp/(P*log(10.0_dp)))
+
+        ! Get dP_i/dP
+        dPi_dP = dPi_dT*dTdP
+
+        ! Convert to df_i/dP, where f_i is mixing ratio of condensible
+        dfi_dP = (1.0_dp/P)*(dPi_dP) - P_sat/P**2.0_dp
+
+        ! Convert to the derivative of the log of the mixing ratio
+        dlog10fi_dP = dfi_dP*(1/(d%f_i_cur(i)*log(10.0_dp)))
+
+        ! The root occurs when the mixing ratio begins to increase in concentration
+        gout(i) = dlog10fi_dP - 1.0e-8_dp
+
+      end block; elseif (d%sp_type(i) == DrySpeciesType) then
         ! Dry species can become condensing species if
         ! they reach saturation.
         gout(i) = d%P_i_cur(i)/P_sat - (1.0_dp + 1.0e-8_dp)
@@ -549,7 +567,6 @@ subroutine mixing_ratios(d, P, T, f_i_layer, f_dry, super_saturated)
     do i = 1,d%sp%ng
       if (d%sp_type(i) == CondensingSpeciesType) then
         f_i_layer(i) = min(d%RH(i)*d%sp%g(i)%sat%sat_pressure(T)/P,1.0_dp)
-        if (f_i_layer(i) == 1.0_dp) super_saturated = .true.
         f_moist = f_moist + f_i_layer(i)
       endif
     enddo
@@ -575,29 +592,6 @@ subroutine mixing_ratios(d, P, T, f_i_layer, f_dry, super_saturated)
 
   end subroutine
 
-  function dT_dP(d, P) result(dTdP)
-    type(AdiabatRCProfileData), intent(inout) :: d
-    real(dp), intent(in) :: P
-    real(dp) :: dTdP
-
-    real(dp) :: P2, P1, deltaP, log10P
-    real(dp) :: T2, T1
-
-    log10P = log10(P)
-    P2 = log10P + log10P*d%epsj
-    P2 = max(min(P2, log10(d%P_surf)),log10(d%P_top))
-    P2 = min(P2, log10(d%P_root))
-    P1 = log10P - log10P*d%epsj
-    P1 = max(min(P1, log10(d%P_surf)),log10(d%P_top))
-    deltaP = P2 - P1
-
-    call d%T%evaluate(P2, T2)
-    call d%T%evaluate(P1, T1)
-
-    dTdP = (T2 - T1)/deltaP
-
-  end function
-
   function general_adiabat_lapse_rate(d, P) result(dlnT_dlnP)
     use clima_eqns, only: heat_capacity_eval
     use clima_eqns_water, only: Rgas_cgs => Rgas

src/clima_saturationdata.f90

@@ -1,3 +1,7 @@
+#:set TYPES = ['real(dp)', 'type(dual)']
+#:set NAMES = ['real', 'dual']
+#:set TYPES_NAMES = list(zip(TYPES, NAMES))
+
 module clima_saturationdata
   use clima_const, only: dp, s_str_len
   implicit none
@@ -28,15 +32,22 @@ module clima_saturationdata
     procedure :: latent_heat_vap
     procedure :: latent_heat_sub
     procedure :: latent_heat
-    procedure :: sat_pressure_crit
-    procedure :: sat_pressure_vap
-    procedure :: sat_pressure_sub
-    procedure :: sat_pressure
+    procedure :: sat_pressure_crit => sat_pressure_crit_real
+    procedure :: sat_pressure_vap => sat_pressure_vap_real
+    procedure :: sat_pressure_sub => sat_pressure_sub_real
+    procedure :: sat_pressure => sat_pressure_real
+    procedure :: sat_pressure_derivative
   end type
   interface SaturationData
     procedure create_SaturationData
   end interface
 
+  #:for NAME in ['sat_pressure_vap','sat_pressure_sub','sat_pressure','integral_fcn','sat_pressure_crit']
+  interface ${NAME}$
+    module procedure :: ${NAME}$_real, ${NAME}$_dual
+  end interface
+  #:endfor
+
 contains
 
   !> Latent heat of condensation above the critical point.
@@ -79,63 +90,98 @@ function latent_heat(self, T) result(L)
     endif
   end function
 
+  #:for TYPE1, NAME in TYPES_NAMES
   !> Saturation pressure of a super-critical gas.
   !> This is non-physical. Its to approximate transitions in atmospheres
   !> between super-critical and sub-critical gases.
-  function sat_pressure_crit(self, T) result(p_sat)
+  function sat_pressure_crit_${NAME}$(self, T) result(p_sat)
+    #:if NAME == 'dual'
+    use forwarddiff
+    #:endif
     use clima_const, only: Rgas
     class(SaturationData), intent(inout) :: self
-    real(dp), intent(in) :: T !! K
-    real(dp) :: p_sat !! dynes/cm2
-    real(dp) :: tmp
+    ${TYPE1}$, intent(in) :: T !! K
+    ${TYPE1}$ :: p_sat !! dynes/cm2
+    ${TYPE1}$ :: tmp
     tmp = (integral_fcn(self%a_v, self%b_v, self%T_critical) - integral_fcn(self%a_v, self%b_v, self%T_ref)) + &
           (integral_fcn(self%a_c, self%b_c, T) - integral_fcn(self%a_c, self%b_c, self%T_critical))
     p_sat = self%P_ref*exp((self%mu/Rgas)*(tmp))
   end function
 
   !> Saturation pressure over liquid
-  function sat_pressure_vap(self, T) result(p_sat)
+  function sat_pressure_vap_${NAME}$(self, T) result(p_sat)
+    #:if NAME == 'dual'
+    use forwarddiff
+    #:endif
     use clima_const, only: Rgas
     class(SaturationData), intent(inout) :: self
-    real(dp), intent(in) :: T !! K
-    real(dp) :: p_sat !! dynes/cm2
-    real(dp) :: tmp
+    ${TYPE1}$, intent(in) :: T !! K
+    ${TYPE1}$ :: p_sat !! dynes/cm2
+    ${TYPE1}$ :: tmp
     tmp = integral_fcn(self%a_v, self%b_v, T) - integral_fcn(self%a_v, self%b_v, self%T_ref)
     p_sat = self%P_ref*exp((self%mu/Rgas)*(tmp))
   end function
 
   !> Saturation pressure over solid
-  function sat_pressure_sub(self, T) result(p_sat)
+  function sat_pressure_sub_${NAME}$(self, T) result(p_sat)
+    #:if NAME == 'dual'
+    use forwarddiff
+    #:endif
     use clima_const, only: Rgas
     class(SaturationData), intent(inout) :: self
-    real(dp), intent(in) :: T !! K
-    real(dp) :: p_sat !! dynes/cm2
-    real(dp) :: tmp
+    ${TYPE1}$, intent(in) :: T !! K
+    ${TYPE1}$ :: p_sat !! dynes/cm2
+    ${TYPE1}$ :: tmp
     tmp = (integral_fcn(self%a_v, self%b_v, self%T_triple) - integral_fcn(self%a_v, self%b_v, self%T_ref)) + &
           (integral_fcn(self%a_s, self%b_s, T) - integral_fcn(self%a_s, self%b_s, self%T_triple))
     p_sat = self%P_ref*exp((self%mu/Rgas)*(tmp))
   end function
 
   !> Saturation pressure over liquid or solid
-  function sat_pressure(self, T) result(p_sat)
+  function sat_pressure_${NAME}$(self, T) result(p_sat)
+    #:if NAME == 'dual'
+    use forwarddiff
+    #:endif
     class(SaturationData), intent(inout) :: self
-    real(dp), intent(in) :: T !! K
-    real(dp) :: p_sat !! dynes/cm2
+    ${TYPE1}$, intent(in) :: T !! K
+    ${TYPE1}$ :: p_sat !! dynes/cm2
     if (T >= self%T_critical) then
-      p_sat = self%sat_pressure_crit(T) ! non-physical
+      p_sat = sat_pressure_crit(self, T) ! non-physical
     elseif (T > self%T_triple .and. T < self%T_critical) then
-      p_sat = self%sat_pressure_vap(T)
+      p_sat = sat_pressure_vap(self, T)
     else ! (T <= T_triple) then
-      p_sat = self%sat_pressure_sub(T)
+      p_sat = sat_pressure_sub(self, T)
     endif
   end function
 
   !> This is $\int L/T^2 dT$
-  function integral_fcn(A, B, T) result(res)
-    real(dp), intent(in) :: A, B, T !! K
-    real(dp) :: res
+  function integral_fcn_${NAME}$(A, B, T) result(res)
+    #:if NAME == 'dual'
+    use forwarddiff
+    #:endif
+    real(dp), intent(in) :: A, B
+    ${TYPE1}$, intent(in) :: T !! K
+    ${TYPE1}$ :: res
     res = -A/T + B*log(T)
   end function
+  #:endfor
+
+  !> Compute the derivative of the SVP function.
+  function sat_pressure_derivative(self, T) result(dPdT)
+    use forwarddiff, only: derivative
+    class(SaturationData), intent(inout) :: self
+    real(dp), intent(in) :: T !! K
+    real(dp) :: dPdT !! dP/dT where P is dynes/cm^2 and T is K
+    real(dp) :: p_sat
+    call derivative(fcn, T, p_sat, dPdT)
+  contains
+    function fcn(x_) result(res_)
+      use forwarddiff, only: dual
+      type(dual), intent(in) :: x_
+      type(dual) :: res_
+      res_ = sat_pressure_dual(self, x_)
+    end function
+  end function
 
   function create_SaturationData(s, name, filename, err) result(sat)
     use fortran_yaml_c_types, only: type_dictionary, type_error
@@ -255,6 +301,9 @@ function create_SaturationData(s, name, filename, err) result(sat)
               'Computed saturation vapor pressures are negative or nan.'
         return
       endif
+
+      P(2) = sat%sat_pressure_derivative(sat%T_triple+1.0e-8_dp)
+
     endblock
 
   end function

src/dependencies/CMakeLists.txt

@@ -1,4 +1,14 @@
 
+CPMAddPackage(
+  NAME forwarddiff
+  VERSION 0.1.2
+  OPTIONS 
+    "BUILD_EXECUTABLE OFF"
+  GITHUB_REPOSITORY "Nicholaswogan/forwarddiff"
+  GIT_TAG "v0.1.2"
+  EXCLUDE_FROM_ALL ON
+)
+
 # fortran-yaml-c
 CPMAddPackage(
   NAME fortran-yaml-c
@@ -29,15 +39,18 @@ CPMAddPackage(
 )
 
 # finterp
-CPMAddPackage(
-  NAME finterp
-  VERSION 1.3.0
-  GITHUB_REPOSITORY "jacobwilliams/finterp"
-  GIT_TAG "1.3.0"
-  DOWNLOAD_ONLY ON
+# CPMAddPackage(
+#   NAME finterp
+#   VERSION 1.3.0
+#   GITHUB_REPOSITORY "jacobwilliams/finterp"
+#   GIT_TAG "1.3.0"
+#   DOWNLOAD_ONLY ON
+# )
+set(fppFiles 
+  linear_interpolation_module.F90
 )
-
-add_library(finterp ${finterp_SOURCE_DIR}/src/linear_interpolation_module.F90 )
+fypp_f90("" "${fppFiles}" outFiles)
+add_library(finterp ${outFiles})
 
 CPMAddPackage(
   NAME dop853