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

Last change on this file since 1397 was 1397, checked in by milmd, 10 years ago

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