1 |
|
|
! $Header$ |
2 |
|
|
module ozonecm_m |
3 |
|
|
|
4 |
|
|
IMPLICIT NONE |
5 |
|
|
|
6 |
|
|
contains |
7 |
|
|
|
8 |
✗✓✗✓
|
3 |
function ozonecm(rlat, paprs,read_climoz, rjour) |
9 |
|
|
|
10 |
|
|
! The ozone climatology is based on an analytic formula which fits the |
11 |
|
|
! Krueger and Mintzner (1976) profile, as well as the variations with |
12 |
|
|
! altitude and latitude of the maximum ozone concentrations and the total |
13 |
|
|
! column ozone concentration of Keating and Young (1986). The analytic |
14 |
|
|
! formula have been established by J.-F. Royer (CRNM, Meteo France), who |
15 |
|
|
! also provided us the code. |
16 |
|
|
|
17 |
|
|
! A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the |
18 |
|
|
! 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976). |
19 |
|
|
|
20 |
|
|
! Keating, G. M. and D. F. Young, 1985: Interim reference models for the |
21 |
|
|
! middle atmosphere, Handbook for MAP, vol. 16, 205-229. |
22 |
|
|
|
23 |
|
|
USE dimphy, only: klon, klev |
24 |
|
|
use assert_m, only: assert |
25 |
|
|
|
26 |
|
|
REAL, INTENT (IN) :: rlat(:) ! (klon) |
27 |
|
|
REAL, INTENT (IN) :: paprs(:, :) ! (klon,klev+1) |
28 |
|
|
REAL, INTENT (IN) :: rjour |
29 |
|
|
INTEGER, INTENT (IN) :: read_climoz |
30 |
|
|
|
31 |
|
|
REAL ozonecm(klon,klev) |
32 |
|
|
! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is |
33 |
|
|
! between interface "k" and interface "k + 1", in kDU. |
34 |
|
|
|
35 |
|
|
! Variables local to the procedure: |
36 |
|
|
|
37 |
|
|
REAL tozon ! equivalent pressure of ozone above interface "k", in Pa |
38 |
|
|
real pi, pl |
39 |
|
|
INTEGER i, k |
40 |
|
|
|
41 |
✗✓ |
3 |
REAL field(klon,klev+1) |
42 |
|
|
! "field(:, k)" is the column-density of ozone between interface |
43 |
|
|
! "k" and the top of the atmosphere (interface "llm + 1"), in kDU. |
44 |
|
|
|
45 |
|
|
real, PARAMETER:: ps=101325. |
46 |
|
|
REAL, parameter:: an = 360., zo3q3 = 4E-8 |
47 |
|
|
REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2 |
48 |
|
|
REAL gms, zslat,zslat2, zsint, zcost, z, ppm, qpm, a |
49 |
|
|
REAL asec, bsec, aprim, zo3a3 |
50 |
|
|
|
51 |
|
|
!---------------------------------------------------------- |
52 |
|
|
|
53 |
✓✓ |
9 |
call assert((/size(rlat), size(paprs, 1)/) == klon, "ozonecm klon") |
54 |
|
3 |
call assert(size(paprs, 2) == klev + 1, "ozonecm klev") |
55 |
|
|
|
56 |
|
|
pi = 4. * atan(1.) |
57 |
✓✓ |
120 |
DO k = 1, klev |
58 |
✓✓ |
116418 |
DO i = 1, klon |
59 |
|
116298 |
zslat = sin(pi / 180. * rlat(i)) |
60 |
|
116298 |
zslat2=zslat*zslat |
61 |
✗✓ |
116298 |
IF (read_climoz==-1) zslat=0. ! Imposing hemispheric symetry |
62 |
|
116298 |
zsint = sin(2 * pi * (rjour + 15.) / an) |
63 |
|
116298 |
zcost = cos(2 * pi * (rjour + 15.) / an) |
64 |
|
|
z = 0.0531 + zsint * (-0.001595+0.009443*zslat) & |
65 |
|
|
+ zcost * (-0.001344-0.00346*zslat) & |
66 |
|
|
+ zslat2 * (.056222 + zslat2 & |
67 |
|
116298 |
* (-.037609+.012248*zsint+.00521*zcost+.008890*zslat)) |
68 |
|
|
zo3a3 = zo3q3/ps/2. |
69 |
|
116298 |
z = z - zo3q3*ps |
70 |
|
|
gms = z |
71 |
|
116298 |
ppm = 800. - 500.*zslat2 - 150.*zcost*zslat |
72 |
|
116298 |
qpm = 1.74E-5 - 7.5E-6*zslat2 - 1.7E-6*zcost*zslat |
73 |
|
116298 |
bsec = 2650. + 5000.*zslat2 |
74 |
|
116298 |
a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.)) |
75 |
|
116298 |
a = a/(bsec**(3./2.)+ppm**(3./2.))**2 |
76 |
|
116298 |
aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a) |
77 |
|
116298 |
aprim = amax1(0., aprim) |
78 |
|
116298 |
asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.)) |
79 |
|
116298 |
asec = amax1(0.0, asec) |
80 |
|
116298 |
aprim = gms - asec/(1.+(bsec/ps)**(3./2.)) |
81 |
|
116298 |
pl = paprs(i, k) |
82 |
|
|
tozon = aprim / (1. + 3. * (ppm / pl)**2) & |
83 |
|
116298 |
+ asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl |
84 |
|
|
! Convert from Pa to kDU: |
85 |
|
116415 |
field(i, k) = tozon / 9.81 / dobson_unit / 1e3 |
86 |
|
|
END DO |
87 |
|
|
END DO |
88 |
|
|
|
89 |
✓✓ |
2985 |
field(:,klev+1) = 0. |
90 |
✓✓✓✓
|
116418 |
forall (k = 1: klev) ozonecm(:,k) = field(:,k) - field(:,k+1) |
91 |
✓✓✓✓
|
116418 |
ozonecm = max(ozonecm, 1e-12) |
92 |
|
|
|
93 |
|
|
! print*,'ozonecm Version2' |
94 |
|
|
|
95 |
|
|
END function ozonecm |
96 |
|
|
|
97 |
|
|
end module ozonecm_m |