subroutine kcmprof_fn(psurf_rcm,qsurf_rcm,Tsurf_rcm,Tstra_rcm,P_rcm,Pl_rcm,z_rcm,T_rcm,q_rcm,m_rcm) use params_h use watercommon_h, only : mH2O use gases_h implicit none ! ---------------------------------------------------------------- ! Purpose: create profiles of T, rho_v, rho_n, Pv and Pn following ! Kasting 1988 ! Authour: Adapted from a code by E. Marcq by R. Wordsworth (2011) ! ---------------------------------------------------------------- #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" integer ilay, nlay parameter (nlay=10000) ! number of vertical layers ! rcm inputs real Tsurf_rcm,Tstra_rcm ! rcm outputs real psurf_rcm,qsurf_rcm real P_rcm(1:nlayermx) real Pl_rcm(1:nlayermx+1) real z_rcm(1:nlayermx) real T_rcm(1:nlayermx),q_rcm(1:nlayermx) real m_rcm(1:nlayermx+1) ! rcm for interpolation (should really use log coords?) !double precision p1,p2,pnew,ilay_rcm double precision lnp1,lnp2,lnpnew real Dp_rcm, dlogp_rcm integer ilay_rcm,ilev_rcm,ifinal_rcm double precision Dz, Dp double precision Ptop, dlogp, Psat_max parameter (Ptop=1.0) ! Pressure at TOA [Pa] !parameter (Psat_max=100000.0) ! Maximum vapour pressure [Pa] parameter (Psat_max=0.0) ! set to zero for dry calculations double precision T(1:nlay) ! temperature [K] double precision Ztab(1:nlay) ! altitude [m] double precision Pv(1:nlay),Pn(1:nlay),P(1:nlay) ! pressure [Pa] double precision rho_v(1:nlay), rho_n(1:nlay) ! density [kg m^-3] double precision a_v(1:nlay) ! = rho_v/rho_n [kg/kg] double precision q_v(1:nlay) ! = rho_v/rho_tot [kg/kg] double precision mtot(1:nlay) ! = (rho_v+rho_n)/(n_v+n_n) [g/mol] integer profil_flag(1:nlay) ! 0 = dry, 1 = moist, 2 = isothermal ! inputs double precision Tsurf ! surface temperature [K] double precision Psurf_v ! surface par. pressure (variable species) [Pa] double precision Psurf_n ! surface par. pressure (incondensible species)[Pa] double precision Ttop ! stratospheric temperature [K] double precision dTdp ! [K/Pa] double precision dPvdp,dPndp ! [Pa/Pa] double precision psat_v ! local Psat_H2O value double precision Tcrit ! Critical temperature [K] ! double precision Tsat ! local Tsat_H2O value double precision rho_vTEMP,rho_nTEMP ! variables necessary for steam.f90 double precision rhol,rhov,nul ! for output double precision vmr logical verbose parameter(verbose=.false.) ! initialise flags profil_flag(:) = 0 !------------------------------- ! assign input variables m_n = dble(mugaz/1000.) ! modify/generalise later?? if(ngasmx.gt.2)then print*,'Only two species possible at present' stop elseif(ngasmx.eq.1)then if(psat_max.gt.0.0)then print*,'Must have Psat_max=0 if no variable species' stop endif print*, 'Assuming pure atmosphere' m_v = 1.0 tcrit = 1000.0 elseif(gnom(ngasmx).eq.'H2O')then m_v = dble(mH2O/1000.) tcrit = 6.47d2 elseif(gnom(ngasmx).eq.'NH3')then m_v = 17.031/1000. tcrit = 4.06d2 elseif(gnom(ngasmx).eq.'CH4')then m_v = 16.04/1000. tcrit = 1.91d2 stop else print*,'Variable gas not recognised!' call abort endif rmn = rc/m_n Psurf_n = dble(psurf_rcm) Ttop = dble(Tstra_rcm) Tsurf = dble(Tsurf_rcm) ! assume vapour saturation at surface (for now) ! if (tsurf > tcrit) then ! print*,'Above critical temperature for ',gnom(2),& ! ', at surface, assuming 10 bar pressure instead' ! Psurf_v = 10*1d5 ! profil_flag(1) = 0 !i! endif psat_v=psat_max if(gnom(ngasmx).eq.'H2O')then call Psat_H2O(tsurf,psat_v) elseif(gnom(ngasmx).eq.'NH3')then call Psat_NH3(tsurf,psat_v) endif ! Moist adiabat unless greater than or equal to psat_max if(psat_v*1d6.lt.psat_max)then Psurf_v = Psat_v*1d6 profil_flag(1) = 1 else Psurf_v = psat_max profil_flag(1) = 0 endif if(verbose)then print*,'Psat_v =',psat_v*1d6 print*,'Tsurf =',Tsurf,' K' print*,'Ttop =',Ttop,' K' print*,'Psurf_v =',Psurf_v,' Pa' print*,'Psurf_n =',Psurf_n,' Pa' print*,'m_n =',m_n,' kg/mol' print*,'m_v =',m_v,' kg/mol' print*,'rc =',rc endif ! define fine pressure grid psurf_rcm = real(Psurf_n+Psurf_v) dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayermx P_rcm(1) = psurf_rcm*exp(dlogp_rcm) do ilay_rcm=1,nlayermx-1 P_rcm(ilay_rcm+1) = P_rcm(ilay_rcm)*exp(dlogp_rcm) enddo Pl_rcm(1) = psurf_rcm do ilev_rcm=2,nlayermx ! log-linear interpolation Pl_rcm(ilev_rcm) = exp( log( P_rcm(ilev_rcm)*P_rcm(ilev_rcm-1) )/2 ) enddo !------------------------------- ! Layer 1 T(1) = Tsurf Pv(1) = Psurf_v Pn(1) = Psurf_n rho_n(1) = m_n*Pn(1)/(Rc*Tsurf) rho_v(1) = m_v*Pv(1)/(Rc*Tsurf) a_v(1) = rho_v(1)/rho_n(1) ! log pressure grid spacing (constant) dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1) call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp) if(verbose)then print*, 'dT/dp ground [K/Pa] =',dTdp endif ! initial delta p, delta z Dp = (Pn(1) + Pv(1))*(exp(dlogp) - 1d0) Dz = -Dp/( g*(rho_n(1) + rho_v(1)) ) !------------------------------- ! Layer 2 T(2) = tsurf + dTdp*Dp Pv(2) = Pv(1) + dPvdp*Dp Pn(2) = Pn(1) + dPndp*Dp rho_n(2) = m_n*Pn(2)/(Rc*T(2)) rho_v(2) = m_v*Pv(2)/(Rc*T(2)) a_v(2) = rho_v(2)/rho_n(2) !------------------------------- ! start vertical ascent Ztab(1) = 0. do ilay=2,nlay-1 ! calculate altitude levels (for diagnostic only) Dz = -Dp/( g*(rho_n(ilay) + rho_v(ilay)) ) Ztab(ilay) = Dz + Ztab(ilay-1) ! 1st assume next layer same as last one profil_flag(ilay) = profil_flag(ilay-1) ! update delta p Dp = (Pn(ilay)+Pv(ilay))*(exp(dlogp) - 1d0) ! intial gradients call to calculate temperature at next level call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),& T(ilay),dTdp,dPvdp,dPndp) T(ilay+1) = T(ilay) + dTdp*Dp ! test for moist adiabat at next level psat_v=psat_max if(gnom(ngasmx).eq.'H2O')then call Psat_H2O(T(ilay+1),psat_v) elseif(gnom(ngasmx).eq.'NH3')then call Psat_NH3(T(ilay+1),psat_v) endif if (psat_v*1d6 .lt. Pv(ilay)+dPvdp*Dp) then profil_flag(ilay)=1 call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),& T(ilay),dTdp,dPvdp,dPndp) endif ! test for stratosphere at next level if (T(ilay+1) .le. Ttop) then profil_flag(ilay)=2 T(ilay+1)=Ttop endif ! calculate pressures at next level Pn(ilay+1) = Pn(ilay) + dPndp*Dp Pv(ilay+1) = Pv(ilay) + dPvdp*Dp if(profil_flag(ilay) .eq. 1)then psat_v=psat_max if(gnom(ngasmx).eq.'H2O')then call Psat_H2O(T(ilay+1),psat_v) elseif(gnom(ngasmx).eq.'NH3')then call Psat_NH3(T(ilay+1),psat_v) endif if(Pv(ilay+1) .lt. psat_v*1e6)then Pv(ilay+1)=psat_v*1d6 endif endif ! calculate gas densities at next level (assume ideal) rho_n(ilay+1) = m_n*Pn(ilay+1)/(rc*T(ilay+1)) select case(profil_flag(ilay)) case(2) ! isothermal rho_v(ilay+1) = rho_v(ilay)/rho_n(ilay)*rho_n(ilay+1) case(1) ! moist ! dont think this is necessary !call psat_est(T(ilay+1),psat_v) ! modify for ammonia!!! rho_v(ilay+1) = m_v*psat_v*1d6/(rc*T(ilay+1)) case(0) ! dry rho_v(ilay+1) = m_v*Pv(ilay+1)/(rc*T(ilay+1)) end select ! print*,'profil_flag=',profil_flag(ilay) ! print*,'T=',T(ilay) ! print*,'a=',rho_v(ilay)/rho_n(ilay) ! if(profil_flag(ilay).eq.1)then ! stop ! endif enddo Ztab(nlay)=Ztab(nlay-1)+Dz !------------------------------- ! save to kcm1d variables ! surface quantities psurf_rcm = Pn(1) + Pv(1) qsurf_rcm = rho_v(1)/(rho_v(1) + rho_n(1)) ! create q_v, mtot for saving do ilay=1,nlay mtot(ilay) = 1d3*(rho_v(ilay) + rho_n(ilay)) / & (rho_v(ilay)/m_v + rho_n(ilay)/m_n) q_v(ilay) = rho_v(ilay)/(rho_v(ilay) + rho_n(ilay)) ! CHECK THIS enddo ! convert to rcm lower-res grid z_rcm(:) = 0.0 T_rcm(:) = 0.0 q_rcm(:) = 0.0 m_rcm(:) = 0.0 m_rcm(1) = real( 1d3*(rho_v(1) + rho_n(1)) / & (rho_v(1)/m_v + rho_n(1)/m_n) ) ilay_rcm=1 do ilay=2,nlay if(ilay_rcm.le.nlayermx)then ! interpolate rcm variables if(Pn(ilay)+Pv(ilay) .lt. P_rcm(ilay_rcm))then if(ilay.eq.1)then print*,'Error in create_profils: Psurf here less than Psurf in RCM!' call abort endif lnp1 = log(Pn(ilay-1)+Pv(ilay-1)) lnp2 = log(Pn(ilay)+Pv(ilay)) lnpnew = dble(log(P_rcm(ilay_rcm))) z_rcm(ilay_rcm) = real(Ztab(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) & + Ztab(ilay)*(lnpnew-lnp1)/(lnp2-lnp1)) T_rcm(ilay_rcm) = real(T(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) & + T(ilay)*(lnpnew-lnp1)/(lnp2-lnp1)) q_rcm(ilay_rcm) = real(q_v(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) & + q_v(ilay)*(lnpnew-lnp1)/(lnp2-lnp1)) m_rcm(ilay_rcm+1) = real(mtot(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) & + mtot(ilay)*(lnpnew-lnp1)/(lnp2-lnp1)) ilay_rcm = ilay_rcm+1 endif endif enddo ifinal_rcm=ilay_rcm-1 if(ifinal_rcm.lt.nlayermx)then if(verbose)then print*,'Interpolation in kcmprof stopped at layer',ilay_rcm,'!' endif do ilay_rcm=ifinal_rcm+1,nlayermx z_rcm(ilay_rcm) = z_rcm(ilay_rcm-1) T_rcm(ilay_rcm) = T_rcm(ilay_rcm-1) q_rcm(ilay_rcm) = q_rcm(ilay_rcm-1) m_rcm(ilay_rcm+1) = m_rcm(ilay_rcm) enddo endif do ilay=2,nlayermx if(T_rcm(ilay).lt.Ttop)then T_rcm(ilay)=Ttop endif enddo return end subroutine kcmprof_fn