source: trunk/libf/phylmd/ozonecm_m.F90 @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 3.3 KB
Line 
1! $Header$
2module ozonecm_m
3
4  IMPLICIT NONE
5
6contains
7
8  function ozonecm(rlat, paprs, 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
30    REAL ozonecm(klon,klev)
31    ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
32    ! between interface "k" and interface "k + 1", in kDU.
33
34    ! Variables local to the procedure:
35
36    REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
37    real pi, pl
38    INTEGER i, k
39
40    REAL field(klon,klev+1)
41    ! "field(:, k)" is the column-density of ozone between interface
42    ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
43
44    real, PARAMETER:: ps=101325.
45    REAL, parameter:: an = 360., zo3q3 = 4E-8
46    REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
47    REAL gms, zslat, zsint, zcost, z, ppm, qpm, a
48    REAL asec, bsec, aprim, zo3a3
49
50    !----------------------------------------------------------
51
52    call assert((/size(rlat), size(paprs, 1)/) == klon, "ozonecm klon")
53    call assert(size(paprs, 2) == klev + 1, "ozonecm klev")
54
55    pi = 4. * atan(1.)
56    DO k = 1, klev
57       DO i = 1, klon
58          zslat = sin(pi / 180. * rlat(i))
59          zsint = sin(2 * pi * (rjour + 15.) / an)
60          zcost = cos(2 * pi * (rjour + 15.) / an)
61          z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
62               + zcost * (-0.001344-0.00346*zslat) &
63               + zslat**2 * (.056222 + zslat**2 &
64               * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))
65          zo3a3 = zo3q3/ps/2.
66          z = z - zo3q3*ps
67          gms = z
68          ppm = 800. - (500.*zslat+150.*zcost)*zslat
69          qpm = 1.74E-5 - (7.5E-6*zslat+1.7E-6*zcost)*zslat
70          bsec = 2650. + 5000.*zslat**2
71          a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
72          a = a/(bsec**(3./2.)+ppm**(3./2.))**2
73          aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
74          aprim = amax1(0., aprim)
75          asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
76          asec = amax1(0.0, asec)
77          aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
78          pl = paprs(i, k)
79          tozon = aprim / (1. + 3. * (ppm / pl)**2) &
80               + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl
81          ! Convert from Pa to kDU:
82          field(i, k) = tozon / 9.81 / dobson_unit / 1e3
83       END DO
84    END DO
85
86    field(:,klev+1) = 0.
87    forall (k = 1: klev) ozonecm(:,k) = field(:,k) - field(:,k+1)
88    ozonecm = max(ozonecm, 1e-12)
89
90  END function ozonecm
91
92end module ozonecm_m
Note: See TracBrowser for help on using the repository browser.