File: rayigc.f90

package info (click to toggle)
code-saturne 4.3.3%2Brepack-1
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 77,992 kB
  • sloc: ansic: 281,257; f90: 122,305; python: 56,490; makefile: 3,915; xml: 3,285; cpp: 3,183; sh: 1,139; lex: 176; yacc: 101; sed: 16
file content (201 lines) | stat: -rw-r--r-- 6,503 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
!-------------------------------------------------------------------------------

! This file is part of Code_Saturne, a general-purpose CFD tool.
!
! Copyright (C) 1998-2016 EDF S.A.
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 2 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.

!-------------------------------------------------------------------------------
!> \file rayigc.f90
!> \brief 1D Radiative scheme - IR CO2 + O3 absorbtion

!> \brief Compute carbonic gaz and ozone absorbtion in infrared
!>
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role
!______________________________________________________________________________!
!> \param[in]       zbas    ground level altitude
!> \param[in]       zz      height above ground level
!> \param[in]       pz      pressure normalized by ground level pressure
!> \param[in]       zzp     intermediate altitude for ozone
!> \param[in]       pzp     corresponding pressure for zzp level
!> \param[out]      xa      CO2 + O3 absorption
!> \param[out]      xda     differential absorption for CO2 + O3
!> \param[in]       q       effective concentration for absorption by water
!>                          vapor
!> \param[in]       u       water vapor optical depth (zz, zzp)
!> \param[in]       tco2    temperature for high level
!> \param[in]       ro      air density
!_______________________________________________________________________________

subroutine rayigc (zbas,zz,pz,zzp,pzp,xa,xda,q,u,tco2,ro)

!===============================================================================
! Module files
!===============================================================================

use paramx
use numvar
use optcal
use cstphy
use entsor
use parall
use period
use ppppar
use ppthch
use ppincl
use atincl

!===============================================================================

implicit none

!... External variables

double precision zbas
double precision zz,pz,zzp,pzp,xa,xda,q,u,tco2,ro

!... Local variables

double precision x1,x2,x3,x4,x1c,x2c,x3c,x4c,y1,y2
double precision tauv,dtauv,xx,exn,exnp1,conco2
double precision uco2,duco2,ao,dao
double precision uo3,duo3,xo1,xo2,xo3,xo4
double precision yo1,yo2,ao3,dao3

data x1,x2,x3,x4/1.33,-0.4572,0.26,0.286/
data x1c,x2c,x3c,x4c/-0.00982,0.0676,0.421,0.01022/
data y1,y2/0.0581,0.0546/
data xo1,xo2,xo3,xo4/0.209,7.e-5,0.436,-0.00321/
data yo1,yo2/0.0749,0.0212/

!===============================================================================

!    Concentration of CO2
conco2 = 3.5d-2

!>   1-Computation of th2o within the range 15mu of Co2
!   ------------------------------------------
if(u.le.20.d0) then
  tauv = x1 + x2*(u + x4)**x3
  dtauv = ro*q*x2*x3*(u + x4)**(x3 - 1.d0)
else
  tauv = 0.33d0 - 0.2754d0*(log10(u) - 1.3011d0)
  dtauv = -0.2754d0/log(10.d0)*ro*q/u
endif

!>   2-Computation of the optical depth for Co2
!   -------------------------------------------
xx = 1.d0 - 0.0065d0*(zz - zbas)/288.15d0
exn = 0.75d0
exnp1 = exn+1.d0
uco2 = -conco2*288.15d0/0.0065d0/(5.31d0*exnp1)*(pz**exnp1                &
     - pzp**exnp1)*(tkelvi/tco2)**(exn/2.d0)

if(uco2.lt.0.d0) uco2 = -uco2
duco2 = conco2*pz**exnp1/xx
duco2 = duco2*(tkelvi/tco2)**(exn/2.d0)
if(uco2.le.1.d0) then
  ao = x1c + x2c*(uco2 + x4c)**x3c
  dao = duco2*x2c*x3c*(uco2 + x4c)**(x3c - 1.d0)
else
  ao = y1 + y2*log10(uco2)
  dao = y2/log(10.d0)*duco2/uco2
endif

!>   3-Computation of the optical depth for O3
!   ---------------------------------------
uo3 = abs(rayuoz(zz) - rayuoz(zzp))
duo3 = raydoz(zz)
if(uo3.le.0.01d0) then
  ao3 = xo1*(uo3 + xo2)**xo3 + xo4
  dao3 = duo3*xo1*(uo3 + xo2)**(xo3 - 1.d0)
else
  ao3 = yo1 + yo2*log10(uo3)
  dao3 = yo2*duo3/log(10.d0)/uo3
endif

!>   4- Compuation of the total absorption (Ozone and co2)
!   -----------------------------------------------
xa = tauv*ao + ao3
xda = tauv*dao + dtauv*ao + dao3

return

!   5- aditionnal functions
!   ------------------------

contains

! 5.1. computes ozones concentration for the altitude zh
!------------------------------------------------------
!> \brief Internal function -
!> computes ozones concentration for the altitude zh
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role
!______________________________________________________________________________!
!> \param[in]       zh      altitude
!_______________________________________________________________________________

function rayuoz(zh)

implicit none
double precision, intent(in) :: zh  ! absolute altitude
double precision ::  rayuoz
double precision ::  a, b, c

a = 0.4d0
b = 20000.d0
c = 5000.d0

rayuoz = a*(1.d0 + exp(-b/c))/(1.d0 + exp((zh - b)/c))

end function rayuoz


! 5.2. computes dO3/dz for the altitude zh
!-------------------------------------------------------
!> \brief Internal function -
!> computes dO3/dz for the altitude zh
!-------------------------------------------------------------------------------
! Arguments
!______________________________________________________________________________.
!  mode           name          role
!______________________________________________________________________________!
!> \param[in]       zh      altitude
!_______________________________________________________________________________

function raydoz(zh)

implicit none
double precision, intent(in) :: zh  ! absolute altitude
double precision ::  raydoz
double precision ::  a, b, c

a = 0.4d0
b = 20000.d0
c = 5000.d0

raydoz = -a/c*(exp((zh - b)/c))*(1.d0 + exp(-b/c))                        &
       / ((1.d0 + exp((zh - b)/c))**2.d0)

end function raydoz

end subroutine rayigc