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

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

Generic GCM:

  • Some code cleanup: turning comcstfi.h into module comcstfi_mod.F90

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