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

Last change on this file since 1351 was 1308, checked in by emillour, 10 years ago

Generic GCM:
Some cleanup to simplify dynamics/physics interactions by getting rid
of dimphys.h (i.e. the nlayermx parameter) and minimizing use of
dimension.h in the physics.
EM
  • Property svn:executable set to *
File size: 10.7 KB
Line 
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 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  integer nlayer
24  real Tsurf_rcm,Tstra_rcm
25
26  ! rcm outputs
27  real psurf_rcm,qsurf_rcm
28  real P_rcm(1:nlayer)
29  real Pl_rcm(1:nlayer+1)
30  real z_rcm(1:nlayer)
31  real T_rcm(1:nlayer),q_rcm(1:nlayer)
32  real m_rcm(1:nlayer+1)
33
34  ! rcm for interpolation (should really use log coords?)
35  !double precision p1,p2,pnew,ilay_rcm
36
37  double precision lnp1,lnp2,lnpnew
38  real Dp_rcm, dlogp_rcm
39  integer ilay_rcm,ilev_rcm,ifinal_rcm
40
41  double precision Dz, Dp
42  double precision Ptop, dlogp, Psat_max
43  parameter (Ptop=1.0)                             ! Pressure at TOA [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 rho_vTEMP,rho_nTEMP
66
67  double precision TCO2cond ! for CO2 condensation quasi-hack
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=.true.)
77
78  logical add_Pvar_to_total
79  parameter(add_Pvar_to_total=.true.)
80
81  ! initialise flags
82  profil_flag(:) = 0
83
84  !-------------------------------
85  ! assign input variables
86  m_n     = dble(mugaz/1000.)
87  cp_n    = cpp
88  ! modify/generalise later??
89
90  Psat_max = 1000000.0 ! maximum vapour pressure [Pa]
91                       ! set huge until further notice
92
93  if(vgas.lt.1)then
94     if(psat_max.gt.0.0)then
95        print*,'Must have Psat_max=0 if no variable species'
96        psat_max=0.0
97        !stop
98     endif
99     print*, 'Assuming pure atmosphere'
100     m_v   = 1.0
101     tcrit = 1000.0
102  elseif(trim(gnom(vgas)).eq.'H2O')then
103     m_v   = dble(mH2O/1000.)
104     tcrit = 6.47d2
105  elseif(trim(gnom(vgas)).eq.'NH3')then
106     m_v   = 17.031/1000.
107     tcrit = 4.06d2
108  elseif(trim(gnom(vgas)).eq.'CH4')then
109     m_v   = 16.04/1000.
110     tcrit = 1.91d2
111     stop
112  else
113     print*,'Variable gas not recognised!'
114     call abort
115  endif
116
117  rmn     = rc/m_n
118  Ttop    = dble(Tstra_rcm)
119  Tsurf   = dble(Tsurf_rcm)
120
121  psat_v  = psat_max
122  if(vgas.gt.0)then
123     if(trim(gnom(vgas)).eq.'H2O')then
124        call Psat_H2O(tsurf,psat_v)
125     elseif(trim(gnom(vgas)).eq.'NH3')then
126        call Psat_NH3(tsurf,psat_v)
127     endif
128  endif
129
130  ! Moist adiabat unless greater than or equal to psat_max
131  if(psat_v*1d6.lt.psat_max)then
132    Psurf_v = Psat_v*1d6
133    profil_flag(1) = 1
134  else
135    Psurf_v = psat_max
136    profil_flag(1) = 0
137  endif
138
139  if(add_Pvar_to_total)then
140    Psurf_n = dble(psurf_rcm)
141    psurf_rcm = real(Psurf_n+Psurf_v)
142  else
143    Psurf_n = dble(psurf_rcm) - Psurf_v
144  endif
145
146  ! include relative humidity option
147  !if(satval.lt.1.0)then
148  !   Psurf_v = Psurf_v*satval
149  !   profil_flag(1) = 0     
150  !endif
151
152  if(verbose)then
153     print*,'Psat_v  =',psat_v*1d6
154     print*,'Tsurf   =',Tsurf,' K'
155     print*,'Ttop    =',Ttop,' K'
156     print*,'Psurf_v =',Psurf_v,' Pa'
157     print*,'Psurf_n =',Psurf_n,' Pa'
158     print*,'m_n     =',m_n,' kg/mol'
159     print*,'m_v     =',m_v,' kg/mol'
160     print*,'rc      =',rc
161  endif
162
163  ! define fine pressure grid
164  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayer
165
166  P_rcm(1)  = psurf_rcm*exp(dlogp_rcm)
167  do ilay_rcm=1,nlayer-1
168     P_rcm(ilay_rcm+1) = P_rcm(ilay_rcm)*exp(dlogp_rcm)
169  enddo
170
171  Pl_rcm(1) = psurf_rcm
172  do ilev_rcm=2,nlayer
173     ! log-linear interpolation
174     Pl_rcm(ilev_rcm) = exp( log( P_rcm(ilev_rcm)*P_rcm(ilev_rcm-1) )/2 )
175  enddo
176
177  !-------------------------------
178  ! Layer 1
179  T(1)     = Tsurf
180  Pv(1)    = Psurf_v
181  Pn(1)    = Psurf_n
182  rho_n(1) = m_n*Pn(1)/(Rc*Tsurf)
183  rho_v(1) = m_v*Pv(1)/(Rc*Tsurf)
184  a_v(1)   = rho_v(1)/rho_n(1)
185
186  ! log pressure grid spacing (constant)
187  dlogp = -(log(Pn(1)+Pv(1))-log(ptop))/(nlay-1)
188
189  call gradients_kcm(profil_flag(1),rho_v(1),rho_n(1),Tsurf,dTdp,dPvdp,dPndp)
190  if(verbose)then
191     print*, 'dT/dp ground [K/Pa] =',dTdp
192  endif
193
194  ! initial delta p, delta z
195  Dp = (Pn(1) + Pv(1))*(exp(dlogp) - 1d0)
196  Dz = -Dp/(  g*(rho_n(1) + rho_v(1))  )
197
198  !-------------------------------
199  ! Layer 2
200  T(2)     = tsurf + dTdp*Dp
201  Pv(2)    = Pv(1) + dPvdp*Dp
202  Pn(2)    = Pn(1) + dPndp*Dp
203  rho_n(2) = m_n*Pn(2)/(Rc*T(2))
204  rho_v(2) = m_v*Pv(2)/(Rc*T(2))
205  a_v(2)   = rho_v(2)/rho_n(2)
206
207  !-------------------------------
208  ! start vertical ascent
209  Ztab(1) = 0.
210  do ilay=2,nlay-1
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
231     if(vgas.gt.0)then
232     if(trim(gnom(vgas)).eq.'H2O')then
233        call Psat_H2O(T(ilay+1),psat_v)
234     elseif(trim(gnom(vgas)).eq.'NH3')then
235        call Psat_NH3(T(ilay+1),psat_v)
236     endif
237     endif
238
239     if (psat_v*1d6 .lt. Pv(ilay)+dPvdp*Dp) then
240        profil_flag(ilay)=1
241        call gradients_kcm(profil_flag(ilay),rho_v(ilay),rho_n(ilay),&
242                         T(ilay),dTdp,dPvdp,dPndp)
243     endif
244
245     ! test for stratosphere at next level
246     if (T(ilay+1) .le. Ttop) then
247        profil_flag(ilay)=2
248        T(ilay+1)=Ttop
249     endif
250
251     ! calculate pressures at next level
252     Pn(ilay+1) = Pn(ilay) + dPndp*Dp
253     Pv(ilay+1) = Pv(ilay) + dPvdp*Dp
254
255     if(profil_flag(ilay) .eq. 1)then
256       
257        psat_v=psat_max
258
259        if(vgas.gt.0)then
260        if(trim(gnom(vgas)).eq.'H2O')then
261           call Psat_H2O(T(ilay+1),psat_v)
262        elseif(trim(gnom(vgas)).eq.'NH3')then
263           call Psat_NH3(T(ilay+1),psat_v)
264        endif
265        endif
266
267        if(Pv(ilay+1) .lt. psat_v*1e6)then
268           Pv(ilay+1)=psat_v*1d6
269        endif
270
271     endif
272
273     ! calculate gas densities at next level (assume ideal)
274     rho_n(ilay+1) = m_n*Pn(ilay+1)/(rc*T(ilay+1))
275     select case(profil_flag(ilay))
276     case(2) ! isothermal
277        rho_v(ilay+1) = rho_v(ilay)/rho_n(ilay)*rho_n(ilay+1)
278     case(1) ! moist
279
280        ! dont think this is necessary
281        !call psat_est(T(ilay+1),psat_v)
282        ! modify for ammonia!!!
283
284        rho_v(ilay+1) = m_v*psat_v*1d6/(rc*T(ilay+1))
285     case(0) ! dry
286        rho_v(ilay+1) = m_v*Pv(ilay+1)/(rc*T(ilay+1))
287     end select
288
289  enddo
290
291  Ztab(nlay)=Ztab(nlay-1)+Dz
292
293  !-------------------------------
294  ! save to kcm1d variables
295
296  ! surface quantities
297  psurf_rcm = Pn(1) + Pv(1)
298  qsurf_rcm = rho_v(1)/(rho_v(1) + rho_n(1))
299
300  ! create q_v, mtot for saving
301  do ilay=1,nlay
302     mtot(ilay) = 1d3*(rho_v(ilay) + rho_n(ilay)) / &
303                  (rho_v(ilay)/m_v + rho_n(ilay)/m_n)
304     q_v(ilay)  = rho_v(ilay)/(rho_v(ilay) + rho_n(ilay))
305     ! CHECK THIS
306  enddo
307
308
309  ! convert to rcm lower-res grid
310  z_rcm(:) = 0.0
311  T_rcm(:) = 0.0
312  q_rcm(:) = 0.0
313  m_rcm(:) = 0.0
314
315  m_rcm(1) = real( 1d3*(rho_v(1) + rho_n(1)) / &
316                  (rho_v(1)/m_v + rho_n(1)/m_n) )
317
318  ilay_rcm=1
319  do ilay=2,nlay
320
321     if(ilay_rcm.le.nlayer)then
322     ! interpolate rcm variables
323
324        if(Pn(ilay)+Pv(ilay) .lt. P_rcm(ilay_rcm))then
325
326           if(ilay.eq.1)then
327              print*,'Error in create_profils: Psurf here less than Psurf in RCM!'
328              call abort
329           endif
330
331           lnp1   = log(Pn(ilay-1)+Pv(ilay-1))
332           lnp2   = log(Pn(ilay)+Pv(ilay))
333           lnpnew = dble(log(P_rcm(ilay_rcm)))
334
335           z_rcm(ilay_rcm) = real(Ztab(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
336                             + Ztab(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
337           T_rcm(ilay_rcm) = real(T(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
338                             + T(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
339           q_rcm(ilay_rcm) = real(q_v(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
340                             + q_v(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
341
342           m_rcm(ilay_rcm+1) = real(mtot(ilay-1)*(lnp2-lnpnew)/(lnp2-lnp1) &
343                             + mtot(ilay)*(lnpnew-lnp1)/(lnp2-lnp1))
344
345           ilay_rcm = ilay_rcm+1
346        endif
347
348     endif
349  enddo
350
351  ifinal_rcm=ilay_rcm-1
352  if(ifinal_rcm.lt.nlayer)then
353     if(verbose)then
354        print*,'Interpolation in kcmprof stopped at layer',ilay_rcm,'!'
355     endif
356
357     do ilay_rcm=ifinal_rcm+1,nlayer
358
359        z_rcm(ilay_rcm) = z_rcm(ilay_rcm-1)
360        T_rcm(ilay_rcm) = T_rcm(ilay_rcm-1)
361        q_rcm(ilay_rcm) = q_rcm(ilay_rcm-1)
362        m_rcm(ilay_rcm+1) = m_rcm(ilay_rcm)
363
364     enddo
365  endif
366
367  do ilay=2,nlayer
368     if(T_rcm(ilay).lt.Ttop)then
369        T_rcm(ilay)=Ttop
370     endif
371  enddo
372
373!    CO2 condensation 'haircut' of temperature profile if necessary
374  if(co2cond)then
375     print*,'CO2 condensation haircut - assumes CO2-dominated atmosphere!'
376     do ilay=2,nlayer
377        if(P_rcm(ilay).lt.518000.)then
378           TCO2cond = (-3167.8)/(log(.01*P_rcm(ilay))-23.23) ! Fanale's formula
379        else
380           TCO2cond = 684.2-92.3*log(P_rcm(ilay))+4.32*log(P_rcm(ilay))**2
381           ! liquid-vapour transition (based on CRC handbook 2003 data)
382        endif
383
384        print*,'p=',P_rcm(ilay),', T=',T_rcm(ilay),' Tcond=',TCO2cond
385        if(T_rcm(ilay).lt.TCO2cond)then
386           T_rcm(ilay)=TCO2cond
387        endif
388     enddo
389  endif
390
391  return
392end subroutine kcmprof_fn
Note: See TracBrowser for help on using the repository browser.