trunk/phylmd/ozonecm.f

module ozonecm_m

  IMPLICIT NONE

contains

  function ozonecm(rjour, paprs)

    ! From phylmd/ozonecm.F, version 1.3 2005/06/06 13:16:33

    ! The ozone climatology is based on an analytic formula which fits the
    ! Krueger and Mintzner (1976) profile, as well as the variations with
    ! altitude and latitude of the maximum ozone concentrations and the total
    ! column ozone concentration of Keating and Young (1986). The analytic
    ! formula have been established by J.-F. Royer (CRNM, Meteo France), who
    ! also provided us the code.

    ! A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the
    ! 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).

    ! Keating, G. M. and D. F. Young, 1985: Interim reference models for the
    ! middle atmosphere, Handbook for MAP, vol. 16, 205-229.

    use dimens_m, only: llm
    USE dimphy, ONLY : klon
    use nr_util, only: assert, pi
    use phyetat0_m, only: rlat

    REAL, INTENT (IN) :: rjour
    REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)

    REAL ozonecm(klon, llm)
    ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
    ! between interface "k" and interface "k + 1", in kDU.

    ! Local:

    REAL tozon, pl
    INTEGER i, k

    REAL field(klon, llm+1)
    ! "field(:, k)" is the column-density of ozone between interface
    ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.

    real, PARAMETER:: ps=101325.
    REAL, parameter:: an = 360., unit = 2.1415E-5, zo3q3 = 4E-8
    REAL gms, zslat, zsint, zcost, z, ppm, qpm, a
    REAL asec, bsec, aprim, zo3a3

    !----------------------------------------------------------

    call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")

    DO k = 1, llm
       DO i = 1, klon
          zslat = sin(pi / 180. * rlat(i))
          zsint = sin(2.*pi*(rjour+15.)/an)
          zcost = cos(2.*pi*(rjour+15.)/an)
          z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
               + zcost * (-0.001344-0.00346*zslat) &
               + zslat**2 * (.056222 + zslat**2 &
               * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))
          zo3a3 = zo3q3/ps/2.
          z = z - zo3q3*ps
          gms = z
          ppm = 800. - (500.*zslat+150.*zcost)*zslat
          qpm = 1.74E-5 - (7.5E-6*zslat+1.7E-6*zcost)*zslat
          bsec = 2650. + 5000.*zslat**2
          a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
          a = a/(bsec**(3./2.)+ppm**(3./2.))**2
          aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
          aprim = amax1(0., aprim)
          asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
          asec = amax1(0.0, asec)
          aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
          pl = paprs(i, k)
          tozon = aprim / (1. + 3. * (ppm / pl)**2) &
               + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl
          ! Convert from Pa to kDU:
          field(i, k) = tozon / 9.81 / unit / 1e3
       END DO
    END DO

    field(:, llm+1) = 0.
    forall (k = 1: llm) ozonecm(:, k) = field(:, k) - field(:, k+1)

  END function ozonecm

end module ozonecm_m
1	guez	22	module ozonecm_m
2	guez	3
3	guez	22	IMPLICIT NONE
4	guez	3
5	guez	22	contains
6	guez	3
7	guez	22	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 dimens_m, only: llm
25			USE dimphy, ONLY : klon
26	guez	36	use nr_util, only: assert, pi
27	guez	22	use phyetat0_m, only: rlat
28
29			REAL, INTENT (IN) :: rjour
30			REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)
31
32			REAL ozonecm(klon, llm)
33			! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
34			! between interface "k" and interface "k + 1", in kDU.
35
36	guez	98	! Local:
37	guez	22
38			REAL tozon, pl
39			INTEGER i, k
40
41			REAL field(klon, llm+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., unit = 2.1415E-5, zo3q3 = 4E-8
47			REAL gms, zslat, zsint, zcost, z, ppm, qpm, a
48			REAL asec, bsec, aprim, zo3a3
49
50			!----------------------------------------------------------
51
52			call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")
53
54			DO k = 1, llm
55			DO i = 1, klon
56			zslat = sin(pi / 180. * rlat(i))
57			zsint = sin(2.pi(rjour+15.)/an)
58			zcost = cos(2.pi(rjour+15.)/an)
59			z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
60			+ zcost * (-0.001344-0.00346*zslat) &
61			+ zslat*2 (.056222 + zslat**2 &
62			* (-.037609+.012248zsint+.00521zcost+.008890*zslat))
63			zo3a3 = zo3q3/ps/2.
64			z = z - zo3q3*ps
65			gms = z
66			ppm = 800. - (500.zslat+150.zcost)*zslat
67			qpm = 1.74E-5 - (7.5E-6zslat+1.7E-6zcost)*zslat
68			bsec = 2650. + 5000.zslat*2
69			a = 4.0(bsec)(3./2.)(ppm)*(3./2.)(1.0+(bsec/ps)**(3./2.))
70			a = a/(bsec(3./2.)+ppm(3./2.))**2
71			aprim = (2.666666qpmppm-a*gms)/(1.0-a)
72			aprim = amax1(0., aprim)
73			asec = (gms-aprim)(1.0+(bsec/ps)*(3./2.))
74			asec = amax1(0.0, asec)
75			aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
76			pl = paprs(i, k)
77			tozon = aprim / (1. + 3. * (ppm / pl)**2) &
78			+ asec / (1. + (bsec / pl)*(3./2.)) + zo3a3 pl * pl
79			! Convert from Pa to kDU:
80			field(i, k) = tozon / 9.81 / unit / 1e3
81			END DO
82			END DO
83
84			field(:, llm+1) = 0.
85			forall (k = 1: llm) ozonecm(:, k) = field(:, k) - field(:, k+1)
86
87			END function ozonecm
88
89			end module ozonecm_m