diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 25dcf1e9a..74df88872 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,9 @@ jobs: steps: - uses: actions/checkout@v4 - name: Install dependencies - run: sudo apt-get install -y --no-install-recommends gfortran libhdf5-openmpi-dev libopenmpi-dev libblas-dev liblapack-dev + run: | + sudo apt-get update + sudo apt-get install -y --no-install-recommends gfortran libhdf5-openmpi-dev libopenmpi-dev libblas-dev liblapack-dev - name: Install ED2 run: | cd ED/build @@ -43,7 +45,9 @@ jobs: with: name: ed2-binary - name: Install dependencies - run: sudo apt-get install -y --no-install-recommends gfortran libhdf5-openmpi-dev libopenmpi-dev libblas-dev liblapack-dev + run: | + sudo apt-get update + sudo apt-get install -y --no-install-recommends gfortran libhdf5-openmpi-dev libopenmpi-dev libblas-dev liblapack-dev - name: Run ED2 working-directory: ./EDTS/ run: | diff --git a/BRAMS/run/RAMSIN b/BRAMS/run/RAMSIN index e5749dfae..f020285be 100644 --- a/BRAMS/run/RAMSIN +++ b/BRAMS/run/RAMSIN @@ -1296,6 +1296,30 @@ $MODEL_OPTIONS + + + !------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a ! + ! fairly robust approach and it is continuous at the triple ! + ! point. The main downside is that it makes the estimate slightly ! + ! inconsistent with the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended ! + ! to become the default. It is accurate, and uses the same ! + ! assumptions as the definition of enthalpy and latent heat of ! + ! vaporisation and sublimation, so it is consistent with the ! + ! second law of thermodynamics. The only minor disadvantage is ! + ! that the function is continuous but not differentiable at the ! + ! fusion point. ! + !------------------------------------------------------------------------------------! + ESAT_METHOD = 1, + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! ! IGRNDVAP -- Methods to find the ground -> canopy conductance. ! ! 0. (ED-2.2 default) Modified Lee and Pielke (1992, J. Appl. Meteorol.), adding ! diff --git a/BRAMS/src/ed2/edcp_load_namelist.f90 b/BRAMS/src/ed2/edcp_load_namelist.f90 index 57e848909..76bdf4c60 100644 --- a/BRAMS/src/ed2/edcp_load_namelist.f90 +++ b/BRAMS/src/ed2/edcp_load_namelist.f90 @@ -158,6 +158,14 @@ subroutine read_ednl(iunit,filename) use canopy_layer_coms , only : crown_mod ! ! intent(out) use canopy_radiation_coms, only : icanrad & ! intent(out) , ihrzrad ! ! intent(out) + use consts_coms , only : day_sec ! ! intent(in) + use detailed_coms , only : dt_census & ! intent(out) + , yr1st_census & ! intent(out) + , mon1st_census & ! intent(out) + , min_recruit_dbh & ! intent(out) + , idetailed & ! intent(out) + , patch_keep ! ! intent(out) + use fusion_fission_coms , only : ifusion ! ! intent(out) !----- Coupled ED-BRAMS modules. -------------------------------------------------------! use mem_edcp , only : co2_offset ! ! intent(out) !----- BRAMS modules. ------------------------------------------------------------------! @@ -194,6 +202,7 @@ subroutine read_ednl(iunit,filename) , itimeh & ! intent(in) , imonthh & ! intent(in) , isoilflg ! ! intent(in) + use mem_basic , only : brams_esat_method => esat_method ! ! intent(in) use mem_leaf , only : nslcon & ! intent(in) , isoilcol & ! intent(in) , slz & ! intent(in) @@ -207,7 +216,7 @@ subroutine read_ednl(iunit,filename) , leaf_sldrain => sldrain & ! intent(in) , leaf_ipercol => ipercol & ! intent(in) , leaf_runoff_time => runoff_time ! ! intent(in) - use leaf_coms , only : leaf_ubmin => ubmin & ! intent(in) + , leaf_ubmin => ubmin & ! intent(in) , leaf_ugbmin => ugbmin & ! intent(in) , leaf_ustmin => ustmin & ! intent(in) , leaf_gamm => gamm & ! intent(in) @@ -217,14 +226,6 @@ subroutine read_ednl(iunit,filename) , leaf_leaf_maxwhc => leaf_maxwhc & ! intent(in) , leaf_min_patch_area => min_patch_area ! ! intent(in) use mem_radiate , only : radfrq ! ! intent(in) - use consts_coms , only : day_sec ! ! intent(in) - use detailed_coms , only : dt_census & ! intent(out) - , yr1st_census & ! intent(out) - , mon1st_census & ! intent(out) - , min_recruit_dbh & ! intent(out) - , idetailed & ! intent(out) - , patch_keep ! ! intent(out) - use fusion_fission_coms , only : ifusion ! ! intent(out) implicit none !----- Arguments. ----------------------------------------------------------------------! integer , intent(in) :: iunit ! Namelist unit number @@ -406,7 +407,7 @@ subroutine read_ednl(iunit,filename) ,nzg,nzs,isoilflg,nslcon,isoilcol,slz,slmstr,stgoff,ngrids & ,leaf_ubmin,leaf_ugbmin,leaf_ustmin,leaf_isoilbc,leaf_sldrain & ,leaf_ipercol,leaf_runoff_time,leaf_gamm,leaf_gamh,leaf_tprandtl & - ,leaf_ribmax,leaf_leaf_maxwhc) + ,leaf_ribmax,leaf_leaf_maxwhc,brams_esat_method) !---------------------------------------------------------------------------------------! ! The following variables can be defined in the regular ED2IN file for stand-alone ! @@ -592,7 +593,7 @@ subroutine copy_in_bramsnl(expnme_b,runtype_b,itimez_b,idatez_b,imonthz_b,iyearz ,iclobber_b,nzg_b,nzs_b,isoilflg_b,nslcon_b,isoilcol_b,slz_b & ,slmstr_b,stgoff_b,ngrids_b,ubmin_b,ugbmin_b,ustmin_b,isoilbc_b & ,sldrain_b,ipercol_b,runoff_time_b,gamm_b,gamh_b,tprandtl_b & - ,ribmax_b,leaf_maxwhc_b) + ,ribmax_b,leaf_maxwhc_b,esat_method_b) use ed_misc_coms , only : expnme & ! intent(out) , runtype & ! intent(out) , itimez & ! intent(out) @@ -608,7 +609,8 @@ subroutine copy_in_bramsnl(expnme_b,runtype_b,itimez_b,idatez_b,imonthz_b,iyearz , imonthh & ! intent(out) , iyearh & ! intent(out) , iclobber & ! intent(out) - , radfrq ! ! intent(out) + , radfrq & ! intent(out) + , esat_method ! ! intent(out) use grid_coms , only : centlon & ! intent(out) , centlat & ! intent(out) , deltax & ! intent(out) @@ -698,6 +700,7 @@ subroutine copy_in_bramsnl(expnme_b,runtype_b,itimez_b,idatez_b,imonthz_b,iyearz real , intent(in) :: gamh_b ! Sfc. lyr. Gamma for heat real , intent(in) :: tprandtl_b ! Turbulent Prandtl number real , intent(in) :: leaf_maxwhc_b ! Leaf max. water holding cap. + integer , intent(in) :: esat_method_b ! Method for esat. !---------------------------------------------------------------------------------------! @@ -760,6 +763,8 @@ subroutine copy_in_bramsnl(expnme_b,runtype_b,itimez_b,idatez_b,imonthz_b,iyearz ribmax = ribmax_b leaf_maxwhc = leaf_maxwhc_b isoilcol = isoilcol_b + + esat_method = esat_method_b !---------------------------------------------------------------------------------------! return diff --git a/BRAMS/src/ed2/edcp_mpiutils.F90 b/BRAMS/src/ed2/edcp_mpiutils.F90 index 3b3445b65..4362473f9 100644 --- a/BRAMS/src/ed2/edcp_mpiutils.F90 +++ b/BRAMS/src/ed2/edcp_mpiutils.F90 @@ -42,6 +42,7 @@ subroutine masterput_ednl(mainnum) , thsums_database & ! intent(in) , ivegt_dynamics & ! intent(in) , ibigleaf & ! intent(in) + , esat_method & ! intent(in) , integration_scheme & ! intent(in) , nsub_euler & ! intent(in) , end_time & ! intent(in) @@ -263,6 +264,7 @@ subroutine masterput_ednl(mainnum) call MPI_Bcast(sldrain,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(ivegt_dynamics,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(ibigleaf,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(esat_method,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(integration_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(nsub_euler,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(rk4_tolerance,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) @@ -474,6 +476,7 @@ subroutine nodeget_ednl(master_num) , thsums_database & ! intent(out) , ivegt_dynamics & ! intent(out) , ibigleaf & ! intent(out) + , esat_method & ! intent(out) , integration_scheme & ! intent(out) , nsub_euler & ! intent(out) , end_time & ! intent(out) @@ -695,7 +698,8 @@ subroutine nodeget_ednl(master_num) call MPI_Bcast(isoilbc,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(sldrain,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(ivegt_dynamics,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) - call MPI_Bcast(ibigleaf,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(ibigleaf,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(esat_method,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(integration_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(nsub_euler,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(rk4_tolerance,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) diff --git a/BRAMS/src/init/rdint.f90 b/BRAMS/src/init/rdint.f90 index 056104641..a3fc12fa5 100644 --- a/BRAMS/src/init/rdint.f90 +++ b/BRAMS/src/init/rdint.f90 @@ -1063,7 +1063,8 @@ subroutine read_nl(filename) , maxubtp ! ! intent(out) use mem_basic , only : co2_on & ! intent(out) , co2con & ! intent(out) - , ico2 ! ! intent(out) + , ico2 & ! intent(out) + , esat_method ! ! intent(out) use domain_decomp , only : domain_fname ! ! intent(out) use therm_lib , only : vapour_on & ! intent(out) , cloud_on & ! intent(out) @@ -1141,15 +1142,15 @@ subroutine read_nl(filename) namelist /MODEL_OPTIONS/ naddsc,icorflg,iadvec,iexev,imassflx,ibnd,jbnd,cphas & ,lsflg,nfpt,distim,iswrtyp,ilwrtyp,icumfdbk,radfrq,lonrad & ,npatch,nvegpat,min_patch_area,isfcl,dtleaf,ndtveg,istar & - ,igrndvap,ubmin,ugbmin,ustmin,gamm,gamh,tprandtl,ribmax & - ,leaf_maxwhc,ico2,co2con,nvgcon,pctlcon,nslcon,isoilcol & - ,drtcon,zrough,albedo,seatmp,dthcon,soil_moist & - ,soil_moist_fail,usdata_in,usmodel_in,slz,slmstr,stgoff & - ,isoilbc,sldrain,ipercol,runoff_time,if_urban_canopy & - ,idiffk,ibruvais,ibotflx,ihorgrad,csx,csz,xkhkm,zkhkm & - ,nna,nnb,nnc,akmin,akmax,hgtmin,hgtmax,level,icloud,irain & - ,ipris,isnow,iaggr,igraup,ihail,cparm,rparm,pparm,sparm & - ,aparm,gparm,hparm,gnu + ,esat_method,igrndvap,ubmin,ugbmin,ustmin,gamm,gamh & + ,tprandtl,ribmax,leaf_maxwhc,ico2,co2con,nvgcon,pctlcon & + ,nslcon,isoilcol,drtcon,zrough,albedo,seatmp,dthcon & + ,soil_moist,soil_moist_fail,usdata_in,usmodel_in,slz & + ,slmstr,stgoff,isoilbc,sldrain,ipercol,runoff_time & + ,if_urban_canopy,idiffk,ibruvais,ibotflx,ihorgrad,csx,csz & + ,xkhkm,zkhkm,nna,nnb,nnc,akmin,akmax,hgtmin,hgtmax,level & + ,icloud,irain,ipris,isnow,iaggr,igraup,ihail,cparm,rparm & + ,pparm,sparm,aparm,gparm,hparm,gnu namelist /MODEL_SOUND/ ipsflg,itsflg,irtsflg,iusflg,hs,ps,ts,rts,us,vs,co2s @@ -1811,6 +1812,7 @@ subroutine read_nl(filename) write (unit=*,fmt=*) ' dtleaf =',dtleaf write (unit=*,fmt=*) ' ndtveg =',ndtveg write (unit=*,fmt=*) ' istar =',istar + write (unit=*,fmt=*) ' esat_method =',esat_method write (unit=*,fmt=*) ' igrndvap =',igrndvap write (unit=*,fmt=*) ' ubmin =',ubmin write (unit=*,fmt=*) ' ugbmin =',ugbmin diff --git a/BRAMS/src/io/opspec.f90 b/BRAMS/src/io/opspec.f90 index e4994e6a7..e811150e8 100644 --- a/BRAMS/src/io/opspec.f90 +++ b/BRAMS/src/io/opspec.f90 @@ -1136,6 +1136,13 @@ subroutine opspec3 ifaterr=ifaterr+1 end if + ! Check whether the choice of vapour pressure formula is valid. + if (esat_method < 0 .or. esat_method > 2) then + print *, 'fatal - ESAT_METHOD must be between 0 and 2, and yours is set to ' & + ,esat_method,'...' + ifaterr=ifaterr+1 + end if + ! Check whether the ground vapour method scheme the user chose is fine or not. if (isfcl /= 0 .and. (igrndvap < 0 .or. igrndvap > 5)) then print *, 'fatal - IGRNDVAP must be between 0 and 5, and yours is set to ' & diff --git a/BRAMS/src/lib/rconstants.f90 b/BRAMS/src/lib/rconstants.f90 index 6a9b8d193..d0b6b6a53 100644 --- a/BRAMS/src/lib/rconstants.f90 +++ b/BRAMS/src/lib/rconstants.f90 @@ -35,7 +35,7 @@ Module rconstants !---------------------------------------------------------------------------------------! real, parameter :: srtwo = 1.414213562373095 ! Square root of 2. [ ---] real, parameter :: srthree = 1.732050807568877 ! Square root of 3. [ ---] - real, parameter :: sqrt2o2 = 0.5 * srtwo ! ½ Square root of 2. [ ---] + real, parameter :: sqrt2o2 = 0.5 * srtwo ! 1/2 Square root of 2. [ ---] real, parameter :: srtwoi = 1./srtwo ! 1./ Square root of 2. [ ---] real, parameter :: srthreei = 1./srthree ! 1./ Square root of 3. [ ---] real, parameter :: onethird = 1./3. ! 1/3 [ ---] @@ -48,15 +48,15 @@ Module rconstants !---------------------------------------------------------------------------------------! ! Universal constants ! !---------------------------------------------------------------------------------------! - real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m²/K^4] - real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m²kg/s²/K] + real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m2/K4] + real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m2kg/s2/K] real, parameter :: avogrado = 6.02252e+23 ! Avogrado number [ #/mol] real, parameter :: loschmidt = 2.68719e+25 ! Loschmidt number [ #/m3] real, parameter :: loschcgs = loschmidt*1.e-6 ! Loschmidt number [ #/cm3] - real, parameter :: t00 = 273.15 ! 0°C [ °C] + real, parameter :: t00 = 273.15 ! 0 degC [ K] real, parameter :: rmol = 8.314510 ! Molar gas constant [ J/mol/K] real, parameter :: rmolcgs = rmol*1.e7 ! Molar gas constant [erg/mol/K] - real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m³] + real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m3] real, parameter :: volmoll = volmol*1.e3 ! Molar volume at STP [ L] !---------------------------------------------------------------------------------------! @@ -102,9 +102,9 @@ Module rconstants ! General Earth properties ! !---------------------------------------------------------------------------------------! real, parameter :: vonk = 0.40 ! Von Kármán constant [ ---] - real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s²] - real, parameter :: gcgs = grav * 100. ! Gravity acceleration [ cm/s²] - real, parameter :: gg = .5 * grav ! ½ g [ m/s²] + real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s2] + real, parameter :: gcgs = grav * 100. ! Gravity acceleration [ cm/s2] + real, parameter :: gg = .5 * grav ! ½ g [ m/s2] real, parameter :: erad = 6370997. ! Earth radius [ m] real, parameter :: spcon = pio180*erad ! One degree of latitude [ m] real, parameter :: spconkm = spcon*0.001 ! One degree of latitude [ km] @@ -113,7 +113,7 @@ Module rconstants real, parameter :: ss60 = 1.8663 ! Polar stereo conversion to 60° [ ] real, parameter :: omega = 7.292e-5 ! Earth's rotation speed [ rad/s] real, parameter :: viscos = .15e-4 ! Viscosity coefficient [ ] - real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m²] + real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m2] real, parameter :: p00 = 1.e5 ! Reference pressure [ Pa] real, parameter :: prefsea = 101325. ! Reference sea level pressure [ Pa] real, parameter :: p00i = 1. / p00 ! 1/p00 [ 1/Pa] @@ -127,7 +127,7 @@ Module rconstants ! third edition, Academic Press, Amsterdam, 418pp. (Chapters 3 and 10). ! ! ! ! Air diffusion properties. These properties are temperature-dependent in reality, ! - ! but for simplicity we assume them constants, using the value at 20°C. ! + ! but for simplicity we assume them constants, using the value at 20 degC. ! ! ! ! Thermal diffusivity - Computed from equation on page 32 of MU08; ! ! Kinematic viscosity - Computed from equation on page 32 of MU08; ! @@ -135,9 +135,9 @@ Module rconstants ! 10.11 (MU08). ! ! These terms could be easily made function of temperature in the future if needed be. ! !---------------------------------------------------------------------------------------! - real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m²/s] + real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m2/s] real, parameter :: dth_diff = 0.007 ! Temperature dependency slope [ 1/K] - real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m²/s] + real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m2/s] real, parameter :: dkin_visc = 0.007 ! Temperature dependency slope [ 1/K] !---------------------------------------------------------------------------------------! @@ -187,8 +187,8 @@ Module rconstants !---------------------------------------------------------------------------------------! ! Liquid water properties ! !---------------------------------------------------------------------------------------! - real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m³] - real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m³/kg] + real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m3] + real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cliq = 4.186e3 ! Liquid water specific heat (Cl) [ J/kg/K] real, parameter :: cliqi = 1./cliq ! Inverse of water heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! @@ -198,12 +198,12 @@ Module rconstants !---------------------------------------------------------------------------------------! ! Ice properties ! !---------------------------------------------------------------------------------------! - real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m³] - real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m³/kg] - real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m³] - real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m³/kg] - real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m³] - real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m³/kg] + real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m3] + real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m3/kg] + real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m3] + real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m3/kg] + real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m3] + real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cice = 2.093e3 ! Ice specific heat (Ci) [ J/kg/K] real, parameter :: cicei = 1. / cice ! Inverse of ice heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! @@ -237,6 +237,21 @@ Module rconstants !---------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! + ! These derived variables are used for computing vapour pressure at saturation when ! + ! using the method proposed by Ambaum (2020). ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure over ! + ! water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899. ! + !---------------------------------------------------------------------------------------! + real, parameter :: dcplorv = (cliq - cph2o) / rh2o ! (cl - cp) / Rv [ ------] + real, parameter :: dcpiorv = (cice - cph2o) / rh2o ! (ci - cp) / Rv [ ------] + real, parameter :: lvorvot3 = alvl3 / (rh2o * t3ple) ! Lv0 / ( Rv * T3) [ ------] + real, parameter :: liorvot3 = alvi3 / (rh2o * t3ple) ! Ls0 / ( Rv * T3) [ ------] + !---------------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------------! ! The following variables are useful when defining the derivatives of theta_il. ! ! They correspond to L?(T) - L?' T. ! @@ -356,11 +371,11 @@ Module rconstants real, parameter :: mol_2_umol = 1.e6 ! mol => µmol real, parameter :: umol_2_mol = 1.e-6 ! µmol => mol real, parameter :: umol_2_kgC = 1.20107e-8 ! µmol(CO2) => kg(C) - real, parameter :: Watts_2_Ein = 4.6e-6 ! W/m2 => mol/m²/s - real, parameter :: Ein_2_Watts = 1./Watts_2_Ein ! mol/m²/s => W/m2 + real, parameter :: Watts_2_Ein = 4.6e-6 ! W/m2 => mol/m2/s + real, parameter :: Ein_2_Watts = 1./Watts_2_Ein ! mol/m2/s => W/m2 real, parameter :: kgC_2_umol = 1. / umol_2_kgC ! kg(C) => µmol(CO2) - real, parameter :: kgom2_2_tonoha = 10. ! kg(C)/m² => ton(C)/ha - real, parameter :: tonoha_2_kgom2 = 0.1 ! ton(C)/ha => kg(C)/m² + real, parameter :: kgom2_2_tonoha = 10. ! kg(C)/m2 => ton(C)/ha + real, parameter :: tonoha_2_kgom2 = 0.1 ! ton(C)/ha => kg(C)/m2 real, parameter :: umols_2_kgCyr = umol_2_kgC * yr_sec ! µmol(CO2)/s => kg(C)/yr real, parameter :: kgCday_2_umols = kgC_2_umol / day_sec ! kg(C)/day => µmol(CO2)/s !---------------------------------------------------------------------------------------! @@ -463,6 +478,10 @@ Module rconstants real(kind=8), parameter :: uiliqt38 = dble(uiliqt3 ) real(kind=8), parameter :: dcpvl8 = dble(dcpvl ) real(kind=8), parameter :: dcpvi8 = dble(dcpvi ) + real(kind=8), parameter :: dcplorv8 = dble(dcplorv ) + real(kind=8), parameter :: dcpiorv8 = dble(dcpiorv ) + real(kind=8), parameter :: lvorvot38 = dble(lvorvot3 ) + real(kind=8), parameter :: liorvot38 = dble(liorvot3 ) real(kind=8), parameter :: del_alvl38 = dble(del_alvl3 ) real(kind=8), parameter :: del_alvi38 = dble(del_alvi3 ) real(kind=8), parameter :: tsupercool_liq8 = dble(tsupercool_liq) diff --git a/BRAMS/src/lib/therm_lib.f90 b/BRAMS/src/lib/therm_lib.f90 index 9c7624f01..40cd06af3 100644 --- a/BRAMS/src/lib/therm_lib.f90 +++ b/BRAMS/src/lib/therm_lib.f90 @@ -27,9 +27,6 @@ module therm_lib integer, parameter :: maxlev = 16 ! Maximum # of levels for adaptive ! quadrature methods. - - logical, parameter :: newthermo = .true. ! Use new thermodynamics [T|F] - !---------------------------------------------------------------------------------------! ! This is the "level" variable, that used to be in micphys. Since it affects more the ! ! thermodynamics choices than the microphysics, it was moved to here. ! @@ -131,7 +128,15 @@ module therm_lib ! Kelvin temperature. This expression came from MK05, equation (10). ! !---------------------------------------------------------------------------------------! real(kind=4) function eslf(temp,l1funout,l2funout,ttfunout) - use rconstants , only : t00 ! ! intent(in) + use rconstants , only : t00 & ! intent(in) + , es3ple & ! intent(in) + , t3ple & ! intent(in) + , dcplorv & ! intent(in) + , lvorvot3 & ! intent(in) + , rh2o & ! intent(in) + , lnexp_min & ! intent(in) + , lnexp_max ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=4), intent(in) :: temp ! Temperature [ K] @@ -149,30 +154,59 @@ real(kind=4) function eslf(temp,l1funout,l2funout,ttfunout) !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-80.,temp-t00) + eslf = cll(0) + x * (cll(1) + x * (cll(2) + x * (cll(3) + x * (cll(4) & + + x * (cll(5) + x * (cll(6) + x * (cll(7) + x * cll(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! + if (present(l1funout)) l1funout = eslf + if (present(l2funout)) l2funout = eslf + if (present(ttfunout)) ttfunout = eslf + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! l1fun = l01_10(0) + l01_10(1)/temp + l01_10(2)*log(temp) + l01_10(3) * temp l2fun = l02_10(0) + l02_10(1)/temp + l02_10(2)*log(temp) + l02_10(3) * temp ttfun = tanh(ttt_10(1) * (temp - ttt_10(2))) eslf = exp(l1fun + ttfun*l2fun) !---------------------------------------------------------------------------------! + !----- These are useful for the derivative calculation. --------------------------! if (present(l1funout)) l1funout = l1fun if (present(l2funout)) l2funout = l2fun if (present(ttfunout)) ttfunout = ttfun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-80.,temp-t00) - eslf = cll(0) + x * (cll(1) + x * (cll(2) + x * (cll(3) + x * (cll(4) & - + x * (cll(5) + x * (cll(6) + x * (cll(7) + x * cll(8)) ) ) ) ) ) ) !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + l1fun = min(lnexp_max,max(lnexp_min, lvorvot3 - alvl(temp) / (rh2o * temp))) + eslf = es3ple * (t3ple/temp) ** dcplorv * exp(l1fun) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(l1funout)) l1funout = eslf if (present(l2funout)) l2funout = eslf if (present(ttfunout)) ttfunout = eslf - end if + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! return end function eslf @@ -190,7 +224,15 @@ end function eslf ! Kelvin temperature, based on MK05 equation (7). ! !---------------------------------------------------------------------------------------! real(kind=4) function esif(temp,iifunout) - use rconstants , only : t00 ! ! intent(in) + use rconstants , only : t00 & ! intent(in) + , es3ple & ! intent(in) + , t3ple & ! intent(in) + , dcpiorv & ! intent(in) + , lsorvot3 & ! intent(in) + , rh2o & ! intent(in) + , lnexp_min & ! intent(in) + , lnexp_max ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=4), intent(in) :: temp ! Temperature [ K] @@ -202,27 +244,53 @@ real(kind=4) function esif(temp,iifunout) !------------------------------------------------------------------------------------! + !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x=max(-80.,temp-t00) + esif = cii(0) + x * (cii(1) + x * (cii(2) + x * (cii(3) + x * (cii(4) & + + x * (cii(5) + x * (cii(6) + x * (cii(7) + x * cii(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! - !----- Updated method, using MK05 ------------------------------------------------! + !----- Dummy values for optional function in case they are requested. ------------! + if (present(iifunout)) iifunout=esif + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! iifun = iii_7(0) + iii_7(1)/temp + iii_7(2) * log(temp) + iii_7(3) * temp esif = exp(iifun) !---------------------------------------------------------------------------------! + !----- This is useful for the derivative calculation. ------==--------------------! if (present(iifunout)) iifunout=iifun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x=max(-80.,temp-t00) - esif = cii(0) + x * (cii(1) + x * (cii(2) + x * (cii(3) + x * (cii(4) & - + x * (cii(5) + x * (cii(6) + x * (cii(7) + x * cii(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + iifun = min(lnexp_max,max(lnexp_min, lsorvot3 - alvi(temp) / (rh2o * temp))) + esif = es3ple * (t3ple/temp) ** dcpiorv * exp(iifun) !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(iifunout)) iifunout=esif - end if + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -708,7 +776,9 @@ end function rhovsil ! pressure with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function eslfp(temp) - use rconstants , only : t00 ! ! intent(in) + use rconstants , only : t00 & ! intent(in) + , rh2o ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=4), intent(in) :: temp @@ -725,21 +795,37 @@ real(kind=4) function eslfp(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-80.,temp-t00) + eslfp = dll(0) + x * (dll(1) + x * (dll(2) + x * (dll(3) + x * (dll(4) & + + x * (dll(5) + x * (dll(6) + x * (dll(7) + x * dll(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! esl = eslf(temp,l2funout=l2fun,ttfunout=ttfun) l1prime = -l01_10(1)/(temp*temp) + l01_10(2)/temp + l01_10(3) l2prime = -l02_10(1)/(temp*temp) + l02_10(2)/temp + l02_10(3) ttprime = ttt_10(1)*(1.-ttfun*ttfun) eslfp = esl * (l1prime + l2prime*ttfun + l2fun*ttprime) - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-80.,temp-t00) - eslfp = dll(0) + x * (dll(1) + x * (dll(2) + x * (dll(3) + x * (dll(4) & - + x * (dll(5) + x * (dll(6) + x * (dll(7) + x * dll(8)) ) ) ) ) ) ) - end if + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esl = eslf(temp) + eslfp = esl * alvl(temp) / (rh2o * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! @@ -759,7 +845,9 @@ end function eslfp ! with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function esifp(temp) - use rconstants , only : t00 ! ! intent(in) + use rconstants , only : t00 & ! intent(in) + , rh2o ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=4), intent(in) :: temp @@ -772,21 +860,38 @@ real(kind=4) function esifp(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! - esi = esif(temp) - iiprime = -iii_7(1)/(temp*temp) + iii_7(2)/temp + iii_7(3) - esifp = esi * iiprime - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! x=max(-80.,temp-t00) esifp = dii(0) + x * (dii(1) + x * (dii(2) + x * (dii(3) + x * (dii(4) & + x * (dii(5) + x * (dii(6) + x * (dii(7) + x * dii(8)) ) ) ) ) ) ) - end if + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! + esi = esif(temp) + iiprime = -iii_7(1)/(temp*temp) + iii_7(2)/temp + iii_7(3) + esifp = esi * iiprime + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esi = esif(temp) + esifp = esi * alvi(temp) / (rh2o * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! + return end function esifp !=======================================================================================! @@ -2249,7 +2354,7 @@ real(kind=4) function vpdefil(pres,temp,humi,is_shv,useice) if (frozen) then psat = esif(temp) else - psat = esif(temp) + psat = eslf(temp) end if !------------------------------------------------------------------------------------! @@ -2722,9 +2827,9 @@ end function tl2uint !=======================================================================================! !=======================================================================================! - ! This function computes the extensive internal energy of water [J/m²] or [ J/m³], ! - ! given the temperature [K], the heat capacity of the "dry" part [J/m²/K] or [J/m³/K], ! - ! water mass [ kg/m²] or [ kg/m³], and liquid fraction [---]. ! + ! This function computes the extensive internal energy of water [J/m2] or [ J/m3], ! + ! given the temperature [K], the heat capacity of the "dry" part [J/m2/K] or [J/m3/K], ! + ! water mass [ kg/m2] or [ kg/m3], and liquid fraction [---]. ! !---------------------------------------------------------------------------------------! real(kind=4) function cmtl2uext(dryhcap,wmass,temp,fliq) use rconstants , only : cice & ! intent(in) @@ -2733,8 +2838,8 @@ real(kind=4) function cmtl2uext(dryhcap,wmass,temp,fliq) implicit none !----- Arguments. -------------------------------------------------------------------! - real(kind=4), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m²/K] or [J/m³/K] - real(kind=4), intent(in) :: wmass ! Water mass [ kg/m²] or [ kg/m³] + real(kind=4), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m2/K] or [J/m3/K] + real(kind=4), intent(in) :: wmass ! Water mass [ kg/m2] or [ kg/m3] real(kind=4), intent(in) :: temp ! Temperature [ K] real(kind=4), intent(in) :: fliq ! Liquid fraction (0-1) [ ---] !------------------------------------------------------------------------------------! @@ -2940,10 +3045,7 @@ end function alvi real(kind=4) function theta_iceliq(exner,temp,rliq,rice) use rconstants , only : alvl3 & ! intent(in) , alvi3 & ! intent(in) - , cpdry & ! intent(in) - , ttripoli & ! intent(in) - , htripoli & ! intent(in) - , htripolii ! ! intent(in) + , cpdry ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -2962,36 +3064,13 @@ real(kind=4) function theta_iceliq(exner,temp,rliq,rice) !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use the ! - ! latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - !------------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - else - qq = alvl3 * rliq + alvi3 * rice - end if + !----- Find the latent heat enthalpy. -----------------------------------------------! + qq = alvl(temp) * rliq + alvi(temp) * rice !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Solve the thermodynamics. For the new thermodynamics we don't approximate ! - ! the exponential to a linear function, nor do we impose temperature above the thre- ! - ! shold from Tripoli and Cotton (1981). ! - !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Decide how to compute, based on temperature. ------------------------------! - theta_iceliq = hh * exp(-qq / hh) / exner - !---------------------------------------------------------------------------------! - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli) then - theta_iceliq = hh * hh / (exner * ( hh + qq)) - else - theta_iceliq = hh * htripoli / (exner * ( htripoli + qq)) - end if - !---------------------------------------------------------------------------------! - end if + !----- Find theta_iceliq. -----------------------------------------------------------! + theta_iceliq = hh * exp(-qq / hh) / exner !------------------------------------------------------------------------------------! return @@ -3015,9 +3094,6 @@ real(kind=4) function dthetail_dt(condconst,thil,exner,pres,temp,rliq,ricein) , dcpvi & ! intent(in) , dcpvl & ! intent(in) , cpdry & ! intent(in) - , ttripoli & ! intent(in) - , htripoli & ! intent(in) - , htripolii & ! intent(in) , t3ple ! ! intent(in) implicit none @@ -3031,8 +3107,8 @@ real(kind=4) function dthetail_dt(condconst,thil,exner,pres,temp,rliq,ricein) real(kind=4), intent(in), optional :: ricein ! Ice mixing ratio [ kg/kg] !----- Local variables --------------------------------------------------------------! real(kind=4) :: rice ! Ice mixing ratio or 0. [ kg/kg] - real(kind=4) :: ldrst ! L × d(rs)/dT × T [ J/kg] - real(kind=4) :: rdlt ! r × d(L)/dT × T [ J/kg] + real(kind=4) :: ldrst ! L * d(rs)/dT * T [ J/kg] + real(kind=4) :: rdlt ! r * d(L)/dT * T [ J/kg] real(kind=4) :: hh ! Sensible heat enthalpy [ J/kg] real(kind=4) :: qq ! Latent heat enthalpy [ J/kg] logical :: thereisice ! Is ice present [ ---] @@ -3070,43 +3146,25 @@ real(kind=4) function dthetail_dt(condconst,thil,exner,pres,temp,rliq,ricein) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use ! - ! the latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - ! The term r × d(L)/dT × T is computed only when we use the new thermodynamics. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - rdlt = (dcpvl * rliq + dcpvi * rice ) * temp - else - qq = alvl3 * rliq + alvi3 * rice - rdlt = 0.0 - end if + !----- Find the latent heat enthalpy. --------------------------------------------! + qq = alvl(temp) * rliq + alvi(temp) * rice + rdlt = (dcpvl * rliq + dcpvi * rice ) * temp !---------------------------------------------------------------------------------! !---------------------------------------------------------------------------------! - ! This is the term L×[d(rs)/dt]×T. L may be either the vapourisation or ! + ! This is the term Lx[d(rs)/dt]xT. L may be either the vapourisation or ! ! sublimation latent heat, depending on the temperature and whether we are consi- ! - ! dering ice or not. We still need to check whether latent heat is a function of ! - ! temperature or not. Also, if condensation mixing ratio is constant, then this ! - ! term will be always zero. ! + ! dering ice or not. If condensation mixing ratio is constant, then this term ! + ! will be always zero. ! !---------------------------------------------------------------------------------! if (condconst) then ldrst = 0. elseif (thereisice .and. temp < t3ple) then - if (newthermo) then - ldrst = alvi3 * rsifp(pres,temp) * temp - else - ldrst = alvi(temp) * rsifp(pres,temp) * temp - end if + ldrst = alvi(temp) * rsifp(pres,temp) * temp else - if (newthermo) then - ldrst = alvl3 * rslfp(pres,temp) * temp - else - ldrst = alvl(temp) * rslfp(pres,temp) * temp - end if + ldrst = alvl(temp) * rslfp(pres,temp) * temp end if !---------------------------------------------------------------------------------! end if @@ -3117,17 +3175,7 @@ real(kind=4) function dthetail_dt(condconst,thil,exner,pres,temp,rliq,ricein) !------------------------------------------------------------------------------------! ! Find the condensed phase consistent with the thermodynamics used. ! !------------------------------------------------------------------------------------! - if (newthermo) then - dthetail_dt = thil * ( 1. + (ldrst + qq - rdlt ) / hh ) / temp - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli) then - dthetail_dt = thil * ( 1. + (ldrst + qq) / (hh+qq) ) / temp - else - dthetail_dt = thil * ( 1. + ldrst / (htripoli + alvl3 * rliq) ) / temp - end if - !---------------------------------------------------------------------------------! - end if + dthetail_dt = thil * ( 1. + (ldrst + qq - rdlt ) / hh ) / temp !------------------------------------------------------------------------------------! return @@ -3144,9 +3192,6 @@ end function dthetail_dt !=======================================================================================! ! This function computes temperature from the ice-liquid water potential temperature ! ! in Kelvin, Exner function in J/kg/K, and liquid and ice mixing ratios in kg/kg. ! - ! For now t1stguess is used only to decide whether I should use the complete case or ! - ! the 253 K to reduce the error on neglecting the changes on latent heat due to temper- ! - ! ature. ! !---------------------------------------------------------------------------------------! real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) use rconstants , only : cpdry & ! intent(in) @@ -3155,9 +3200,7 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) , alvl3 & ! intent(in) , alvi3 & ! intent(in) , t00 & ! intent(in) - , t3ple & ! intent(in) - , ttripoli & ! intent(in) - , htripolii ! ! intent(in) + , t3ple ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real(kind=4), intent(in) :: thil ! Ice-liquid water potential temp. [ K] @@ -3167,7 +3210,6 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) real(kind=4), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=4), intent(in) :: t1stguess ! 1st. guess for temperature [ K] !----- Local variables. -------------------------------------------------------------! - real(kind=4) :: til ! Ice liquid temperature [ K] real(kind=4) :: deriv ! Function derivative real(kind=4) :: fun ! Function for which we seek a root. real(kind=4) :: funa ! Smallest guess function @@ -3190,20 +3232,6 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) thil2temp = cpdryi * thil * exner return !---------------------------------------------------------------------------------! - elseif (.not. newthermo) then - !---------------------------------------------------------------------------------! - ! There is condensation but we are using the old thermodynamics, which can be ! - ! solved analytically. ! - !---------------------------------------------------------------------------------! - til = cpdryi * thil * exner - if (t1stguess > ttripoli) then - thil2temp = 0.5 & - * (til + sqrt(til * (til + cpdryi4 * (alvl3 * rliq + alvi3 * rice)))) - else - thil2temp = til * ( 1. + (alvl3 * rliq + alvi3 * rice) * htripolii) - end if - return - !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -3373,11 +3401,11 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) write (unit=*,fmt='(a,1x,f12.4)' ) 'Exner [J/kg/K] =',exner write (unit=*,fmt='(a,1x,f12.4)' ) 'rliq [ g/kg] =',1000.*rliq write (unit=*,fmt='(a,1x,f12.4)' ) 'rice [ g/kg] =',1000.*rice - write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ °C] =',t1stguess-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ degC] =',t1stguess-t00 write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome (downdraft values).' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ °C] =',tempa-t00 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ °C] =',tempz-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ degC] =',tempa-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ degC] =',tempz-t00 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',fun write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -3407,34 +3435,26 @@ end function thil2temp ! saturation mixing ratio [kg/kg], keeping pressure constant. This is based on the ! ! ice-liquid potential temperature equation. ! !---------------------------------------------------------------------------------------! - real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) + real(kind=4) function dtempdrs(temp,rliq,rice,rconmin) use rconstants , only : alvl3 & ! intent(in) , alvi3 & ! intent(in) , dcpvl & ! intent(in) , dcpvi & ! intent(in) , cpdry & ! intent(in) - , cpdryi & ! intent(in) - , ttripoli & ! intent(in) - , htripolii ! ! intent(in) + , cpdryi ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real(kind=4), intent(in) :: exner ! Exner function [ J/kg/K] real(kind=4), intent(in) :: thil ! Ice-liquid potential temperature (*) [ K] real(kind=4), intent(in) :: temp ! Temperature [ K] real(kind=4), intent(in) :: rliq ! Liquid mixing ratio [ kg/kg] real(kind=4), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=4), intent(in) :: rconmin ! Min. non-zero condensate mix. ratio [ kg/kg] - !------------------------------------------------------------------------------------! - ! (*) Thil is not used in this formulation but it may be used should you opt by ! - ! other ways to compute theta_il, so don't remove this argument. ! - !------------------------------------------------------------------------------------! !----- Local variables --------------------------------------------------------------! real(kind=4) :: qq ! Enthalpy -- latent heat [ J/kg] real(kind=4) :: qpt ! d(qq)/dT * T [ J/kg] real(kind=4) :: hh ! Enthalpy -- sensible heat [ J/kg] real(kind=4) :: rcon ! Condensate mixing ratio [ kg/kg] - real(kind=4) :: til ! Ice-liquid temperature [ K] !------------------------------------------------------------------------------------! @@ -3453,15 +3473,10 @@ real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! ! Find the enthalpy associated with latent heat and its derivative ! - ! correction. This is dependent on the thermodynamics used. ! + ! correction. ! !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - qpt = dcpvl * rliq + dcpvi * rice - else - qq = alvl3 * rliq + alvi3 * rice - qpt = 0.0 - end if + qq = alvl(temp) * rliq + alvi(temp) * rice + qpt = dcpvl * rliq + dcpvi * rice !---------------------------------------------------------------------------------! @@ -3470,21 +3485,8 @@ real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Decide how to compute, based on the thermodynamics method. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - dtempdrs = - temp * qq / ( rcon * (hh + qq - qpt) ) - else - til = cpdryi * thil * exner - !----- Decide how to compute, based on temperature. ---------------------------! - if (temp > ttripoli) then - dtempdrs = - til * qq / ( rcon * cpdry * (2.*temp-til)) - else - dtempdrs = - til * qq * htripolii / rcon - end if - !------------------------------------------------------------------------------! - end if + !----- Find the temperature derivative wrt saturation mixing ratio. --------------! + dtempdrs = - temp * qq / ( rcon * (hh + qq - qpt) ) !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -3770,7 +3772,7 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) !----- Local variables for iterative method -----------------------------------------! real(kind=4) :: pvap ! Sat. vapour pressure real(kind=4) :: theta ! Potential temperature - real(kind=4) :: deriv ! Function derivative + real(kind=4) :: deriv ! Function derivative real(kind=4) :: funnow ! Function for which we seek a root. real(kind=4) :: funa ! Smallest guess function real(kind=4) :: funz ! Largest guess function @@ -3785,8 +3787,6 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) logical :: converged ! Convergence handle logical :: zside ! Side checker for Regula Falsi logical :: frozen ! Will use ice thermodynamics - !----- Local constants. -------------------------------------------------------------! - logical , parameter :: debug = .false. !------------------------------------------------------------------------------------! @@ -3808,18 +3808,16 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! - if (debug) then - write (unit=36,fmt='(a,1x,i5,1x,2(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & - 'NEWTON: it=',itn,'tlclz=',tlclz-t00,'pvap=',0.01*pvap,'fun=',funnow & - ,'deriv=',deriv - end if + !write (unit=36,fmt='(a,1x,i5,1x,2(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & + ! 'NEWTON: it=',itn,'tlclz=',tlclz-t00,'pvap=',0.01*pvap,'fun=',funnow & + ! ,'deriv=',deriv !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! @@ -3903,11 +3896,9 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! - if (debug) then - write (unit=36,fmt='(a,1x,i5,1x,2(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & - '2NGGSS: tt=',itb,'tlclz=',tlclz-t00,'pvap=',0.01*pvap,'fun=',funz & - ,'delta=',delta - end if + !write (unit=36,fmt='(a,1x,i5,1x,2(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & + ! '2NGGSS: tt=',itb,'tlclz=',tlclz-t00,'pvap=',0.01*pvap,'fun=',funz & + ! ,'delta=',delta !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! @@ -3939,7 +3928,7 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) write (unit=*,fmt='(2(a,1x,es14.7,1x))') 'theta =',theta ,'delta =',delta write (unit=*,fmt='(2(a,1x,es14.7,1x))') 'tlcla =',tlcla ,'funa =',funa write (unit=*,fmt='(2(a,1x,es14.7,1x))') 'tlclz =',tlclz ,'funz =',funz - call abort_run ('Failed finding the second guess for regula falsi' & + call fatal_error('Failed finding the second guess for regula falsi' & ,'thetaeiv2thil','therm_lib.f90') end if end if @@ -3956,10 +3945,8 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! - if (debug) then - write (unit=36,fmt='(a,1x,i5,1x,3(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & - 'ANSWER: itb=',itn,'tlcl=',tlcl-t00,'eslcl=',0.01*pvap & - ,'thil=',thetaeiv2thil,'funa=',funa,'funz=',funz - write (unit=36,fmt='(a)') '---------------------------------------------------' - write (unit=36,fmt='(a)') ' ' - end if + !write (unit=36,fmt='(a,1x,i5,1x,3(a,1x,f11.4,1x),2(a,1x,es11.4,1x))') & + ! 'ANSWER: itb=',itn,'tlcl=',tlcl-t00,'eslcl=',0.01*pvap & + ! ,'thil=',thetaeiv2thil,'funa=',funa,'funz=',funz + !write (unit=36,fmt='(a)') '-------------------------------------------------------' + !write (unit=36,fmt='(a)') ' ' !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! else @@ -4014,9 +3999,9 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) write (unit=*,fmt='(a,1x,i12)') ' ITERATIONS :',itb write (unit=*,fmt='(a,1x,f12.5)') ' PVAP [ hPa]:',pvap write (unit=*,fmt='(a,1x,f12.5)') ' THETA [ K]:',theta - write (unit=*,fmt='(a,1x,f12.5)') ' TLCL [ °C]:',tlcl-t00 - write (unit=*,fmt='(a,1x,f12.5)') ' TLCLA [ °C]:',tlcla-t00 - write (unit=*,fmt='(a,1x,f12.5)') ' TLCLZ [ °C]:',tlclz-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCL [ degC]:',tlcl-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCLA [ degC]:',tlcla-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCLZ [ degC]:',tlclz-t00 write (unit=*,fmt='(a,1x,es12.5)') ' FUNA [ K]:',funa write (unit=*,fmt='(a,1x,es12.5)') ' FUNZ [ K]:',funz write (unit=*,fmt='(a,1x,es12.5)') ' DERIV [ ---]:',deriv @@ -4025,11 +4010,12 @@ real(kind=4) function thetaeiv2thil(theiv,pres,rtot,useice) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(60a1)') ('-',ii=1,60) - call abort_run ('TLCL didn''t converge, qgave up!' & + call fatal_error('TLCL didn''t converge, qgave up!' & ,'thetaeiv2thil','therm_lib.f90') end if return + end function thetaeiv2thil !=======================================================================================! !=======================================================================================! @@ -4099,7 +4085,7 @@ subroutine thetaeivs2temp(theivs,pres,theta,temp,rsat,useice) !------------------------------------------------------------------------------------! - ! The 1st. guess, no idea, guess 0°C. ! + ! The 1st. guess, no idea, guess 0 degC. ! !------------------------------------------------------------------------------------! tempz = t00 theta = tempz * exnernormi @@ -4555,13 +4541,13 @@ subroutine lcl_il(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a,1x,f12.4)' ) 'theta_il [ K] =',thil write (unit=*,fmt='(a,1x,f12.4)' ) 'Pressure [ hPa] =',0.01*pres - write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ °C] =',temp-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ degC] =',temp-t00 write (unit=*,fmt='(a,1x,f12.4)' ) 'rtot [ g/kg] =',1000.*rtot write (unit=*,fmt='(a,1x,f12.4)' ) 'rvap [ g/kg] =',1000.*rvap write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome.' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ °C] =',tlcla-t00 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ °C] =',tlclz-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ degC] =',tlcla-t00 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ degC] =',tlclz-t00 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',funnow write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -4569,7 +4555,7 @@ subroutine lcl_il(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a,1x,es12.4)') 'toler [ ----] =',toler write (unit=*,fmt='(a,1x,es12.4)') 'error [ ----] =' & ,abs(tlclz-tlcla)/tlclz - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ °C] =',tlcl + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ degC] =',tlcl call abort_run ('TLCL didn''t converge, gave up!','lcl_il','therm_lib.f90') end if return @@ -4618,7 +4604,7 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) !----- Local variables --------------------------------------------------------------! real(kind=4) :: tempa ! Lower bound for regula falsi iteration real(kind=4) :: tempz ! Upper bound for regula falsi iteration - real(kind=4) :: t1stguess ! Book keeping temperature 1st guess + real(kind=4) :: t1stguess ! Book keeping temperature 1st guess real(kind=4) :: fun1st ! Book keeping 1st guess function real(kind=4) :: funa ! Function evaluation at tempa real(kind=4) :: funz ! Function evaluation at tempz @@ -4632,8 +4618,6 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) logical :: zside ! Aux. Flag, for two purposes: ! 1. Found a 2nd guess for regula falsi. ! 2. I retained the "zside" (T/F) - !----- Local constants. -------------------------------------------------------------! - logical , parameter :: debug = .false. !------------------------------------------------------------------------------------! t1stguess = temp @@ -4683,12 +4667,7 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! - if (debug) then - write (unit=46,fmt='(a,1x,i5,1x,6(a,1x,f11.4,1x),a,1x,es11.4,1x)') & - 'NEWTON: it=',itn,'temp=',tempz-t00,'rsat=',1000.*rsat,'rliq=',1000.*rliq & - ,'rice=',1000.*rice,'rvap=',1000.*rvap,'fun=',funnow,'deriv=',deriv - end if + !write (unit=46,fmt='(a,1x,i5,1x,6(a,1x,f11.4,1x),a,1x,es11.4,1x)') & + ! 'NEWTON: it=',itn,'temp=',tempz-t00,'rsat=',1000.*rsat,'rliq=',1000.*rliq & + ! ,'rice=',1000.*rice,'rvap=',1000.*rvap,'fun=',funnow,'deriv=',deriv !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! - + converged = abs(tempa-tempz) < toler*tempz !---------------------------------------------------------------------------------! ! Convergence. The temperature will be the mid-point between tempa and tempz. ! @@ -4815,14 +4790,14 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) !---------------------------------------------------------------------------------! else if (abs(funnow-funa) < 100.*toler*tempa) then - delta = 0.5 + delta = 100.*toler*tempa else - delta = max(abs(funa)*abs((tempz-tempa)/(funnow-funa)),0.5) + delta = max(abs(funa)*abs((tempz-tempa)/(funnow-funa)),100.*toler*tempa) end if tempz = tempa + delta funz = funa !----- Just to enter at least once. The 1st time tempz=tempa-2*delta ----------! - zside = .false. + zside = .false. zgssloop: do itb=1,maxfpo tempz = tempa + real((-1)**itb * (itb+3)/2) * delta rsat = max(toodry,rslif(pres,tempz)) @@ -4856,7 +4831,7 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) write (unit=*,fmt='(a,1x,f12.5)') ' DELTA [ K]:',delta write (unit=*,fmt='(60a1)') ('-',ii=1,60) - call abort_run ('Failed finding the second guess for regula falsi' & + call fatal_error('Failed finding the second guess for regula falsi' & ,'thil2tqall','therm_lib.f90') end if end if @@ -4892,12 +4867,10 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> Output: ' write (unit=*,fmt='(a,1x,i12)') ' ITERATIONS :',itb - write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ °C]:',temp-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ degC]:',temp-t00 write (unit=*,fmt='(a,1x,f12.5)') ' RVAP [ g/kg]:',1000.*rvap write (unit=*,fmt='(a,1x,f12.5)') ' RLIQ [ g/kg]:',1000.*rliq write (unit=*,fmt='(a,1x,f12.5)') ' RICE [ g/kg]:',1000.*rice - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ °C]:',tempa-t00 - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ °C]:',tempz-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ degC]:',tempa-t00 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ degC]:',tempz-t00 write (unit=*,fmt='(a,1x,es12.5)') ' FUNA [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' FUNZ [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' DERIV [ ---]:',deriv @@ -4986,20 +4959,18 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) write (unit=*,fmt='(a,1x,es12.5)') ' ERR_Z [ ---]:',abs(temp-tempz)/temp write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(60a1)') ('-',ii=1,60) - call abort_run ('Failed finding equilibrium, I gave up!','thil2tqall' & + call fatal_error('Failed finding equilibrium, I gave up!','thil2tqall' & ,'therm_lib.f90') end if !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> maxfpo & ! intent(in) , maxit4 => maxit & ! intent(in) , maxlev4 => maxlev & ! intent(in) - , newthermo4 => newthermo & ! intent(in) , fthva_rp4 => fthva_rp ! ! intent(in) !---------------------------------------------------------------------------------------! @@ -35,8 +34,6 @@ module therm_lib8 integer, parameter :: maxlev = maxlev4 ! Maximum # of levels for adaptive ! quadrature methods. - logical, parameter :: newthermo = newthermo4 ! Use new thermodynamics [T|F] - !---------------------------------------------------------------------------------------! ! This is the "level" variable, that used to be in micphys. Since it affects more the ! ! thermodynamics choices than the microphysics, it was moved to here. ! @@ -135,10 +132,18 @@ module therm_lib8 !=======================================================================================! !=======================================================================================! ! This function calculates the liquid saturation vapour pressure as a function of ! - ! Kelvin temperature. This expression came from MK05, equation (10). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function eslf8(temp,l1funout,l2funout,ttfunout) - use rconstants , only : t008 ! ! intent(in) + use rconstants , only : t008 & ! intent(in) + , es3ple8 & ! intent(in) + , t3ple8 & ! intent(in) + , dcplorv8 & ! intent(in) + , lvorvot38 & ! intent(in) + , rh2o8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=8), intent(in) :: temp ! Temperature [ K] @@ -147,39 +152,66 @@ real(kind=8) function eslf8(temp,l1funout,l2funout,ttfunout) real(kind=8), intent(out), optional :: ttfunout ! Interpolation function real(kind=8), intent(out), optional :: l2funout ! Function for low temperatures !----- Local variables. -------------------------------------------------------------! - real(kind=8) :: l1fun ! - real(kind=8) :: ttfun ! - real(kind=8) :: l2fun ! - real(kind=8) :: x ! + real(kind=8) :: l1fun ! + real(kind=8) :: ttfun ! + real(kind=8) :: l2fun ! + real(kind=8) :: x ! !------------------------------------------------------------------------------------! !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !----- Original method, using polynomial fit (FWC92) -----------------------------! + x = max(-8.0d1,temp-t008) + eslf8 = cll8(0) + x * (cll8(1) + x * (cll8(2) + x * (cll8(3) + x * (cll8(4) & + + x * (cll8(5) + x * (cll8(6) + x * (cll8(7) + x * cll8(8)) )))))) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! + if (present(l1funout)) l1funout = eslf8 + if (present(l2funout)) l2funout = eslf8 + if (present(ttfunout)) ttfunout = eslf8 + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! l1fun = l01_108(0) + l01_108(1)/temp + l01_108(2)*log(temp) + l01_108(3) * temp l2fun = l02_108(0) + l02_108(1)/temp + l02_108(2)*log(temp) + l02_108(3) * temp ttfun = tanh(ttt_108(1) * (temp - ttt_108(2))) eslf8 = exp(l1fun + ttfun*l2fun) !---------------------------------------------------------------------------------! + !----- These are useful for the derivative calculation. --------------------------! if (present(l1funout)) l1funout = l1fun if (present(l2funout)) l2funout = l2fun if (present(ttfunout)) ttfunout = ttfun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.0d1,temp-t008) - eslf8 = cll8(0) + x * (cll8(1) + x * (cll8(2) + x * (cll8(3) + x * (cll8(4) & - + x * (cll8(5) + x * (cll8(6) + x * (cll8(7) + x * cll8(8)) )))))) !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + l1fun = min(lnexp_max8,max(lnexp_min8, lvorvot38 - alvl8(temp) / (rh2o8 * temp))) + eslf8 = es3ple8 * (t3ple8/temp) ** dcplorv8 * exp(l1fun) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(l1funout)) l1funout = eslf8 if (present(l2funout)) l2funout = eslf8 if (present(ttfunout)) ttfunout = eslf8 - end if + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! return end function eslf8 @@ -194,10 +226,18 @@ end function eslf8 !=======================================================================================! !=======================================================================================! ! This function calculates the ice saturation vapour pressure as a function of ! - ! Kelvin temperature, based on MK05 equation (7). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function esif8(temp,iifunout) - use rconstants , only : t008 ! ! intent(in) + use rconstants , only : t008 & ! intent(in) + , es3ple8 & ! intent(in) + , t3ple8 & ! intent(in) + , dcpiorv8 & ! intent(in) + , lsorvot38 & ! intent(in) + , rh2o8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=8), intent(in) :: temp ! Temperature [ K] @@ -209,27 +249,51 @@ real(kind=8) function esif8(temp,iifunout) !------------------------------------------------------------------------------------! + !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then + select case (esat_method) + case (0) + !----- Original method, using polynomial fit (FWC92) -----------------------------! + x = max(-8.d1,temp-t008) + esif8 = cii8(0) + x * (cii8(1) + x * (cii8(2) + x * (cii8(3) + x * (cii8(4) & + + x * (cii8(5) + x * (cii8(6) + x * (cii8(7) + x * cii8(8)))))))) + !---------------------------------------------------------------------------------! - !----- Updated method, using MK05 ------------------------------------------------! + !----- Dummy values for optional function in case they are requested. ------------! + if (present(iifunout)) iifunout=esif8 + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! iifun = iii_78(0) + iii_78(1)/temp + iii_78(2) * log(temp) + iii_78(3) * temp esif8 = exp(iifun) !---------------------------------------------------------------------------------! + !----- This is useful for the derivative calculation. ------==--------------------! if (present(iifunout)) iifunout=iifun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.d1,temp-t008) - esif8 = cii8(0) + x * (cii8(1) + x * (cii8(2) + x * (cii8(3) + x * (cii8(4) & - + x * (cii8(5) + x * (cii8(6) + x * (cii8(7) + x * cii8(8)))))))) + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + iifun = min(lnexp_max8,max(lnexp_min8, lsorvot38 - alvi8(temp) / (rh2o8 * temp))) + esif8 = es3ple8 * (t3ple8/temp) ** dcpiorv8 * exp(iifun) !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(iifunout)) iifunout=esif8 - end if + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -715,7 +779,9 @@ end function rhovsil8 ! pressure with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function eslfp8(temp) - use rconstants , only : t008 ! ! intent(in) + use rconstants , only : t008 & ! intent(in) + , rh2o8 ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=8), intent(in) :: temp @@ -732,24 +798,39 @@ real(kind=8) function eslfp8(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-8.d1,temp-t008) + eslfp8 = dll8(0) + x * (dll8(1) + x * (dll8(2) + x * (dll8(3) + x * (dll8(4) & + + x * (dll8(5) + x * (dll8(6) + x * (dll8(7) + x * dll8(8)))))))) + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! esl = eslf8(temp,l2funout=l2fun,ttfunout=ttfun) l1prime = -l01_108(1)/(temp*temp) + l01_108(2)/temp + l01_108(3) l2prime = -l02_108(1)/(temp*temp) + l02_108(2)/temp + l02_108(3) ttprime = ttt_108(1)*(1.d0 - ttfun*ttfun) eslfp8 = esl * (l1prime + l2prime*ttfun + l2fun*ttprime) - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.d1,temp-t008) - eslfp8 = dll8(0) + x * (dll8(1) + x * (dll8(2) + x * (dll8(3) + x * (dll8(4) & - + x * (dll8(5) + x * (dll8(6) + x * (dll8(7) + x * dll8(8)))))))) - end if + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esl = eslf8(temp) + eslfp8 = esl * alvl8(temp) / (rh2o8 * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! - return end function eslfp8 !=======================================================================================! @@ -766,7 +847,9 @@ end function eslfp8 ! with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function esifp8(temp) - use rconstants , only : t008 ! ! intent(in) + use rconstants , only : t008 & ! intent(in) + , rh2o8 ! ! intent(in) + use mem_basic , only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=8), intent(in) :: temp @@ -779,19 +862,35 @@ real(kind=8) function esifp8(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! - esi = esif8(temp) - iiprime = -iii_78(1)/(temp*temp) + iii_78(2)/temp + iii_78(3) - esifp8 = esi * iiprime - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! x = max(-8.d1,temp-t008) esifp8 = dii8(0) + x * (dii8(1) + x * (dii8(2) + x * (dii8(3) + x * (dii8(4) & + x * (dii8(5) + x * (dii8(6) + x * (dii8(7) + x * dii8(8)))))))) - end if + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! + esi = esif8(temp) + iiprime = -iii_78(1)/(temp*temp) + iii_78(2)/temp + iii_78(3) + esifp8 = esi * iiprime + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esi = esif8(temp) + esifp8 = esi * alvi8(temp) / (rh2o8 * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -2255,7 +2354,7 @@ real(kind=8) function vpdefil8(pres,temp,humi,is_shv,useice) if (frozen) then psat = esif8(temp) else - psat = esif8(temp) + psat = eslf8(temp) end if !------------------------------------------------------------------------------------! @@ -2728,9 +2827,9 @@ end function tl2uint8 !=======================================================================================! !=======================================================================================! - ! This function computes the internal energy of water [J/m²] or [ J/m³], given the ! - ! temperature [K], the heat capacity of the "dry" part [J/m²/K] or [J/m³/K], water mass ! - ! [ kg/m²] or [ kg/m³], and liquid fraction [---]. ! + ! This function computes the internal energy of water [J/m2] or [ J/m3], given the ! + ! temperature [K], the heat capacity of the "dry" part [J/m2/K] or [J/m3/K], water mass ! + ! [ kg/m2] or [ kg/m3], and liquid fraction [---]. ! !---------------------------------------------------------------------------------------! real(kind=8) function cmtl2uext8(dryhcap,wmass,temp,fliq) use rconstants , only : cice8 & ! intent(in) @@ -2739,8 +2838,8 @@ real(kind=8) function cmtl2uext8(dryhcap,wmass,temp,fliq) implicit none !----- Arguments. -------------------------------------------------------------------! - real(kind=8), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m²/K] or [J/m³/K] - real(kind=8), intent(in) :: wmass ! Mass [ kg/m²] or [ kg/m³] + real(kind=8), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m2/K] or [J/m3/K] + real(kind=8), intent(in) :: wmass ! Mass [ kg/m2] or [ kg/m3] real(kind=8), intent(in) :: temp ! Temperature [ K] real(kind=8), intent(in) :: fliq ! Liquid fraction (0-1) [ ---] !------------------------------------------------------------------------------------! @@ -2945,10 +3044,7 @@ end function alvi8 real(kind=8) function theta_iceliq8(exner,temp,rliq,rice) use rconstants , only : alvl38 & ! intent(in) , alvi38 & ! intent(in) - , cpdry8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , cpdry8 ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -2967,36 +3063,13 @@ real(kind=8) function theta_iceliq8(exner,temp,rliq,rice) !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use the ! - ! latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - !------------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - else - qq = alvl38 * rliq + alvi38 * rice - end if + !----- Find the latent heat enthalpy. -----------------------------------------------! + qq = alvl8(temp) * rliq + alvi8(temp) * rice !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Solve the thermodynamics. For the new thermodynamics we don't approximate ! - ! the exponential to a linear function, nor do we impose temperature above the thre- ! - ! shold from Tripoli and Cotton (1981). ! - !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Decide how to compute, based on temperature. ------------------------------! - theta_iceliq8 = hh * exp(-qq / hh) / exner - !---------------------------------------------------------------------------------! - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli8) then - theta_iceliq8 = hh * hh / (exner * ( hh + qq)) - else - theta_iceliq8 = hh * htripoli8 / (exner * ( htripoli8 + qq)) - end if - !---------------------------------------------------------------------------------! - end if + !----- Find theta_iceliq. -----------------------------------------------------------! + theta_iceliq8 = hh * exp(-qq / hh) / exner !------------------------------------------------------------------------------------! return @@ -3020,9 +3093,6 @@ real(kind=8) function dthetail_dt8(condconst,thil,exner,pres,temp,rliq,ricein) , dcpvi8 & ! intent(in) , dcpvl8 & ! intent(in) , cpdry8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripoli8 & ! intent(in) - , htripolii8 & ! intent(in) , t3ple8 ! ! intent(in) implicit none @@ -3036,8 +3106,8 @@ real(kind=8) function dthetail_dt8(condconst,thil,exner,pres,temp,rliq,ricein) real(kind=8), intent(in), optional :: ricein ! Ice mixing ratio [ kg/kg] !----- Local variables --------------------------------------------------------------! real(kind=8) :: rice ! Ice mixing ratio or 0. [ kg/kg] - real(kind=8) :: ldrst ! L × d(rs)/dT × T [ J/kg] - real(kind=8) :: rdlt ! r × d(L)/dT × T [ J/kg] + real(kind=8) :: ldrst ! L * d(rs)/dT * T [ J/kg] + real(kind=8) :: rdlt ! r * d(L)/dT * T [ J/kg] real(kind=8) :: hh ! Sensible heat enthalpy [ J/kg] real(kind=8) :: qq ! Latent heat enthalpy [ J/kg] logical :: thereisice ! Is ice present [ ---] @@ -3075,42 +3145,24 @@ real(kind=8) function dthetail_dt8(condconst,thil,exner,pres,temp,rliq,ricein) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use ! - ! the latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - ! The term r × d(L)/dT × T is computed only when we use the new thermodynamics. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - rdlt = (dcpvl8 * rliq + dcpvi8 * rice ) * temp - else - qq = alvl38 * rliq + alvi38 * rice - rdlt = 0.d0 - end if + !----- Find the latent heat enthalpy. --------------------------------------------! + qq = alvl8(temp) * rliq + alvi8(temp) * rice + rdlt = (dcpvl8 * rliq + dcpvi8 * rice ) * temp !---------------------------------------------------------------------------------! !---------------------------------------------------------------------------------! - ! This is the term L×[d(rs)/dt]×T. L may be either the vapourisation or ! + ! This is the term Lx[d(rs)/dt]xT. L may be either the vapourisation or ! ! sublimation latent heat, depending on the temperature and whether we are consi- ! - ! dering ice or not. We still need to check whether latent heat is a function of ! - ! temperature or not. Also, if condensation mixing ratio is constant, then this ! - ! term will be always zero. ! + ! dering ice or not. If condensation mixing ratio is constant, then this term ! + ! will be always zero. ! !---------------------------------------------------------------------------------! if (condconst) then ldrst = 0.d0 elseif (thereisice .and. temp < t3ple8) then - if (newthermo) then - ldrst = alvi38 * rsifp8(pres,temp) * temp - else - ldrst = alvi8(temp) * rsifp8(pres,temp) * temp - end if + ldrst = alvi8(temp) * rsifp8(pres,temp) * temp else - if (newthermo) then - ldrst = alvl38 * rslfp8(pres,temp) * temp - else - ldrst = alvl8(temp) * rslfp8(pres,temp) * temp - end if + ldrst = alvl8(temp) * rslfp8(pres,temp) * temp end if !---------------------------------------------------------------------------------! end if @@ -3121,17 +3173,7 @@ real(kind=8) function dthetail_dt8(condconst,thil,exner,pres,temp,rliq,ricein) !------------------------------------------------------------------------------------! ! Find the condensed phase consistent with the thermodynamics used. ! !------------------------------------------------------------------------------------! - if (newthermo) then - dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq - rdlt ) / hh ) / temp - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli8) then - dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq) / (hh+qq) ) / temp - else - dthetail_dt8 = thil * ( 1.d0 + ldrst / (htripoli8 + alvl38 * rliq) ) / temp - end if - !---------------------------------------------------------------------------------! - end if + dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq - rdlt ) / hh ) / temp !------------------------------------------------------------------------------------! return @@ -3159,9 +3201,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) , alvl38 & ! intent(in) , alvi38 & ! intent(in) , t008 & ! intent(in) - , t3ple8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , t3ple8 ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real(kind=8), intent(in) :: thil ! Ice-liquid water potential temp. [ K] @@ -3171,7 +3211,6 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) real(kind=8), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=8), intent(in) :: t1stguess ! 1st. guess for temperature [ K] !----- Local variables. -------------------------------------------------------------! - real(kind=8) :: til ! Ice liquid temperature [ K] real(kind=8) :: deriv ! Function derivative real(kind=8) :: fun ! Function for which we seek a root. real(kind=8) :: funa ! Smallest guess function @@ -3194,21 +3233,6 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) thil2temp8 = cpdryi8 * thil * exner return !---------------------------------------------------------------------------------! - elseif (.not. newthermo) then - !---------------------------------------------------------------------------------! - ! There is condensation but we are using the old thermodynamics, which can be ! - ! solved analytically. ! - !---------------------------------------------------------------------------------! - til = cpdryi8 * thil * exner - if (t1stguess > ttripoli8) then - thil2temp8 = 5.d-1 & - * (til + sqrt( til & - * (til + cpdryi48 * (alvl38 * rliq + alvi38 * rice)))) - else - thil2temp8 = til * ( 1.d0 + (alvl38 * rliq + alvi38 * rice) * htripolii8) - end if - return - !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -3229,7 +3253,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) !------------------------------------------------------------------------------------! tempz = t1stguess fun = theta_iceliq8(exner,tempz,rliq,rice) - deriv = dthetail_dt8(.true.,fun,exner,pres,tempz,rliq,rice) + deriv = dthetail_dt8(.true.,fun,pres,tempz,rliq,rice) fun = fun - thil !------------------------------------------------------------------------------------! @@ -3244,7 +3268,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) thil2temp8 = tempz converged = .true. return - else + else tempa = tempz funa = fun end if @@ -3257,7 +3281,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) funa = fun tempz = tempa - fun/deriv fun = theta_iceliq8(exner,tempz,rliq,rice) - deriv = dthetail_dt8(.true.,fun,exner,pres,tempz,rliq,rice) + deriv = dthetail_dt8(.true.,fun,pres,tempz,rliq,rice) fun = fun - thil !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! @@ -3293,7 +3317,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) else if (abs(fun-funa) < toler8 * tempa) then delta = 1.d2 * toler8 * tempa - else + else delta = max(abs(funa * (tempz-tempa)/(fun-funa)),1.d2 * toler8 * tempa) end if tempz = tempa + delta @@ -3313,7 +3337,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) write (unit=*,fmt='(2(a,1x,es14.7,1x))') 'tempa=',tempa,'funa=',funa write (unit=*,fmt='(2(a,1x,es14.7,1x))') 'tempz=',tempz,'funz=',funz write (unit=*,fmt='(1(a,1x,es14.7,1x))') 'delta=',delta - call abort_run ('Failed finding the second guess for regula falsi' & + call fatal_error('Failed finding the second guess for regula falsi' & ,'thil2temp8','therm_lib8.f90') end if end if @@ -3378,11 +3402,11 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) write (unit=*,fmt='(a,1x,f12.4)' ) 'Exner [J/kg/K] =',exner write (unit=*,fmt='(a,1x,f12.4)' ) 'rliq [ g/kg] =',1.d3*rliq write (unit=*,fmt='(a,1x,f12.4)' ) 'rice [ g/kg] =',1.d3*rice - write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ °C] =',t1stguess-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ degC] =',t1stguess-t008 write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome (downdraft values).' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ °C] =',tempa-t008 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ °C] =',tempz-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ degC] =',tempa-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ degC] =',tempz-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',fun write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -3392,7 +3416,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) , abs(thil2temp8-tempa)/thil2temp8 write (unit=*,fmt='(a,1x,f12.4)' ) 'thil2temp8 [ K] =',thil2temp8 - call abort_run ('Temperature didn''t converge, giving up!!!' & + call fatal_error('Temperature didn''t converge, giving up!!!' & ,'thil2temp8','therm_lib8.f90') end if @@ -3412,34 +3436,25 @@ end function thil2temp8 ! saturation mixing ratio [kg/kg], keeping pressure constant. This is based on the ! ! ice-liquid potential temperature equation. ! !---------------------------------------------------------------------------------------! - real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) + real(kind=8) function dtempdrs8(temp,rliq,rice,rconmin) use rconstants , only : alvl38 & ! intent(in) , alvi38 & ! intent(in) , dcpvl8 & ! intent(in) , dcpvi8 & ! intent(in) , cpdry8 & ! intent(in) - , cpdryi8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , cpdryi8 ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real(kind=8), intent(in) :: exner ! Exner function [ J/kg/K] - real(kind=8), intent(in) :: thil ! Ice-liquid potential temperature (*) [ K] real(kind=8), intent(in) :: temp ! Temperature [ K] real(kind=8), intent(in) :: rliq ! Liquid mixing ratio [ kg/kg] real(kind=8), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=8), intent(in) :: rconmin ! Min. non-zero condensate mix. ratio [ kg/kg] - !------------------------------------------------------------------------------------! - ! (*) Thil is not used in this formulation but it may be used should you opt by ! - ! other ways to compute theta_il, so don't remove this argument. ! - !------------------------------------------------------------------------------------! !----- Local variables --------------------------------------------------------------! real(kind=8) :: qq ! Enthalpy -- latent heat [ J/kg] real(kind=8) :: qpt ! d(qq)/dT * T [ J/kg] real(kind=8) :: hh ! Enthalpy -- sensible heat [ J/kg] real(kind=8) :: rcon ! Condensate mixing ratio [ kg/kg] - real(kind=8) :: til ! Ice-liquid temperature [ K] !------------------------------------------------------------------------------------! @@ -3458,15 +3473,10 @@ real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! ! Find the enthalpy associated with latent heat and its derivative ! - ! correction. This is dependent on the thermodynamics used. ! + ! correction. ! !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - qpt = dcpvl8 * rliq + dcpvi8 * rice - else - qq = alvl38 * rliq + alvi38 * rice - qpt = 0.d0 - end if + qq = alvl8(temp) * rliq + alvi8(temp) * rice + qpt = dcpvl8 * rliq + dcpvi8 * rice !---------------------------------------------------------------------------------! @@ -3475,21 +3485,8 @@ real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Decide how to compute, based on the thermodynamics method. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - dtempdrs8 = - temp * qq / ( rcon * (hh + qq - qpt) ) - else - til = cpdryi8 * thil * exner - !----- Decide how to compute, based on temperature. ---------------------------! - if (temp > ttripoli8) then - dtempdrs8 = - til * qq / ( rcon * cpdry8 * (2.*temp-til)) - else - dtempdrs8 = - til * qq * htripolii8 / rcon - end if - !------------------------------------------------------------------------------! - end if + !----- Find the temperature derivative wrt saturation mixing ratio. --------------! + dtempdrs8 = - temp * qq / ( rcon * (hh + qq - qpt) ) !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -4002,9 +3999,9 @@ real(kind=8) function thetaeiv2thil8(theiv,pres,rtot,useice) write (unit=*,fmt='(a,1x,i12)') ' ITERATIONS :',itb write (unit=*,fmt='(a,1x,f12.5)') ' PVAP [ hPa]:',pvap write (unit=*,fmt='(a,1x,f12.5)') ' THETA [ K]:',theta - write (unit=*,fmt='(a,1x,f12.5)') ' TLCL [ °C]:',tlcl-t008 - write (unit=*,fmt='(a,1x,f12.5)') ' TLCLA [ °C]:',tlcla-t008 - write (unit=*,fmt='(a,1x,f12.5)') ' TLCLZ [ °C]:',tlclz-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCL [ degC]:',tlcl-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCLA [ degC]:',tlcla-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TLCLZ [ degC]:',tlclz-t008 write (unit=*,fmt='(a,1x,es12.5)') ' FUNA [ K]:',funa write (unit=*,fmt='(a,1x,es12.5)') ' FUNZ [ K]:',funz write (unit=*,fmt='(a,1x,es12.5)') ' DERIV [ ---]:',deriv @@ -4087,7 +4084,7 @@ subroutine thetaeivs2temp8(theivs,pres,theta,temp,rsat,useice) !------------------------------------------------------------------------------------! - ! The 1st. guess, no idea, guess 0°C. ! + ! The 1st. guess, no idea, guess 0 degC. ! !------------------------------------------------------------------------------------! tempz = t008 theta = tempz * exnernormi @@ -4543,13 +4540,13 @@ subroutine lcl_il8(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a,1x,f12.4)' ) 'theta_il [ K] =',thil write (unit=*,fmt='(a,1x,f12.4)' ) 'Pressure [ hPa] =',1.d-2*pres - write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ °C] =',temp-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ degC] =',temp-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'rtot [ g/kg] =',1.d3*rtot write (unit=*,fmt='(a,1x,f12.4)' ) 'rvap [ g/kg] =',1.d3*rvap write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome.' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ °C] =',tlcla-t008 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ °C] =',tlclz-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ degC] =',tlcla-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ degC] =',tlclz-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',funnow write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -4557,7 +4554,7 @@ subroutine lcl_il8(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a,1x,es12.4)') 'toler8 [ ----] =',toler8 write (unit=*,fmt='(a,1x,es12.4)') 'error [ ----] =' & ,abs(tlclz-tlcla)/tlclz - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ °C] =',tlcl + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ degC] =',tlcl call abort_run ('TLCL didn''t converge, gave up!','lcl_il8','therm_lib8.f90') end if return @@ -4590,8 +4587,7 @@ subroutine thil2tqall8(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) , cpdryi8 & ! intent(in) , t008 & ! intent(in) , toodry8 & ! intent(in) - , t3ple8 & ! intent(in) - , ttripoli8 ! ! intent(in) + , t3ple8 ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! @@ -4949,12 +4945,12 @@ subroutine thil2tqall8(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' -> Output: ' write (unit=*,fmt='(a,1x,i12)') ' ITERATIONS :',itb - write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ °C]:',temp-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ degC]:',temp-t008 write (unit=*,fmt='(a,1x,f12.5)') ' RVAP [ g/kg]:',1.d3*rvap write (unit=*,fmt='(a,1x,f12.5)') ' RLIQ [ g/kg]:',1.d3*rliq write (unit=*,fmt='(a,1x,f12.5)') ' RICE [ g/kg]:',1.d3*rice - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ °C]:',tempa-t008 - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ °C]:',tempz-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ degC]:',tempa-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ degC]:',tempz-t008 write (unit=*,fmt='(a,1x,es12.5)') ' FUNA [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' FUNZ [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' DERIV [ ---]:',deriv @@ -5283,8 +5279,8 @@ end subroutine uint2tl8 !=======================================================================================! !=======================================================================================! ! This subroutine computes the temperature (Kelvin) and liquid fraction from ! - ! extensive internal energy (J/m² or J/m³), water mass (kg/m² or kg/m³), and heat ! - ! capacity (J/m²/K or J/m³/K). ! + ! extensive internal energy (J/m2 or J/m3), water mass (kg/m2 or kg/m3), and heat ! + ! capacity (J/m2/K or J/m3/K). ! !---------------------------------------------------------------------------------------! subroutine uextcm2tl8(uext,wmass,dryhcap,temp,fliq) use rconstants , only : cliqi8 & ! intent(in) @@ -5297,19 +5293,19 @@ subroutine uextcm2tl8(uext,wmass,dryhcap,temp,fliq) , tsupercool_liq8 ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real(kind=8), intent(in) :: uext ! Extensive internal energy [ J/m²] or [ J/m³] - real(kind=8), intent(in) :: wmass ! Water mass [ kg/m²] or [ kg/m³] - real(kind=8), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m²/K] or [J/m³/K] + real(kind=8), intent(in) :: uext ! Extensive internal energy [ J/m2] or [ J/m3] + real(kind=8), intent(in) :: wmass ! Water mass [ kg/m2] or [ kg/m3] + real(kind=8), intent(in) :: dryhcap ! Heat cap. of "dry" part [J/m2/K] or [J/m3/K] real(kind=8), intent(out) :: temp ! Temperature [ K] real(kind=8), intent(out) :: fliq ! Liquid fraction (0-1) [ ---] !----- Local variable ---------------------------------------------------------------! - real(kind=8) :: uefroz ! qw of ice at triple pt. [ J/m²] or [ J/m³] - real(kind=8) :: uemelt ! qw of liq. at triple pt. [ J/m²] or [ J/m³] + real(kind=8) :: uefroz ! qw of ice at triple pt. [ J/m2] or [ J/m3] + real(kind=8) :: uemelt ! qw of liq. at triple pt. [ J/m2] or [ J/m3] !------------------------------------------------------------------------------------! - !----- Convert melting heat to J/m² or J/m³ -----------------------------------------! + !----- Convert melting heat to J/m2 or J/m3 -----------------------------------------! uefroz = (dryhcap + wmass * cice8) * t3ple8 uemelt = uefroz + wmass * alli8 !------------------------------------------------------------------------------------! diff --git a/BRAMS/src/memory/mem_basic.f90 b/BRAMS/src/memory/mem_basic.f90 index 9a315af74..8fa09a37a 100644 --- a/BRAMS/src/memory/mem_basic.f90 +++ b/BRAMS/src/memory/mem_basic.f90 @@ -48,6 +48,29 @@ module mem_basic !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + integer :: esat_method + !---------------------------------------------------------------------------------------! + + !----- This variable is just an alias to set CO2 or not. -------------------------------! logical :: co2_on !---------------------------------------------------------------------------------------! diff --git a/BRAMS/src/micro/mic_misc.f90 b/BRAMS/src/micro/mic_misc.f90 index 230086fa6..ec0082275 100644 --- a/BRAMS/src/micro/mic_misc.f90 +++ b/BRAMS/src/micro/mic_misc.f90 @@ -416,7 +416,7 @@ subroutine fill_thermovars(m1,i,j,flpw,thp,btheta,pp,rtp,rv,wp,dn0,pi0,micro) tairstr(k) = tair(k) pottemp(k) = extemp2theta(exner(k),tair(k)) !----- The A1 term in Walko et al. (2000) according to the new thermodynamics. ------! - sa1(k) = (-1.) * dtempdrs(exner(k),thil(k),tairstr(k),rliq(k),rice(k),1.e-12) + sa1(k) = (-1.) * dtempdrs(tairstr(k),rliq(k),rice(k),1.e-12) !------------------------------------------------------------------------------------! ! Ice-vapour equivalent potential temperature. This variable is conserved even ! @@ -1032,7 +1032,7 @@ subroutine x02(lcat) rx(k,lcat) = rx(k,lcat) - rmelt - ricetor6 rx(k,6) = rx(k,6) + rmelt + ricetor6 qr(k,6) = qr(k,6) + rmelt * uiliqt3 + ricetor6 * uiicet3 - qx(k,lcat) = uiicet3 !---- All water is gone, leave ice at 0°C only --------! + qx(k,lcat) = uiicet3 !---- All water is gone, leave ice at 0C only ---------! qr(k,lcat) = qx(k,lcat) * rx(k,lcat) if (rx(k,6) > rxmin(6)) qx(k,6) = qr(k,6) / rx(k,6) diff --git a/BRAMS/src/mpi/mpass_init.F90 b/BRAMS/src/mpi/mpass_init.F90 index 76e524e79..f00e0d787 100644 --- a/BRAMS/src/mpi/mpass_init.F90 +++ b/BRAMS/src/mpi/mpass_init.F90 @@ -422,6 +422,7 @@ subroutine masterput_nl(master_num) call MPI_Bcast(DTLEAF,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(NDTVEG,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(ISTAR,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(ESAT_METHOD,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(IGRNDVAP,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(ICO2,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) @@ -1302,6 +1303,7 @@ subroutine nodeget_nl call MPI_Bcast(DTLEAF,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(NDTVEG,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(ISTAR,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(ESAT_METHOD,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(IGRNDVAP,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(ICO2,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) diff --git a/ED/build/make/include.mk.macos b/ED/build/make/include.mk.macos index 9d0f940d1..452e33ce4 100644 --- a/ED/build/make/include.mk.macos +++ b/ED/build/make/include.mk.macos @@ -100,8 +100,8 @@ CMACH=MACOS # MOD_EXT -- Suffix for external modules. This should be typically set to mod. # #------------------------------------------------------------------------------------------# FC_TYPE=GNU -F_COMP=$(HOMEBREW_PREFIX)/bin/gfortran-15 -C_COMP=$(HOMEBREW_PREFIX)/bin/gcc-15 +F_COMP=$(HOMEBREW_PREFIX)/bin/gfortran-16 +C_COMP=$(HOMEBREW_PREFIX)/bin/gcc-16 LOADER=$(F_COMP) LIBS= MOD_EXT=mod @@ -119,16 +119,20 @@ MOD_EXT=mod ifeq ($(KIND_COMP),$(filter $(KIND_COMP), A B C D)) F_OPTS= -O0 -ffree-line-length-none -g -fimplicit-none -Wall -finit-real=snan \ -finit-integer=-2147483648 -ffpe-trap=invalid,zero,overflow,underflow \ - -fcheck=all -frecursive -fsignaling-nans -Werror -fopenmp -fbacktrace -static + -fcheck=all -fsignaling-nans -Werror -fopenmp -fbacktrace -static \ + -Wl,-stack_size,0x20000000 C_OPTS= -O0 -DLITTLE -g -fbacktrace -static LOADER_OPTS= -O0 -ffree-line-length-none -g -fimplicit-none -Wall -finit-real=snan \ -finit-integer=-2147483648 -ffpe-trap=invalid,zero,overflow,underflow \ - -fcheck=all -frecursive -fsignaling-nans -Werror -fopenmp -fbacktrace + -fcheck=all -fsignaling-nans -Werror -fopenmp -fbacktrace \ + -Wl,-stack_size,0x20000000 #---------------------------------------------------------------------------------------# else ifeq ($(KIND_COMP),E) - F_OPTS= -O3 -ffree-line-length-none -frecursive -fopenmp -fbacktrace -static + F_OPTS= -O3 -ffree-line-length-none -fopenmp -fbacktrace -static \ + -Wl,-stack_size,0x20000000 C_OPTS= -O3 -DLITTLE -g -fbacktrace -static - LOADER_OPTS= -O3 -ffree-line-length-none -frecursive -fopenmp -fbacktrace + LOADER_OPTS= -O3 -ffree-line-length-none -fopenmp -fbacktrace \ + -Wl,-stack_size,0x20000000 #---------------------------------------------------------------------------------------# else $(error Option KIND_COMP provided ($(KIND_COMP)) is invalid)." diff --git a/ED/run/ED2IN b/ED/run/ED2IN index b7df194cd..d4c967a6e 100644 --- a/ED/run/ED2IN +++ b/ED/run/ED2IN @@ -795,6 +795,28 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + NL%ESAT_METHOD = 1 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. (Deprecated) Euler step. The fastest, but it has only a ! diff --git a/ED/src/io/ed_load_namelist.f90 b/ED/src/io/ed_load_namelist.f90 index cd181df0b..e99f1b219 100644 --- a/ED/src/io/ed_load_namelist.f90 +++ b/ED/src/io/ed_load_namelist.f90 @@ -170,6 +170,7 @@ subroutine copy_nl(copy_type) , radfrq & ! intent(out) , ivegt_dynamics & ! intent(out) , ibigleaf & ! intent(out) + , esat_method & ! intent(out) , integration_scheme & ! intent(out) , nsub_euler & ! intent(out) , ffilout & ! intent(out) @@ -372,6 +373,7 @@ subroutine copy_nl(copy_type) ivegt_dynamics = nl%ivegt_dynamics ibigleaf = nl%ibigleaf + esat_method = nl%esat_method integration_scheme = nl%integration_scheme nsub_euler = nl%nsub_euler rk4_tolerance = nl%rk4_tolerance diff --git a/ED/src/io/ed_opspec.F90 b/ED/src/io/ed_opspec.F90 index 2cff44de0..aa105f36c 100644 --- a/ED/src/io/ed_opspec.F90 +++ b/ED/src/io/ed_opspec.F90 @@ -1175,6 +1175,7 @@ subroutine ed_opspec_misc , ied_init_mode & ! intent(in) , ivegt_dynamics & ! intent(in) , ibigleaf & ! intent(in) + , esat_method & ! intent(in) , integration_scheme & ! intent(in) , nsub_euler & ! intent(in) , iallom & ! intent(in) @@ -1678,6 +1679,14 @@ subroutine ed_opspec_misc ifaterr = ifaterr +1 end if + if (esat_method < 0 .or. esat_method > 2) then + write (reason,fmt='(a,1x,i4,a)') & + 'Invalid ESAT_METHOD, it must be between 0 and 2. Yours is set to' & + ,esat_method,'...' + call opspec_fatal(reason,'opspec_misc') + ifaterr = ifaterr +1 + end if + !---------------------------------------------------------------------------------------! ! Integration scheme can be only 0 (Euler) or 1 (4th order Runge-Kutta). The ! ! branch thermodynamics is currently working only with Runge-Kutta, so we won't allow ! diff --git a/ED/src/io/ed_xml_config.f90 b/ED/src/io/ed_xml_config.f90 index 33a1d72ac..bee088fda 100644 --- a/ED/src/io/ed_xml_config.f90 +++ b/ED/src/io/ed_xml_config.f90 @@ -268,6 +268,8 @@ recursive subroutine read_ed_xml_config(filename) if(texist) ivegt_dynamics = ival call getConfigINT ('ibigleaf','misc',i,ival,texist) if(texist) ibigleaf = ival + call getConfigINT ('esat_method','misc',i,ival,texist) + if(texist) esat_method = ival call getConfigINT ('integration_scheme','misc',i,ival,texist) if(texist) integration_scheme = ival @@ -1939,6 +1941,7 @@ subroutine write_ed_xml_config call putConfigSTRING("history_out_filepath",sfilout) call putConfigINT("ivegt_dynamics",ivegt_dynamics) call putConfigINT("ibigleaf",ibigleaf) + call putConfigINT("esat_method",esat_method) call putConfigINT("integration_scheme",integration_scheme) call libxml2f90_ll_closetag("misc") diff --git a/ED/src/memory/consts_coms.F90 b/ED/src/memory/consts_coms.F90 index 7fe069102..58be21598 100644 --- a/ED/src/memory/consts_coms.F90 +++ b/ED/src/memory/consts_coms.F90 @@ -143,6 +143,10 @@ Module consts_coms , b_uiliqt3 => uiliqt3 & ! intent(in) , b_dcpvl => dcpvl & ! intent(in) , b_dcpvi => dcpvi & ! intent(in) + , b_dcplorv => dcplorv & ! intent(in) + , b_dcpiorv => dcpiorv & ! intent(in) + , b_lvorvot3 => lvorvot3 & ! intent(in) + , b_lsorvot3 => lsorvot3 & ! intent(in) , b_del_alvl3 => del_alvl3 & ! intent(in) , b_del_alvi3 => del_alvi3 & ! intent(in) , b_tsupercool_liq => tsupercool_liq & ! intent(in) @@ -318,6 +322,10 @@ Module consts_coms real , parameter :: uiliqt3 = b_uiliqt3 real , parameter :: dcpvl = b_dcpvl real , parameter :: dcpvi = b_dcpvi + real , parameter :: dcplorv = b_dcplorv + real , parameter :: dcpiorv = b_dcpiorv + real , parameter :: lvorvot3 = b_lvorvot3 + real , parameter :: lsorvot3 = b_lsorvot3 real , parameter :: del_alvl3 = b_del_alvl3 real , parameter :: del_alvi3 = b_del_alvi3 real , parameter :: tsupercool_liq = b_tsupercool_liq @@ -396,7 +404,7 @@ Module consts_coms !---------------------------------------------------------------------------------------! real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m2/K4] real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m2 kg/s2/K] - real, parameter :: t00 = 273.15 ! 0 degC [ degC] + real, parameter :: t00 = 273.15 ! 0 degC [ K] real, parameter :: rmol = 8.314510 ! Molar gas constant [ J/mol/K] real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m3] real, parameter :: volmoll = volmol*1e3 ! Molar volume at STP [ L] @@ -567,6 +575,20 @@ Module consts_coms !---------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! + ! These derived variables are used for computing vapour pressure at saturation when ! + ! using the method proposed by Ambaum (2020). ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure over ! + ! water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899. ! + !---------------------------------------------------------------------------------------! + real, parameter :: dcplorv = (cliq - cph2o) / rh2o ! (cl - cp) / Rv [ ------] + real, parameter :: dcpiorv = (cice - cph2o) / rh2o ! (ci - cp) / Rv [ ------] + real, parameter :: lvorvot3 = alvl3 / (rh2o * t3ple) ! Lv0 / ( Rv * T3) [ ------] + real, parameter :: lsorvot3 = alvi3 / (rh2o * t3ple) ! Ls0 / ( Rv * T3) [ ------] + !---------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! ! The following variables are useful when defining the derivatives of theta_il. ! ! They correspond to L?(T) - L?' T. ! @@ -782,6 +804,10 @@ Module consts_coms real(kind=8), parameter :: uiliqt38 = dble(uiliqt3 ) real(kind=8), parameter :: dcpvl8 = dble(dcpvl ) real(kind=8), parameter :: dcpvi8 = dble(dcpvi ) + real(kind=8), parameter :: dcplorv8 = dble(dcplorv ) + real(kind=8), parameter :: dcpiorv8 = dble(dcpiorv ) + real(kind=8), parameter :: lvorvot38 = dble(lvorvot3 ) + real(kind=8), parameter :: lsorvot38 = dble(lsorvot3 ) real(kind=8), parameter :: del_alvl38 = dble(del_alvl3 ) real(kind=8), parameter :: del_alvi38 = dble(del_alvi3 ) real(kind=8), parameter :: tsupercool_liq8 = dble(tsupercool_liq) diff --git a/ED/src/memory/ed_misc_coms.f90 b/ED/src/memory/ed_misc_coms.f90 index f6e496671..cfdae6f6d 100644 --- a/ED/src/memory/ed_misc_coms.f90 +++ b/ED/src/memory/ed_misc_coms.f90 @@ -165,6 +165,29 @@ module ed_misc_coms + + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + integer :: esat_method + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. Euler step. The fastest, but it doesn't estimate ! diff --git a/ED/src/memory/ename_coms.f90 b/ED/src/memory/ename_coms.f90 index e4f7fc55d..9ad61c520 100644 --- a/ED/src/memory/ename_coms.f90 +++ b/ED/src/memory/ename_coms.f90 @@ -146,6 +146,7 @@ module ename_coms !----- Options for model dynamics. --------------------------------------------------! integer :: ivegt_dynamics integer :: ibigleaf + integer :: esat_method integer :: integration_scheme integer :: nsub_euler integer :: plant_hydro_scheme @@ -393,6 +394,7 @@ subroutine init_ename_vars(enl) enl%ivegt_dynamics = undef_integer enl%ibigleaf = undef_integer + enl%esat_method = undef_integer enl%integration_scheme = undef_integer enl%nsub_euler = undef_integer enl%plant_hydro_scheme = undef_integer diff --git a/ED/src/mpi/ed_mpass_init.F90 b/ED/src/mpi/ed_mpass_init.F90 index 24e37d3f1..a25bfe798 100644 --- a/ED/src/mpi/ed_mpass_init.F90 +++ b/ED/src/mpi/ed_mpass_init.F90 @@ -139,6 +139,8 @@ subroutine ed_masterput_nl(par_run) , ied_init_mode & ! intent(in) , thsums_database & ! intent(in) , ivegt_dynamics & ! intent(in) + , ibigleaf & ! intent(in) + , esat_method & ! intent(in) , integration_scheme & ! intent(in) , nsub_euler & ! intent(in) , end_time & ! intent(in) @@ -458,6 +460,8 @@ subroutine ed_masterput_nl(par_run) call MPI_Bcast(centlon,maxgrds,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(ivegt_dynamics,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(ibigleaf,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(esat_method,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(integration_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(nsub_euler,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(rk4_tolerance,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) @@ -1304,6 +1308,8 @@ subroutine ed_nodeget_nl , ied_init_mode & ! intent(out) , thsums_database & ! intent(out) , ivegt_dynamics & ! intent(out) + , ibigleaf & ! intent(out) + , esat_method & ! intent(out) , integration_scheme & ! intent(out) , nsub_euler & ! intent(out) , end_time & ! intent(out) @@ -1627,6 +1633,8 @@ subroutine ed_nodeget_nl call MPI_Bcast(centlon,maxgrds,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(ivegt_dynamics,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(ibigleaf,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(esat_method,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(integration_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(nsub_euler,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(rk4_tolerance,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) diff --git a/ED/src/utils/therm_lib.f90 b/ED/src/utils/therm_lib.f90 index 06197d672..9b293f3ee 100644 --- a/ED/src/utils/therm_lib.f90 +++ b/ED/src/utils/therm_lib.f90 @@ -27,7 +27,6 @@ module therm_lib ! crashing, for other methods. integer , parameter :: maxlev = 16 ! Maximum # of levels for adap- ! tive quadrature methods. - logical , parameter :: newthermo = .true. ! Use new thermodynamics [T|F] !---------------------------------------------------------------------------------------! @@ -133,10 +132,18 @@ module therm_lib !=======================================================================================! !=======================================================================================! ! This function calculates the liquid saturation vapour pressure as a function of ! - ! Kelvin temperature. This expression came from MK05, equation (10). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function eslf(temp,l1funout,l2funout,ttfunout) - use consts_coms, only : t00 ! ! intent(in) + use consts_coms , only : t00 & ! intent(in) + , es3ple & ! intent(in) + , t3ple & ! intent(in) + , dcplorv & ! intent(in) + , lvorvot3 & ! intent(in) + , rh2o & ! intent(in) + , lnexp_min & ! intent(in) + , lnexp_max ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=4), intent(in) :: temp ! Temperature [ K] @@ -154,30 +161,59 @@ real(kind=4) function eslf(temp,l1funout,l2funout,ttfunout) !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-80.,temp-t00) + eslf = cll(0) + x * (cll(1) + x * (cll(2) + x * (cll(3) + x * (cll(4) & + + x * (cll(5) + x * (cll(6) + x * (cll(7) + x * cll(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! + if (present(l1funout)) l1funout = eslf + if (present(l2funout)) l2funout = eslf + if (present(ttfunout)) ttfunout = eslf + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! l1fun = l01_10(0) + l01_10(1)/temp + l01_10(2)*log(temp) + l01_10(3) * temp l2fun = l02_10(0) + l02_10(1)/temp + l02_10(2)*log(temp) + l02_10(3) * temp ttfun = tanh(ttt_10(1) * (temp - ttt_10(2))) eslf = exp(l1fun + ttfun*l2fun) !---------------------------------------------------------------------------------! + !----- These are useful for the derivative calculation. --------------------------! if (present(l1funout)) l1funout = l1fun if (present(l2funout)) l2funout = l2fun if (present(ttfunout)) ttfunout = ttfun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-80.,temp-t00) - eslf = cll(0) + x * (cll(1) + x * (cll(2) + x * (cll(3) + x * (cll(4) & - + x * (cll(5) + x * (cll(6) + x * (cll(7) + x * cll(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + l1fun = min(lnexp_max,max(lnexp_min, lvorvot3 - alvl(temp) / (rh2o * temp))) + eslf = es3ple * (t3ple/temp) ** dcplorv * exp(l1fun) !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(l1funout)) l1funout = eslf if (present(l2funout)) l2funout = eslf if (present(ttfunout)) ttfunout = eslf - end if + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! return end function eslf @@ -192,10 +228,18 @@ end function eslf !=======================================================================================! !=======================================================================================! ! This function calculates the ice saturation vapour pressure as a function of ! - ! Kelvin temperature, based on MK05 equation (7). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function esif(temp,iifunout) - use consts_coms, only : t00 ! ! intent(in) + use consts_coms , only : t00 & ! intent(in) + , es3ple & ! intent(in) + , t3ple & ! intent(in) + , dcpiorv & ! intent(in) + , lsorvot3 & ! intent(in) + , rh2o & ! intent(in) + , lnexp_min & ! intent(in) + , lnexp_max ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=4), intent(in) :: temp ! Temperature [ K] @@ -207,27 +251,53 @@ real(kind=4) function esif(temp,iifunout) !------------------------------------------------------------------------------------! + !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x=max(-80.,temp-t00) + esif = cii(0) + x * (cii(1) + x * (cii(2) + x * (cii(3) + x * (cii(4) & + + x * (cii(5) + x * (cii(6) + x * (cii(7) + x * cii(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! - !----- Updated method, using MK05 ------------------------------------------------! + !----- Dummy values for optional function in case they are requested. ------------! + if (present(iifunout)) iifunout=esif + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! iifun = iii_7(0) + iii_7(1)/temp + iii_7(2) * log(temp) + iii_7(3) * temp esif = exp(iifun) !---------------------------------------------------------------------------------! + !----- This is useful for the derivative calculation. ------==--------------------! if (present(iifunout)) iifunout=iifun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x=max(-80.,temp-t00) - esif = cii(0) + x * (cii(1) + x * (cii(2) + x * (cii(3) + x * (cii(4) & - + x * (cii(5) + x * (cii(6) + x * (cii(7) + x * cii(8)) ) ) ) ) ) ) !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + iifun = min(lnexp_max,max(lnexp_min, lsorvot3 - alvi(temp) / (rh2o * temp))) + esif = es3ple * (t3ple/temp) ** dcpiorv * exp(iifun) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(iifunout)) iifunout=esif - end if + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -713,7 +783,9 @@ end function rhovsil ! pressure with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function eslfp(temp) - use consts_coms, only : t00 ! ! intent(in) + use consts_coms , only : t00 & ! intent(in) + , rh2o ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=4), intent(in) :: temp @@ -730,21 +802,37 @@ real(kind=4) function eslfp(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-80.,temp-t00) + eslfp = dll(0) + x * (dll(1) + x * (dll(2) + x * (dll(3) + x * (dll(4) & + + x * (dll(5) + x * (dll(6) + x * (dll(7) + x * dll(8)) ) ) ) ) ) ) + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! esl = eslf(temp,l2funout=l2fun,ttfunout=ttfun) l1prime = -l01_10(1)/(temp*temp) + l01_10(2)/temp + l01_10(3) l2prime = -l02_10(1)/(temp*temp) + l02_10(2)/temp + l02_10(3) ttprime = ttt_10(1)*(1.-ttfun*ttfun) eslfp = esl * (l1prime + l2prime*ttfun + l2fun*ttprime) - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-80.,temp-t00) - eslfp = dll(0) + x * (dll(1) + x * (dll(2) + x * (dll(3) + x * (dll(4) & - + x * (dll(5) + x * (dll(6) + x * (dll(7) + x * dll(8)) ) ) ) ) ) ) - end if + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esl = eslf(temp) + eslfp = esl * alvl(temp) / (rh2o * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! @@ -764,7 +852,9 @@ end function eslfp ! with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=4) function esifp(temp) - use consts_coms, only : t00 ! ! intent(in) + use consts_coms , only : t00 & ! intent(in) + , rh2o ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=4), intent(in) :: temp @@ -777,19 +867,35 @@ real(kind=4) function esifp(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! - esi = esif(temp) - iiprime = -iii_7(1)/(temp*temp) + iii_7(2)/temp + iii_7(3) - esifp = esi * iiprime - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! x=max(-80.,temp-t00) esifp = dii(0) + x * (dii(1) + x * (dii(2) + x * (dii(3) + x * (dii(4) & + x * (dii(5) + x * (dii(6) + x * (dii(7) + x * dii(8)) ) ) ) ) ) ) - end if + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! + esi = esif(temp) + iiprime = -iii_7(1)/(temp*temp) + iii_7(2)/temp + iii_7(3) + esifp = esi * iiprime + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esi = esif(temp) + esifp = esi * alvi(temp) / (rh2o * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -2943,13 +3049,9 @@ end function alvi ! function [J/kg/K], temperature [K], and liquid and ice mixing ratios [kg/kg]. ! !---------------------------------------------------------------------------------------! real(kind=4) function theta_iceliq(exner,temp,rliq,rice) - use consts_coms, only : alvl3 & ! intent(in) - , alvi3 & ! intent(in) - , cpdry & ! intent(in) - , ttripoli & ! intent(in) - , htripoli & ! intent(in) - , htripolii ! ! intent(in) - + use consts_coms , only : alvl3 & ! intent(in) + , alvi3 & ! intent(in) + , cpdry ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! real(kind=4), intent(in) :: exner ! Exner function [ J/kg/K] @@ -2967,36 +3069,13 @@ real(kind=4) function theta_iceliq(exner,temp,rliq,rice) !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use the ! - ! latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - !------------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - else - qq = alvl3 * rliq + alvi3 * rice - end if + !----- Find the latent heat enthalpy. -----------------------------------------------! + qq = alvl(temp) * rliq + alvi(temp) * rice !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Solve the thermodynamics. For the new thermodynamics we don't approximate ! - ! the exponential to a linear function, nor do we impose temperature above the thre- ! - ! shold from Tripoli and Cotton (1981). ! - !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Decide how to compute, based on temperature. ------------------------------! - theta_iceliq = hh * exp(-qq / hh) / exner - !---------------------------------------------------------------------------------! - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli) then - theta_iceliq = hh * hh / (exner * ( hh + qq)) - else - theta_iceliq = hh * htripoli / (exner * ( htripoli + qq)) - end if - !---------------------------------------------------------------------------------! - end if + !----- Find theta_iceliq. -----------------------------------------------------------! + theta_iceliq = hh * exp(-qq / hh) / exner !------------------------------------------------------------------------------------! return @@ -3015,15 +3094,12 @@ end function theta_iceliq ! to temperature, useful in iterative methods. ! !---------------------------------------------------------------------------------------! real(kind=4) function dthetail_dt(condconst,thil,pres,temp,rliq,ricein) - use consts_coms, only : alvl3 & ! intent(in) - , alvi3 & ! intent(in) - , dcpvi & ! intent(in) - , dcpvl & ! intent(in) - , cpdry & ! intent(in) - , ttripoli & ! intent(in) - , htripoli & ! intent(in) - , htripolii & ! intent(in) - , t3ple ! ! intent(in) + use consts_coms , only : alvl3 & ! intent(in) + , alvi3 & ! intent(in) + , dcpvi & ! intent(in) + , dcpvl & ! intent(in) + , cpdry & ! intent(in) + , t3ple ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -3074,18 +3150,9 @@ real(kind=4) function dthetail_dt(condconst,thil,pres,temp,rliq,ricein) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use ! - ! the latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - ! The term r x d(L)/dT x T is computed only when we use the new thermodynamics. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - rdlt = (dcpvl * rliq + dcpvi * rice ) * temp - else - qq = alvl3 * rliq + alvi3 * rice - rdlt = 0.0 - end if + !----- Find the latent heat enthalpy. --------------------------------------------! + qq = alvl(temp) * rliq + alvi(temp) * rice + rdlt = (dcpvl * rliq + dcpvi * rice ) * temp !---------------------------------------------------------------------------------! @@ -3093,24 +3160,15 @@ real(kind=4) function dthetail_dt(condconst,thil,pres,temp,rliq,ricein) !---------------------------------------------------------------------------------! ! This is the term Lx[d(rs)/dt]xT. L may be either the vapourisation or ! ! sublimation latent heat, depending on the temperature and whether we are consi- ! - ! dering ice or not. We still need to check whether latent heat is a function of ! - ! temperature or not. Also, if condensation mixing ratio is constant, then this ! - ! term will be always zero. ! + ! dering ice or not. If condensation mixing ratio is constant, then this term ! + ! will be always zero. ! !---------------------------------------------------------------------------------! if (condconst) then ldrst = 0. elseif (thereisice .and. temp < t3ple) then - if (newthermo) then - ldrst = alvi3 * rsifp(pres,temp) * temp - else - ldrst = alvi(temp) * rsifp(pres,temp) * temp - end if + ldrst = alvi(temp) * rsifp(pres,temp) * temp else - if (newthermo) then - ldrst = alvl3 * rslfp(pres,temp) * temp - else - ldrst = alvl(temp) * rslfp(pres,temp) * temp - end if + ldrst = alvl(temp) * rslfp(pres,temp) * temp end if !---------------------------------------------------------------------------------! end if @@ -3121,17 +3179,7 @@ real(kind=4) function dthetail_dt(condconst,thil,pres,temp,rliq,ricein) !------------------------------------------------------------------------------------! ! Find the condensed phase consistent with the thermodynamics used. ! !------------------------------------------------------------------------------------! - if (newthermo) then - dthetail_dt = thil * ( 1. + (ldrst + qq - rdlt ) / hh ) / temp - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli) then - dthetail_dt = thil * ( 1. + (ldrst + qq) / (hh+qq) ) / temp - else - dthetail_dt = thil * ( 1. + ldrst / (htripoli + alvl3 * rliq) ) / temp - end if - !---------------------------------------------------------------------------------! - end if + dthetail_dt = thil * ( 1. + (ldrst + qq - rdlt ) / hh ) / temp !------------------------------------------------------------------------------------! return @@ -3148,9 +3196,6 @@ end function dthetail_dt !=======================================================================================! ! This function computes temperature from the ice-liquid water potential temperature ! ! in Kelvin, Exner function in J/kg/K, and liquid and ice mixing ratios in kg/kg. ! - ! For now t1stguess is used only to decide whether I should use the complete case or ! - ! the 253 K to reduce the error on neglecting the changes on latent heat due to temper- ! - ! ature. ! !---------------------------------------------------------------------------------------! real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) use consts_coms, only : cpdry & ! intent(in) @@ -3159,9 +3204,7 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) , alvl3 & ! intent(in) , alvi3 & ! intent(in) , t00 & ! intent(in) - , t3ple & ! intent(in) - , ttripoli & ! intent(in) - , htripolii ! ! intent(in) + , t3ple ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real(kind=4), intent(in) :: thil ! Ice-liquid water potential temp. [ K] @@ -3171,7 +3214,6 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) real(kind=4), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=4), intent(in) :: t1stguess ! 1st. guess for temperature [ K] !----- Local variables. -------------------------------------------------------------! - real(kind=4) :: til ! Ice liquid temperature [ K] real(kind=4) :: deriv ! Function derivative real(kind=4) :: fun ! Function for which we seek a root. real(kind=4) :: funa ! Smallest guess function @@ -3194,20 +3236,6 @@ real(kind=4) function thil2temp(thil,exner,pres,rliq,rice,t1stguess) thil2temp = cpdryi * thil * exner return !---------------------------------------------------------------------------------! - elseif (.not. newthermo) then - !---------------------------------------------------------------------------------! - ! There is condensation but we are using the old thermodynamics, which can be ! - ! solved analytically. ! - !---------------------------------------------------------------------------------! - til = cpdryi * thil * exner - if (t1stguess > ttripoli) then - thil2temp = 0.5 & - * (til + sqrt(til * (til + cpdryi4 * (alvl3 * rliq + alvi3 * rice)))) - else - thil2temp = til * ( 1. + (alvl3 * rliq + alvi3 * rice) * htripolii) - end if - return - !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -3411,34 +3439,25 @@ end function thil2temp ! saturation mixing ratio [kg/kg], keeping pressure constant. This is based on the ! ! ice-liquid potential temperature equation. ! !---------------------------------------------------------------------------------------! - real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) + real(kind=4) function dtempdrs(temp,rliq,rice,rconmin) use consts_coms, only : alvl3 & ! intent(in) , alvi3 & ! intent(in) , dcpvl & ! intent(in) , dcpvi & ! intent(in) , cpdry & ! intent(in) - , cpdryi & ! intent(in) - , ttripoli & ! intent(in) - , htripolii ! ! intent(in) + , cpdryi ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real(kind=4), intent(in) :: exner ! Exner function [ J/kg/K] - real(kind=4), intent(in) :: thil ! Ice-liquid potential temperature (*) [ K] real(kind=4), intent(in) :: temp ! Temperature [ K] real(kind=4), intent(in) :: rliq ! Liquid mixing ratio [ kg/kg] real(kind=4), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=4), intent(in) :: rconmin ! Min. non-zero condensate mix. ratio [ kg/kg] - !------------------------------------------------------------------------------------! - ! (*) Thil is not used in this formulation but it may be used should you opt by ! - ! other ways to compute theta_il, so don't remove this argument. ! - !------------------------------------------------------------------------------------! !----- Local variables --------------------------------------------------------------! real(kind=4) :: qq ! Enthalpy -- latent heat [ J/kg] real(kind=4) :: qpt ! d(qq)/dT * T [ J/kg] real(kind=4) :: hh ! Enthalpy -- sensible heat [ J/kg] real(kind=4) :: rcon ! Condensate mixing ratio [ kg/kg] - real(kind=4) :: til ! Ice-liquid temperature [ K] !------------------------------------------------------------------------------------! @@ -3457,15 +3476,10 @@ real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! ! Find the enthalpy associated with latent heat and its derivative ! - ! correction. This is dependent on the thermodynamics used. ! + ! correction. ! !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl(temp) * rliq + alvi(temp) * rice - qpt = dcpvl * rliq + dcpvi * rice - else - qq = alvl3 * rliq + alvi3 * rice - qpt = 0.0 - end if + qq = alvl(temp) * rliq + alvi(temp) * rice + qpt = dcpvl * rliq + dcpvi * rice !---------------------------------------------------------------------------------! @@ -3474,21 +3488,8 @@ real(kind=4) function dtempdrs(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Decide how to compute, based on the thermodynamics method. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - dtempdrs = - temp * qq / ( rcon * (hh + qq - qpt) ) - else - til = cpdryi * thil * exner - !----- Decide how to compute, based on temperature. ---------------------------! - if (temp > ttripoli) then - dtempdrs = - til * qq / ( rcon * cpdry * (2.*temp-til)) - else - dtempdrs = - til * qq * htripolii / rcon - end if - !------------------------------------------------------------------------------! - end if + !----- Find the temperature derivative wrt saturation mixing ratio. --------------! + dtempdrs = - temp * qq / ( rcon * (hh + qq - qpt) ) !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -4968,7 +4969,7 @@ subroutine thil2tqall(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) !write (unit=46,fmt='(a,1x,i5,1x,6(a,1x,f11.4,1x))') & ! 'ANSWER: it=',itb,'funf=',funnow,'temp=',temp-t00 & ! ,'rsat=',1000.*rsat,'rliq=',1000.*rliq,'rice=',1000.*rice,'rvap=',1000.*rvap - !write (unit=46,fmt='(a)') '----------------------------------------------------------' + !write (unit=46,fmt='(a)') '---------------------------------------------------------! !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>! !><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> maxfpo & ! intent(in) , maxit4 => maxit & ! intent(in) , maxlev4 => maxlev & ! intent(in) - , newthermo4 => newthermo & ! intent(in) , level4 => level & ! intent(in) , vapour_on4 => vapour_on & ! intent(in) , cloud_on4 => cloud_on & ! intent(in) @@ -39,7 +38,6 @@ module therm_lib8 ! crashing, for other methods. integer , parameter :: maxlev = maxlev4 ! Maximum # of levels for adaptive ! quadrature methods. - logical , parameter :: newthermo = newthermo4 ! Use new thermodynamics [T|F] !---------------------------------------------------------------------------------------! @@ -145,10 +143,18 @@ module therm_lib8 !=======================================================================================! !=======================================================================================! ! This function calculates the liquid saturation vapour pressure as a function of ! - ! Kelvin temperature. This expression came from MK05, equation (10). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function eslf8(temp,l1funout,l2funout,ttfunout) - use consts_coms, only : t008 ! ! intent(in) + use consts_coms , only : t008 & ! intent(in) + , es3ple8 & ! intent(in) + , t3ple8 & ! intent(in) + , dcplorv8 & ! intent(in) + , lvorvot38 & ! intent(in) + , rh2o8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=8), intent(in) :: temp ! Temperature [ K] @@ -166,30 +172,57 @@ real(kind=8) function eslf8(temp,l1funout,l2funout,ttfunout) !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !----- Original method, using polynomial fit (FWC92) -----------------------------! + x = max(-8.0d1,temp-t008) + eslf8 = cll8(0) + x * (cll8(1) + x * (cll8(2) + x * (cll8(3) + x * (cll8(4) & + + x * (cll8(5) + x * (cll8(6) + x * (cll8(7) + x * cll8(8)) )))))) + !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! + if (present(l1funout)) l1funout = eslf8 + if (present(l2funout)) l2funout = eslf8 + if (present(ttfunout)) ttfunout = eslf8 + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! l1fun = l01_108(0) + l01_108(1)/temp + l01_108(2)*log(temp) + l01_108(3) * temp l2fun = l02_108(0) + l02_108(1)/temp + l02_108(2)*log(temp) + l02_108(3) * temp ttfun = tanh(ttt_108(1) * (temp - ttt_108(2))) eslf8 = exp(l1fun + ttfun*l2fun) !---------------------------------------------------------------------------------! + !----- These are useful for the derivative calculation. --------------------------! if (present(l1funout)) l1funout = l1fun if (present(l2funout)) l2funout = l2fun if (present(ttfunout)) ttfunout = ttfun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.0d1,temp-t008) - eslf8 = cll8(0) + x * (cll8(1) + x * (cll8(2) + x * (cll8(3) + x * (cll8(4) & - + x * (cll8(5) + x * (cll8(6) + x * (cll8(7) + x * cll8(8)) )))))) + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + l1fun = min(lnexp_max8,max(lnexp_min8, lvorvot38 - alvl8(temp) / (rh2o8 * temp))) + eslf8 = es3ple8 * (t3ple8/temp) ** dcplorv8 * exp(l1fun) !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(l1funout)) l1funout = eslf8 if (present(l2funout)) l2funout = eslf8 if (present(ttfunout)) ttfunout = eslf8 - end if + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! return end function eslf8 @@ -204,10 +237,18 @@ end function eslf8 !=======================================================================================! !=======================================================================================! ! This function calculates the ice saturation vapour pressure as a function of ! - ! Kelvin temperature, based on MK05 equation (7). ! + ! Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function esif8(temp,iifunout) - use consts_coms, only : t008 ! ! intent(in) + use consts_coms , only : t008 & ! intent(in) + , es3ple8 & ! intent(in) + , t3ple8 & ! intent(in) + , dcpiorv8 & ! intent(in) + , lsorvot38 & ! intent(in) + , rh2o8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !----- Required arguments. ----------------------------------------------------------! real(kind=8), intent(in) :: temp ! Temperature [ K] @@ -219,27 +260,51 @@ real(kind=8) function esif8(temp,iifunout) !------------------------------------------------------------------------------------! + !------------------------------------------------------------------------------------! - ! Choose between the old and the new thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then + select case (esat_method) + case (0) + !----- Original method, using polynomial fit (FWC92) -----------------------------! + x = max(-8.d1,temp-t008) + esif8 = cii8(0) + x * (cii8(1) + x * (cii8(2) + x * (cii8(3) + x * (cii8(4) & + + x * (cii8(5) + x * (cii8(6) + x * (cii8(7) + x * cii8(8)))))))) + !---------------------------------------------------------------------------------! - !----- Updated method, using MK05 ------------------------------------------------! + !----- Dummy values for optional function in case they are requested. ------------! + if (present(iifunout)) iifunout=esif8 + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! iifun = iii_78(0) + iii_78(1)/temp + iii_78(2) * log(temp) + iii_78(3) * temp esif8 = exp(iifun) !---------------------------------------------------------------------------------! + !----- This is useful for the derivative calculation. ------==--------------------! if (present(iifunout)) iifunout=iifun - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.d1,temp-t008) - esif8 = cii8(0) + x * (cii8(1) + x * (cii8(2) + x * (cii8(3) + x * (cii8(4) & - + x * (cii8(5) + x * (cii8(6) + x * (cii8(7) + x * cii8(8)))))))) + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method, which is thermodynamically more consistent with the other ! + ! ED-2.2 assumptions. ! + ! ! + ! Ambaum, M. H. P., 2020: Accurate, simple equation for saturated vapour pressure ! + ! over water and ice. Quart. J. Royal Meteorol. Soc., 146 (733), 4252-4258, ! + ! doi:10.1002/qj.3899 (hereafter A20). ! + !---------------------------------------------------------------------------------! + iifun = min(lnexp_max8,max(lnexp_min8, lsorvot38 - alvi8(temp) / (rh2o8 * temp))) + esif8 = es3ple8 * (t3ple8/temp) ** dcpiorv8 * exp(iifun) !---------------------------------------------------------------------------------! + + !----- Dummy values for optional functions in case they are requested. -----------! if (present(iifunout)) iifunout=esif8 - end if + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -725,7 +790,9 @@ end function rhovsil8 ! pressure with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function eslfp8(temp) - use consts_coms, only : t008 ! ! intent(in) + use consts_coms , only : t008 & ! intent(in) + , rh2o8 ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=8), intent(in) :: temp @@ -742,24 +809,39 @@ real(kind=8) function eslfp8(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! + x = max(-8.d1,temp-t008) + eslfp8 = dll8(0) + x * (dll8(1) + x * (dll8(2) + x * (dll8(3) + x * (dll8(4) & + + x * (dll8(5) + x * (dll8(6) + x * (dll8(7) + x * dll8(8)))))))) + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 10). ! + !---------------------------------------------------------------------------------! esl = eslf8(temp,l2funout=l2fun,ttfunout=ttfun) l1prime = -l01_108(1)/(temp*temp) + l01_108(2)/temp + l01_108(3) l2prime = -l02_108(1)/(temp*temp) + l02_108(2)/temp + l02_108(3) ttprime = ttt_108(1)*(1.d0 - ttfun*ttfun) eslfp8 = esl * (l1prime + l2prime*ttfun + l2fun*ttprime) - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! - x = max(-8.d1,temp-t008) - eslfp8 = dll8(0) + x * (dll8(1) + x * (dll8(2) + x * (dll8(3) + x * (dll8(4) & - + x * (dll8(5) + x * (dll8(6) + x * (dll8(7) + x * dll8(8)))))))) - end if + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esl = eslf8(temp) + eslfp8 = esl * alvl8(temp) / (rh2o8 * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! - return end function eslfp8 !=======================================================================================! @@ -776,7 +858,9 @@ end function eslfp8 ! with respect to temperature as a function of Kelvin temperature. ! !---------------------------------------------------------------------------------------! real(kind=8) function esifp8(temp) - use consts_coms, only : t008 ! ! intent(in) + use consts_coms , only : t008 & ! intent(in) + , rh2o8 ! ! intent(in) + use ed_misc_coms, only : esat_method ! ! intent(in) implicit none !------ Arguments. ------------------------------------------------------------------! real(kind=8), intent(in) :: temp @@ -789,19 +873,35 @@ real(kind=8) function esifp8(temp) !------------------------------------------------------------------------------------! - ! Decide which function to use, based on the thermodynamics. ! + ! Choose the thermodynamic approach. ! !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Updated method, using MK05 ------------------------------------------------! - esi = esif8(temp) - iiprime = -iii_78(1)/(temp*temp) + iii_78(2)/temp + iii_78(3) - esifp8 = esi * iiprime - else - !----- Original method, using polynomial fit (FWC92) -----------------------------! + select case (esat_method) + case (0) + !---------------------------------------------------------------------------------! + ! Original ED2.0 method, polynomial fit based on FWC92. ! + !---------------------------------------------------------------------------------! x = max(-8.d1,temp-t008) esifp8 = dii8(0) + x * (dii8(1) + x * (dii8(2) + x * (dii8(3) + x * (dii8(4) & + x * (dii8(5) + x * (dii8(6) + x * (dii8(7) + x * dii8(8)))))))) - end if + !---------------------------------------------------------------------------------! + case (1) + !---------------------------------------------------------------------------------! + ! ED-2.2 method, use MK05 (equation 7). ! + !---------------------------------------------------------------------------------! + esi = esif8(temp) + iiprime = -iii_78(1)/(temp*temp) + iii_78(2)/temp + iii_78(3) + esifp8 = esi * iiprime + !---------------------------------------------------------------------------------! + case (2) + !---------------------------------------------------------------------------------! + ! Use A20's method. The derivative is by definition the Clausius-Clapeyron ! + ! equation. This is the key advantage of this approach as it makes the ! + ! formulation consistent with the second law of thermodynamics. ! + !---------------------------------------------------------------------------------! + esi = esif8(temp) + esifp8 = esi * alvi8(temp) / (rh2o8 * temp * temp) + !---------------------------------------------------------------------------------! + end select !------------------------------------------------------------------------------------! return @@ -2954,10 +3054,7 @@ end function alvi8 real(kind=8) function theta_iceliq8(exner,temp,rliq,rice) use consts_coms, only : alvl38 & ! intent(in) , alvi38 & ! intent(in) - , cpdry8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , cpdry8 ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -2976,36 +3073,13 @@ real(kind=8) function theta_iceliq8(exner,temp,rliq,rice) !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use the ! - ! latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - !------------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - else - qq = alvl38 * rliq + alvi38 * rice - end if + !----- Find the latent heat enthalpy. -----------------------------------------------! + qq = alvl8(temp) * rliq + alvi8(temp) * rice !------------------------------------------------------------------------------------! - !------------------------------------------------------------------------------------! - ! Solve the thermodynamics. For the new thermodynamics we don't approximate ! - ! the exponential to a linear function, nor do we impose temperature above the thre- ! - ! shold from Tripoli and Cotton (1981). ! - !------------------------------------------------------------------------------------! - if (newthermo) then - !----- Decide how to compute, based on temperature. ------------------------------! - theta_iceliq8 = hh * exp(-qq / hh) / exner - !---------------------------------------------------------------------------------! - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli8) then - theta_iceliq8 = hh * hh / (exner * ( hh + qq)) - else - theta_iceliq8 = hh * htripoli8 / (exner * ( htripoli8 + qq)) - end if - !---------------------------------------------------------------------------------! - end if + !----- Find theta_iceliq. -----------------------------------------------------------! + theta_iceliq8 = hh * exp(-qq / hh) / exner !------------------------------------------------------------------------------------! return @@ -3029,9 +3103,6 @@ real(kind=8) function dthetail_dt8(condconst,thil,pres,temp,rliq,ricein) , dcpvi8 & ! intent(in) , dcpvl8 & ! intent(in) , cpdry8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripoli8 & ! intent(in) - , htripolii8 & ! intent(in) , t3ple8 ! ! intent(in) implicit none @@ -3083,42 +3154,24 @@ real(kind=8) function dthetail_dt8(condconst,thil,pres,temp,rliq,ricein) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Find the latent heat enthalpy. If using the old thermodynamics, we use ! - ! the latent heat at T = T3ple, otherwise we use the temperature-dependent one. ! - ! The term r � d(L)/dT � T is computed only when we use the new thermodynamics. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - rdlt = (dcpvl8 * rliq + dcpvi8 * rice ) * temp - else - qq = alvl38 * rliq + alvi38 * rice - rdlt = 0.d0 - end if + !----- Find the latent heat enthalpy. --------------------------------------------! + qq = alvl8(temp) * rliq + alvi8(temp) * rice + rdlt = (dcpvl8 * rliq + dcpvi8 * rice ) * temp !---------------------------------------------------------------------------------! !---------------------------------------------------------------------------------! - ! This is the term L�[d(rs)/dt]�T. L may be either the vapourisation or ! + ! This is the term Lx[d(rs)/dt]xT. L may be either the vapourisation or ! ! sublimation latent heat, depending on the temperature and whether we are consi- ! - ! dering ice or not. We still need to check whether latent heat is a function of ! - ! temperature or not. Also, if condensation mixing ratio is constant, then this ! - ! term will be always zero. ! + ! dering ice or not. If condensation mixing ratio is constant, then this term ! + ! will be always zero. ! !---------------------------------------------------------------------------------! if (condconst) then ldrst = 0.d0 elseif (thereisice .and. temp < t3ple8) then - if (newthermo) then - ldrst = alvi38 * rsifp8(pres,temp) * temp - else - ldrst = alvi8(temp) * rsifp8(pres,temp) * temp - end if + ldrst = alvi8(temp) * rsifp8(pres,temp) * temp else - if (newthermo) then - ldrst = alvl38 * rslfp8(pres,temp) * temp - else - ldrst = alvl8(temp) * rslfp8(pres,temp) * temp - end if + ldrst = alvl8(temp) * rslfp8(pres,temp) * temp end if !---------------------------------------------------------------------------------! end if @@ -3129,17 +3182,7 @@ real(kind=8) function dthetail_dt8(condconst,thil,pres,temp,rliq,ricein) !------------------------------------------------------------------------------------! ! Find the condensed phase consistent with the thermodynamics used. ! !------------------------------------------------------------------------------------! - if (newthermo) then - dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq - rdlt ) / hh ) / temp - else - !----- Decide how to compute, based on temperature. ------------------------------! - if (temp > ttripoli8) then - dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq) / (hh+qq) ) / temp - else - dthetail_dt8 = thil * ( 1.d0 + ldrst / (htripoli8 + alvl38 * rliq) ) / temp - end if - !---------------------------------------------------------------------------------! - end if + dthetail_dt8 = thil * ( 1.d0 + (ldrst + qq - rdlt ) / hh ) / temp !------------------------------------------------------------------------------------! return @@ -3156,9 +3199,6 @@ end function dthetail_dt8 !=======================================================================================! ! This function computes temperature from the ice-liquid water potential temperature ! ! in Kelvin, Exner function in J/kg/K, and liquid and ice mixing ratios in kg/kg. ! - ! For now t1stguess is used only to decide whether I should use the complete case or ! - ! the 253 K to reduce the error on neglecting the changes on latent heat due to temper- ! - ! ature. ! !---------------------------------------------------------------------------------------! real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) use consts_coms, only : cpdry8 & ! intent(in) @@ -3167,9 +3207,7 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) , alvl38 & ! intent(in) , alvi38 & ! intent(in) , t008 & ! intent(in) - , t3ple8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , t3ple8 ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real(kind=8), intent(in) :: thil ! Ice-liquid water potential temp. [ K] @@ -3179,7 +3217,6 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) real(kind=8), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=8), intent(in) :: t1stguess ! 1st. guess for temperature [ K] !----- Local variables. -------------------------------------------------------------! - real(kind=8) :: til ! Ice liquid temperature [ K] real(kind=8) :: deriv ! Function derivative real(kind=8) :: fun ! Function for which we seek a root. real(kind=8) :: funa ! Smallest guess function @@ -3202,21 +3239,6 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) thil2temp8 = cpdryi8 * thil * exner return !---------------------------------------------------------------------------------! - elseif (.not. newthermo) then - !---------------------------------------------------------------------------------! - ! There is condensation but we are using the old thermodynamics, which can be ! - ! solved analytically. ! - !---------------------------------------------------------------------------------! - til = cpdryi8 * thil * exner - if (t1stguess > ttripoli8) then - thil2temp8 = 5.d-1 & - * (til + sqrt( til & - * (til + cpdryi48 * (alvl38 * rliq + alvi38 * rice)))) - else - thil2temp8 = til * ( 1.d0 + (alvl38 * rliq + alvi38 * rice) * htripolii8) - end if - return - !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -3386,11 +3408,11 @@ real(kind=8) function thil2temp8(thil,exner,pres,rliq,rice,t1stguess) write (unit=*,fmt='(a,1x,f12.4)' ) 'Exner [J/kg/K] =',exner write (unit=*,fmt='(a,1x,f12.4)' ) 'rliq [ g/kg] =',1.d3*rliq write (unit=*,fmt='(a,1x,f12.4)' ) 'rice [ g/kg] =',1.d3*rice - write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ �C] =',t1stguess-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 't1stguess [ degC] =',t1stguess-t008 write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome (downdraft values).' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ �C] =',tempa-t008 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ �C] =',tempz-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempa [ degC] =',tempa-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tempz [ degC] =',tempz-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',fun write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -3420,34 +3442,25 @@ end function thil2temp8 ! saturation mixing ratio [kg/kg], keeping pressure constant. This is based on the ! ! ice-liquid potential temperature equation. ! !---------------------------------------------------------------------------------------! - real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) + real(kind=8) function dtempdrs8(temp,rliq,rice,rconmin) use consts_coms, only : alvl38 & ! intent(in) , alvi38 & ! intent(in) , dcpvl8 & ! intent(in) , dcpvi8 & ! intent(in) , cpdry8 & ! intent(in) - , cpdryi8 & ! intent(in) - , ttripoli8 & ! intent(in) - , htripolii8 ! ! intent(in) + , cpdryi8 ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real(kind=8), intent(in) :: exner ! Exner function [ J/kg/K] - real(kind=8), intent(in) :: thil ! Ice-liquid potential temperature (*) [ K] real(kind=8), intent(in) :: temp ! Temperature [ K] real(kind=8), intent(in) :: rliq ! Liquid mixing ratio [ kg/kg] real(kind=8), intent(in) :: rice ! Ice mixing ratio [ kg/kg] real(kind=8), intent(in) :: rconmin ! Min. non-zero condensate mix. ratio [ kg/kg] - !------------------------------------------------------------------------------------! - ! (*) Thil is not used in this formulation but it may be used should you opt by ! - ! other ways to compute theta_il, so don't remove this argument. ! - !------------------------------------------------------------------------------------! !----- Local variables --------------------------------------------------------------! real(kind=8) :: qq ! Enthalpy -- latent heat [ J/kg] real(kind=8) :: qpt ! d(qq)/dT * T [ J/kg] real(kind=8) :: hh ! Enthalpy -- sensible heat [ J/kg] real(kind=8) :: rcon ! Condensate mixing ratio [ kg/kg] - real(kind=8) :: til ! Ice-liquid temperature [ K] !------------------------------------------------------------------------------------! @@ -3466,15 +3479,10 @@ real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! ! Find the enthalpy associated with latent heat and its derivative ! - ! correction. This is dependent on the thermodynamics used. ! + ! correction. ! !---------------------------------------------------------------------------------! - if (newthermo) then - qq = alvl8(temp) * rliq + alvi8(temp) * rice - qpt = dcpvl8 * rliq + dcpvi8 * rice - else - qq = alvl38 * rliq + alvi38 * rice - qpt = 0.d0 - end if + qq = alvl8(temp) * rliq + alvi8(temp) * rice + qpt = dcpvl8 * rliq + dcpvi8 * rice !---------------------------------------------------------------------------------! @@ -3483,21 +3491,8 @@ real(kind=8) function dtempdrs8(exner,thil,temp,rliq,rice,rconmin) !---------------------------------------------------------------------------------! - !---------------------------------------------------------------------------------! - ! Decide how to compute, based on the thermodynamics method. ! - !---------------------------------------------------------------------------------! - if (newthermo) then - dtempdrs8 = - temp * qq / ( rcon * (hh + qq - qpt) ) - else - til = cpdryi8 * thil * exner - !----- Decide how to compute, based on temperature. ---------------------------! - if (temp > ttripoli8) then - dtempdrs8 = - til * qq / ( rcon * cpdry8 * (2.*temp-til)) - else - dtempdrs8 = - til * qq * htripolii8 / rcon - end if - !------------------------------------------------------------------------------! - end if + !----- Find the temperature derivative wrt saturation mixing ratio. --------------! + dtempdrs8 = - temp * qq / ( rcon * (hh + qq - qpt) ) !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! @@ -4551,13 +4546,13 @@ subroutine lcl_il8(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a,1x,f12.4)' ) 'theta_il [ K] =',thil write (unit=*,fmt='(a,1x,f12.4)' ) 'Pressure [ hPa] =',1.d-2*pres - write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ �C] =',temp-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'Temperature [ degC] =',temp-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'rtot [ g/kg] =',1.d3*rtot write (unit=*,fmt='(a,1x,f12.4)' ) 'rvap [ g/kg] =',1.d3*rvap write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' Last iteration outcome.' - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ �C] =',tlcla-t008 - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ �C] =',tlclz-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcla [ degC] =',tlcla-t008 + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlclz [ degC] =',tlclz-t008 write (unit=*,fmt='(a,1x,f12.4)' ) 'fun [ K] =',funnow write (unit=*,fmt='(a,1x,f12.4)' ) 'funa [ K] =',funa write (unit=*,fmt='(a,1x,f12.4)' ) 'funz [ K] =',funz @@ -4565,7 +4560,7 @@ subroutine lcl_il8(thil,pres,temp,rtot,rvap,tlcl,plcl,dzlcl,useice) write (unit=*,fmt='(a,1x,es12.4)') 'toler8 [ ----] =',toler8 write (unit=*,fmt='(a,1x,es12.4)') 'error [ ----] =' & ,abs(tlclz-tlcla)/tlclz - write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ �C] =',tlcl + write (unit=*,fmt='(a,1x,f12.4)' ) 'tlcl [ degC] =',tlcl call fatal_error('TLCL didn''t converge, gave up!','lcl_il8','therm_lib8.f90') end if return @@ -4598,8 +4593,7 @@ subroutine thil2tqall8(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) , cpdryi8 & ! intent(in) , t008 & ! intent(in) , toodry8 & ! intent(in) - , t3ple8 & ! intent(in) - , ttripoli8 ! ! intent(in) + , t3ple8 ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! @@ -4957,12 +4951,12 @@ subroutine thil2tqall8(thil,exner,pres,rtot,rliq,rice,temp,rvap,rsat) write (unit=*,fmt='(a)') ' ' write (unit=*,fmt='(a)') ' -> Output: ' write (unit=*,fmt='(a,1x,i12)') ' ITERATIONS :',itb - write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ �C]:',temp-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMP [ degC]:',temp-t008 write (unit=*,fmt='(a,1x,f12.5)') ' RVAP [ g/kg]:',1.d3*rvap write (unit=*,fmt='(a,1x,f12.5)') ' RLIQ [ g/kg]:',1.d3*rliq write (unit=*,fmt='(a,1x,f12.5)') ' RICE [ g/kg]:',1.d3*rice - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ �C]:',tempa-t008 - write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ �C]:',tempz-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPA [ degC]:',tempa-t008 + write (unit=*,fmt='(a,1x,f12.5)') ' TEMPZ [ degC]:',tempz-t008 write (unit=*,fmt='(a,1x,es12.5)') ' FUNA [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' FUNZ [ K]:',funnow write (unit=*,fmt='(a,1x,es12.5)') ' DERIV [ ---]:',deriv diff --git a/EDTS/Templates/ED2IN-MAIN b/EDTS/Templates/ED2IN-MAIN index f5d62d32f..7671cfdb6 100644 --- a/EDTS/Templates/ED2IN-MAIN +++ b/EDTS/Templates/ED2IN-MAIN @@ -769,6 +769,28 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + NL%ESAT_METHOD = 2 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. (Deprecated) Euler step. The fastest, but it has only a ! diff --git a/EDTS/Templates/ED2IN-TEST b/EDTS/Templates/ED2IN-TEST index 5c7ac833d..4d33c15c2 100644 --- a/EDTS/Templates/ED2IN-TEST +++ b/EDTS/Templates/ED2IN-TEST @@ -778,6 +778,28 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + NL%ESAT_METHOD = 2 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. (Deprecated) Euler step. The fastest, but it has only a ! diff --git a/EDTS/Templates/ED2IN-tonzi b/EDTS/Templates/ED2IN-tonzi index f18e8e84c..aff3bf9eb 100644 --- a/EDTS/Templates/ED2IN-tonzi +++ b/EDTS/Templates/ED2IN-tonzi @@ -721,6 +721,28 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + NL%ESAT_METHOD = 2 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. Euler step. The fastest, but it doesn't estimate ! diff --git a/EDTS/Templates/ED2IN-tonzi.harvest b/EDTS/Templates/ED2IN-tonzi.harvest index def693a1c..eff01e961 100644 --- a/EDTS/Templates/ED2IN-tonzi.harvest +++ b/EDTS/Templates/ED2IN-tonzi.harvest @@ -721,6 +721,28 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ESAT_METHOD -- Which method to use for computing vapour pressure at saturation. ! + ! 0. (Deprecated) Original ED2 formulation, based on Flatau et al. ! + ! (1992). This is not recommended because the function is not ! + ! continuous at T = 0 degC. ! + ! 1. (ED-2.2 default) Based on Murphy and Koop (2005). This is a fairly ! + ! robust approach and it is continuous at the triple point. The main ! + ! downside is that it makes the estimate slightly inconsistent with ! + ! the second law of thermodynamics. ! + ! 2. (ED-2.2 alternative) Based on Ambaum (2020). This is intended to ! + ! become the default. It is accurate, and uses the same assumptions ! + ! as the definition of enthalpy and latent heat of vaporisation and ! + ! sublimation, so it is consistent with the second law of ! + ! thermodynamics. The only minor disadvantage is that the function ! + ! is continuous but not differentiable at the fusion point. ! + !---------------------------------------------------------------------------------------! + NL%ESAT_METHOD = 1 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! INTEGRATION_SCHEME -- The biophysics integration scheme. ! ! 0. Euler step. The fastest, but it doesn't estimate ! diff --git a/EDTS/Templates/ED2IN-umbs.bg b/EDTS/Templates/ED2IN-umbs.bg index 58209b5af..6c13a7b21 100644 --- a/EDTS/Templates/ED2IN-umbs.bg +++ b/EDTS/Templates/ED2IN-umbs.bg @@ -103,6 +103,7 @@ $ED_NL NL%SLDRAIN = 10 NL%IVEGT_DYNAMICS = 1 NL%IBIGLEAF = 0 + NL%ESAT_METHOD = 2 NL%INTEGRATION_SCHEME = 0 NL%NSUB_EULER=40 NL%RK4_TOLERANCE = 0.01