/[lmdze]/trunk/phylmd/ozonecm.f
ViewVC logotype

Contents of /trunk/phylmd/ozonecm.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 278 - (show annotations)
Thu Jul 12 17:53:18 2018 UTC (5 years, 10 months ago) by guez
File size: 3266 byte(s)
Created procedure phyetat0_new to avoid side effects on variables of
module phyetat0_m. Had then to change test_ozonecm: read rlat from
netCDF file instead of specifying it in test_ozonecm.

1 module ozonecm_m
2
3 IMPLICIT NONE
4
5 contains
6
7 function ozonecm(rjour, paprs)
8
9 ! From phylmd/ozonecm.F, version 1.3 2005/06/06 13:16:33
10
11 ! The ozone climatology is based on an analytic formula which fits the
12 ! Krueger and Mintzner (1976) profile, as well as the variations with
13 ! altitude and latitude of the maximum ozone concentrations and the total
14 ! column ozone concentration of Keating and Young (1986). The analytic
15 ! formula have been established by J.-F. Royer (CRNM, Meteo France), who
16 ! also provided us the code.
17
18 ! A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the
19 ! 1976 U. S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).
20
21 ! Keating, G. M. and D. F. Young, 1985: Interim reference models for the
22 ! middle atmosphere, Handbook for MAP, vol. 16, 205-229.
23
24 use nr_util, only: assert, pi
25
26 use dimensions, only: llm
27 use phyetat0_m, only: rlat
28
29 REAL, INTENT (IN) :: rjour
30
31 REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm + 1)
32 ! pression pour chaque inter-couche, en Pa
33
34 REAL ozonecm(size(paprs, 1), llm) ! (klon, llm)
35 ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
36 ! between interface "k" and interface "k + 1", in kDU.
37
38 ! Local:
39
40 REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
41 INTEGER i, k
42
43 REAL field(llm + 1)
44 ! "field(k)" is the column-density of ozone between interface
45 ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
46
47 real, PARAMETER:: ps = 101325.
48 REAL, parameter:: an = 360., zo3q3 = 4E-8
49 real, parameter:: zo3a3 = zo3q3 / ps / 2.
50 REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
51 REAL gms, slat, slat2, sint, cost, ppm, a
52 REAL asec, bsec, aprim
53
54 !----------------------------------------------------------
55
56 call assert(size(paprs, 2) == llm + 1, "ozonecm")
57
58 sint = sin(2 * pi * (rjour + 15.) / an)
59 cost = cos(2 * pi * (rjour + 15.) / an)
60 field(llm + 1) = 0.
61
62 DO i = 1, size(paprs, 1)
63 slat = sin(pi / 180. * rlat(i))
64 slat2 = slat * slat
65 gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &
66 * (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &
67 * (- 0.037609 + 0.012248 * sint + 0.00521 * cost + 0.00889 &
68 * slat)) - zo3q3 * ps
69 ppm = 800. - 500. * slat2 - 150. * cost * slat
70 bsec = 2650. + 5000. * slat2
71 a = 4. * bsec**1.5 * ppm**1.5 * (1. + (bsec / ps)**1.5) &
72 / (bsec**1.5 + ppm**1.5)**2
73 aprim = max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &
74 - 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))
75 asec = max(0., (gms - aprim) * (1. + (bsec / ps)**1.5))
76 aprim = gms - asec / (1. + (bsec / ps)**1.5)
77
78 DO k = 1, llm
79 tozon = aprim / (1. + 3. * (ppm / paprs(i, k))**2) &
80 + asec / (1. + (bsec / paprs(i, k))**1.5) &
81 + zo3a3 * paprs(i, k) * paprs(i, k)
82 ! Convert from Pa to kDU:
83 field(k) = tozon / 9.81 / dobson_unit / 1e3
84 END DO
85
86 forall (k = 1: llm) ozonecm(i, k) = field(k) - field(k + 1)
87 END DO
88
89 ozonecm = max(ozonecm, 1e-12)
90
91 END function ozonecm
92
93 end module ozonecm_m

  ViewVC Help
Powered by ViewVC 1.1.21