source: trunk/LMDZ.GENERIC/libf/phystd/kcmprof_fn.F90 @ 848

Last change on this file since 848 was 716, checked in by rwordsworth, 12 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

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