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

Last change on this file since 374 was 366, checked in by rwordsworth, 13 years ago

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