source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/kcmprof_fn.F90 @ 222

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM
  • 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(trim(gnom(vgas)).eq.'H2O')then
102     m_v   = dble(mH2O/1000.)
103     tcrit = 6.47d2
104  elseif(trim(gnom(vgas)).eq.'NH3')then
105     m_v   = 17.031/1000.
106     tcrit = 4.06d2
107  elseif(trim(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(trim(gnom(vgas)).eq.'H2O')then
123        call Psat_H2O(tsurf,psat_v)
124     elseif(trim(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(trim(gnom(vgas)).eq.'H2O')then
232        call Psat_H2O(T(ilay+1),psat_v)
233     elseif(trim(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(trim(gnom(vgas)).eq.'H2O')then
260           call Psat_H2O(T(ilay+1),psat_v)
261        elseif(trim(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.