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

Last change on this file since 364 was 305, checked in by rwordsworth, 14 years ago

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