source: trunk/LMDZ.PLUTO/libf/phypluto/dyn1d/kcmprof_fn.F90 @ 3572

Last change on this file since 3572 was 3232, checked in by afalco, 12 months ago

Pluto PCM:
1D functional
AF
  • Property svn:executable set to *
File size: 10.0 KB
RevLine 
[3184]1subroutine kcmprof_fn(nlayer,psurf_rcm,qsurf_rcm,Tsurf_rcm,Tstra_rcm,P_rcm,Pl_rcm,z_rcm,T_rcm,q_rcm,m_rcm)
2
3use params_h
4use gases_h
5use comcstfi_mod, only: mugaz, cpp, g
[3232]6use callkeys_mod, only: n2cond
[3184]7implicit none
8
9!     ----------------------------------------------------------------
10!     Purpose: create profiles of T, rho_v, rho_n, Pv and Pn following
11!              Kasting 1988
12!     Authour: Adapted from a code by E. Marcq by R. Wordsworth (2011)
13!     ----------------------------------------------------------------
14
15  integer ilay, nlay
16  parameter (nlay=10000) ! number of vertical layers
17
18  ! rcm inputs
19  integer nlayer
20  real Tsurf_rcm,Tstra_rcm
21
22  ! rcm outputs
23  real psurf_rcm,qsurf_rcm
24  real P_rcm(1:nlayer)
25  real Pl_rcm(1:nlayer+1)
26  real z_rcm(1:nlayer)
27  real T_rcm(1:nlayer),q_rcm(1:nlayer)
28  real m_rcm(1:nlayer+1)
29
30  ! rcm for interpolation (should really use log coords?)
31  !double precision p1,p2,pnew,ilay_rcm
32
33  double precision lnp1,lnp2,lnpnew
34  real Dp_rcm, dlogp_rcm
35  integer ilay_rcm,ilev_rcm,ifinal_rcm
36
37  double precision Dz, Dp
38  double precision Ptop, dlogp, Psat_max
39  parameter (Ptop=1.0)                             ! Pressure at TOA [Pa]
40
41  double precision T(1:nlay)                       ! temperature [K]
42  double precision Ztab(1:nlay)                    ! altitude [m]
43  double precision Pv(1:nlay),Pn(1:nlay),P(1:nlay) ! pressure [Pa]
44  double precision rho_v(1:nlay), rho_n(1:nlay)    ! density [kg m^-3]
45  double precision a_v(1:nlay)                     ! = rho_v/rho_n [kg/kg]
46  double precision q_v(1:nlay)                     ! = rho_v/rho_tot [kg/kg]
47  double precision mtot(1:nlay)                    ! = (rho_v+rho_n)/(n_v+n_n) [g/mol]
48
49  integer profil_flag(1:nlay) ! 0 = dry, 1 = moist, 2 = isothermal
50
51  ! inputs
52  double precision Tsurf   ! surface temperature [K]
53  double precision Psurf_v ! surface par. pressure (variable species) [Pa]
54  double precision Psurf_n ! surface par. pressure (incondensible species)[Pa]
55  double precision Ttop    ! stratospheric temperature [K]
56
57  double precision dTdp        ! [K/Pa]
58  double precision dPvdp,dPndp ! [Pa/Pa]
59  double precision psat_v      ! local Psat_H2O value
60  double precision Tcrit       ! Critical temperature [K]
61  double precision rho_vTEMP,rho_nTEMP
62
63  double precision TCO2cond ! for CO2 condensation quasi-hack
64
65  ! variables necessary for steam.f90
66  double precision rhol,rhov,nul
67
68  ! for output
69  double precision vmr
70
71  logical verbose
72  parameter(verbose=.true.)
73
74  logical add_Pvar_to_total
75  parameter(add_Pvar_to_total=.true.)
76
77  ! initialise flags
78  profil_flag(:) = 0
79
80  !-------------------------------
81  ! assign input variables
82  m_n     = dble(mugaz/1000.)
83  cp_n    = cpp
84  ! modify/generalise later??
85
86  Psat_max = 1000000.0 ! maximum vapour pressure [Pa]
87                       ! set huge until further notice
88
89  if(vgas.lt.1)then
90     if(psat_max.gt.0.0)then
91        print*,'Must have Psat_max=0 if no variable species'
92        psat_max=0.0
93        !stop
94     endif
95     print*, 'Assuming pure atmosphere'
96     m_v   = 1.0
97     tcrit = 1000.0
[3232]98!   elseif(trim(gnom(vgas)).eq.'H2O')then
99   !   m_v   = dble(mH2O/1000.) ! AF24: no water in Pluton
100   !   tcrit = 6.47d2
[3184]101  elseif(trim(gnom(vgas)).eq.'NH3')then
102     m_v   = 17.031/1000.
103     tcrit = 4.06d2
104  elseif(trim(gnom(vgas)).eq.'CH4')then
105     m_v   = 16.04/1000.
106     tcrit = 1.91d2
107     stop
108  else
109     print*,'Variable gas not recognised!'
110     call abort
111  endif
112
113  rmn     = rc/m_n
114  Ttop    = dble(Tstra_rcm)
115  Tsurf   = dble(Tsurf_rcm)
116
117  psat_v  = psat_max
118  if(vgas.gt.0)then
119     if(trim(gnom(vgas)).eq.'H2O')then
120        call Psat_H2O(tsurf,psat_v)
121     elseif(trim(gnom(vgas)).eq.'NH3')then
122        call Psat_NH3(tsurf,psat_v)
123     endif
124  endif
125
126  ! Moist adiabat unless greater than or equal to psat_max
127  if(psat_v*1d6.lt.psat_max)then
128    Psurf_v = Psat_v*1d6
129    profil_flag(1) = 1
130  else
131    Psurf_v = psat_max
132    profil_flag(1) = 0
133  endif
134
135  if(add_Pvar_to_total)then
136    Psurf_n = dble(psurf_rcm)
137    psurf_rcm = real(Psurf_n+Psurf_v)
138  else
139    Psurf_n = dble(psurf_rcm) - Psurf_v
140  endif
141
142  ! include relative humidity option
143  !if(satval.lt.1.0)then
144  !   Psurf_v = Psurf_v*satval
145  !   profil_flag(1) = 0     
146  !endif
147
148  if(verbose)then
149     print*,'Psat_v  =',psat_v*1d6
150     print*,'Tsurf   =',Tsurf,' K'
151     print*,'Ttop    =',Ttop,' K'
152     print*,'Psurf_v =',Psurf_v,' Pa'
153     print*,'Psurf_n =',Psurf_n,' Pa'
154     print*,'m_n     =',m_n,' kg/mol'
155     print*,'m_v     =',m_v,' kg/mol'
156     print*,'rc      =',rc
157  endif
158
159  ! define fine pressure grid
160  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayer
161
162  P_rcm(1)  = psurf_rcm*exp(dlogp_rcm)
163  do ilay_rcm=1,nlayer-1
164     P_rcm(ilay_rcm+1) = P_rcm(ilay_rcm)*exp(dlogp_rcm)
165  enddo
166
167  Pl_rcm(1) = psurf_rcm
168  do ilev_rcm=2,nlayer
169     ! log-linear interpolation
170     Pl_rcm(ilev_rcm) = exp( log( P_rcm(ilev_rcm)*P_rcm(ilev_rcm-1) )/2 )
171  enddo
172
173  !-------------------------------
174  ! Layer 1
175  T(1)     = Tsurf
176  Pv(1)    = Psurf_v
177  Pn(1)    = Psurf_n
178  rho_n(1) = m_n*Pn(1)/(Rc*Tsurf)
179  rho_v(1) = m_v*Pv(1)/(Rc*Tsurf)
180  a_v(1)   = rho_v(1)/rho_n(1)
181
182  ! log pressure grid spacing (constant)
183  dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1)
184 
185  call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp)
186  if(verbose)then
187     print*, 'dT/dp ground [K/Pa] =',dTdp
188  endif
189
190  ! initial delta p, delta z
191  Dp = (Pn(1) + Pv(1))*(exp(dlogp) - 1d0)
192  Dz = -Dp/(  g*(rho_n(1) + rho_v(1))  )
193
194  !-------------------------------
195  ! Layer 2
196  T(2)     = tsurf + dTdp*Dp
197  Pv(2)    = Pv(1) + dPvdp*Dp
198  Pn(2)    = Pn(1) + dPndp*Dp
199  rho_n(2) = m_n*Pn(2)/(Rc*T(2))
200  rho_v(2) = m_v*Pv(2)/(Rc*T(2))
201  a_v(2)   = rho_v(2)/rho_n(2)
202
203  !-------------------------------
204  ! start vertical ascent
205  Ztab(1) = 0.
206  do ilay=2,nlay-1
207
208     ! calculate altitude levels (for diagnostic only)
209     Dz         = -Dp/(  g*(rho_n(ilay) + rho_v(ilay))  )
210     Ztab(ilay) = Dz + Ztab(ilay-1)
211
212     ! 1st assume next layer same as last one
213     profil_flag(ilay) = profil_flag(ilay-1)
214 
215     ! update delta p
216     Dp = (Pn(ilay)+Pv(ilay))*(exp(dlogp) - 1d0)
217
218     ! intial gradients call to calculate temperature at next level
219     call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),&
220                        T(ilay),dTdp,dPvdp,dPndp)
221
222     T(ilay+1) = T(ilay) + dTdp*Dp
223
224     ! test for moist adiabat at next level
225     psat_v=psat_max
226
227     if(vgas.gt.0)then
228     if(trim(gnom(vgas)).eq.'H2O')then
229        call Psat_H2O(T(ilay+1),psat_v)
230     elseif(trim(gnom(vgas)).eq.'NH3')then
231        call Psat_NH3(T(ilay+1),psat_v)
232     endif
233     endif
234
235     if (psat_v*1d6 .lt. Pv(ilay)+dPvdp*Dp) then
236        profil_flag(ilay)=1
237        call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),&
238                         T(ilay),dTdp,dPvdp,dPndp)
239     endif
240
241     ! test for stratosphere at next level
242     if (T(ilay+1) .le. Ttop) then
243        profil_flag(ilay)=2
244        T(ilay+1)=Ttop
245     endif
246
247     ! calculate pressures at next level
248     Pn(ilay+1) = Pn(ilay) + dPndp*Dp
249     Pv(ilay+1) = Pv(ilay) + dPvdp*Dp
250
251     if(profil_flag(ilay) .eq. 1)then
252       
253        psat_v=psat_max
254
255        if(vgas.gt.0)then
256        if(trim(gnom(vgas)).eq.'H2O')then
257           call Psat_H2O(T(ilay+1),psat_v)
258        elseif(trim(gnom(vgas)).eq.'NH3')then
259           call Psat_NH3(T(ilay+1),psat_v)
260        endif
261        endif
262
263        if(Pv(ilay+1) .lt. psat_v*1e6)then
264           Pv(ilay+1)=psat_v*1d6
265        endif
266
267     endif
268
269     ! calculate gas densities at next level (assume ideal)
270     rho_n(ilay+1) = m_n*Pn(ilay+1)/(rc*T(ilay+1))
271     select case(profil_flag(ilay))
272     case(2) ! isothermal
273        rho_v(ilay+1) = rho_v(ilay)/rho_n(ilay)*rho_n(ilay+1)
274     case(1) ! moist
275
276        ! dont think this is necessary
277        !call psat_est(T(ilay+1),psat_v)
278        ! modify for ammonia!!!
279
280        rho_v(ilay+1) = m_v*psat_v*1d6/(rc*T(ilay+1))
281     case(0) ! dry
282        rho_v(ilay+1) = m_v*Pv(ilay+1)/(rc*T(ilay+1))
283     end select
284
285  enddo
286
287  Ztab(nlay)=Ztab(nlay-1)+Dz
288
289  !-------------------------------
290  ! save to kcm1d variables
291
292  ! surface quantities
293  psurf_rcm = Pn(1) + Pv(1)
294  qsurf_rcm = rho_v(1)/(rho_v(1) + rho_n(1))
295
296  ! create q_v, mtot for saving
297  do ilay=1,nlay
298     mtot(ilay) = 1d3*(rho_v(ilay) + rho_n(ilay)) / &
299                  (rho_v(ilay)/m_v + rho_n(ilay)/m_n)
300     q_v(ilay)  = rho_v(ilay)/(rho_v(ilay) + rho_n(ilay))
301     ! CHECK THIS
302  enddo
303
304
305  ! convert to rcm lower-res grid
306  z_rcm(:) = 0.0
307  T_rcm(:) = 0.0
308  q_rcm(:) = 0.0
309  m_rcm(:) = 0.0
310
311  m_rcm(1) = real( 1d3*(rho_v(1) + rho_n(1)) / &
312                  (rho_v(1)/m_v + rho_n(1)/m_n) )
313
314  ilay_rcm=1
315  do ilay=2,nlay
316
317     if(ilay_rcm.le.nlayer)then
318     ! interpolate rcm variables
319
320        if(Pn(ilay)+Pv(ilay) .lt. P_rcm(ilay_rcm))then
321
322           if(ilay.eq.1)then
323              print*,'Error in create_profils: Psurf here less than Psurf in RCM!'
324              call abort
325           endif
326
327           lnp1   = log(Pn(ilay-1)+Pv(ilay-1))
328           lnp2   = log(Pn(ilay)+Pv(ilay))
329           lnpnew = dble(log(P_rcm(ilay_rcm)))
330
331           z_rcm(ilay_rcm) = real(Ztab(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
332                             + Ztab(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
333           T_rcm(ilay_rcm) = real(T(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
334                             + T(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
335           q_rcm(ilay_rcm) = real(q_v(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
336                             + q_v(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
337
338           m_rcm(ilay_rcm+1) = real(mtot(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
339                             + mtot(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
340
341           ilay_rcm = ilay_rcm+1
342        endif
343
344     endif
345  enddo
346
347  ifinal_rcm=ilay_rcm-1
348  if(ifinal_rcm.lt.nlayer)then
349     if(verbose)then
350        print*,'Interpolation in kcmprof stopped at layer',ilay_rcm,'!'
351     endif
352
353     do ilay_rcm=ifinal_rcm+1,nlayer
354
355        z_rcm(ilay_rcm) = z_rcm(ilay_rcm-1)
356        T_rcm(ilay_rcm) = T_rcm(ilay_rcm-1)
357        q_rcm(ilay_rcm) = q_rcm(ilay_rcm-1)
358        m_rcm(ilay_rcm+1) = m_rcm(ilay_rcm)
359
360     enddo
361  endif
362
363  do ilay=2,nlayer
364     if(T_rcm(ilay).lt.Ttop)then
365        T_rcm(ilay)=Ttop
366     endif
367  enddo
368
369  return
370end subroutine kcmprof_fn
Note: See TracBrowser for help on using the repository browser.