trunk/phylmd/ozonecm.f90

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 ozone
    ! models for the middle atmosphere, Handbook for MAP, vol. 16,
    ! 205-229.

    use nr_util, only: assert, twopi, deg_to_rad

    use dimensions, only: llm
    use phyetat0_m, only: rlat

    REAL, INTENT (IN) :: rjour

    REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm + 1)
    ! pression pour chaque inter-couche, en Pa

    REAL ozonecm(size(paprs, 1), llm) ! (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 ! equivalent pressure of ozone above interface "k", in Pa
    INTEGER i, k

    REAL field(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., zo3q3 = 4E-8
    real, parameter:: zo3a3 = zo3q3 / ps / 2.
    REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
    REAL gms, slat, slat2, sint, cost, ppm, a
    REAL asec, bsec, aprim

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

    call assert(size(paprs, 2) == llm + 1, "ozonecm")

    sint = sin(twopi * (rjour + 15.) / an)
    cost = cos(twopi * (rjour + 15.) / an)
    field(llm + 1) = 0.

    DO i = 1, size(paprs, 1)
       slat = sin(deg_to_rad * rlat(i))
       slat2 = slat * slat
       gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &
            * (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &
            * (- 0.037609 + 0.012248 * sint + 0.00521 * cost + 0.00889 &
            * slat)) - zo3q3 * ps
       ppm = 800. - 500. * slat2 - 150. * cost * slat
       bsec = 2650. + 5000. * slat2
       a = 4. * bsec**1.5 * ppm**1.5 * (1. + (bsec / ps)**1.5) &
            / (bsec**1.5 + ppm**1.5)**2
       asec = max(0., (gms - max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &
            - 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))) &
            * (1. + (bsec / ps)**1.5))
       aprim = gms - asec / (1. + (bsec / ps)**1.5)

       DO k = 1, llm
          tozon = aprim / (1. + 3. * (ppm / paprs(i, k))**2) &
               + asec / (1. + (bsec / paprs(i, k))**1.5) &
               + zo3a3 * paprs(i, k) * paprs(i, k)
          ! Convert from Pa to kDU:
          field(k) = tozon / 9.81 / dobson_unit / 1e3
       END DO

       forall (k = 1: llm) ozonecm(i, k) = field(k) - field(k + 1)
    END DO

    ozonecm = max(ozonecm, 1e-12)

  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	guez	315	! 1976 U. S. Standard Atmosphere, J. Geophys. Res., 81, 4477 (1976).
20	guez	22
21	guez	315	! Keating, G. M. and D. F. Young, 1985: Interim reference ozone
22			! models for the middle atmosphere, Handbook for MAP, vol. 16,
23			! 205-229.
24	guez	22
25	guez	315	use nr_util, only: assert, twopi, deg_to_rad
26	guez	278
27	guez	265	use dimensions, only: llm
28	guez	22	use phyetat0_m, only: rlat
29
30			REAL, INTENT (IN) :: rjour
31	guez	118
32	guez	266	REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm + 1)
33	guez	118	! pression pour chaque inter-couche, en Pa
34	guez	22
35	guez	266	REAL ozonecm(size(paprs, 1), llm) ! (klon, llm)
36	guez	22	! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
37			! between interface "k" and interface "k + 1", in kDU.
38
39	guez	98	! Local:
40	guez	22
41	guez	118	REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
42	guez	22	INTEGER i, k
43
44	guez	266	REAL field(llm + 1)
45	guez	118	! "field(k)" is the column-density of ozone between interface
46	guez	22	! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
47
48	guez	118	real, PARAMETER:: ps = 101325.
49			REAL, parameter:: an = 360., zo3q3 = 4E-8
50			real, parameter:: zo3a3 = zo3q3 / ps / 2.
51			REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
52			REAL gms, slat, slat2, sint, cost, ppm, a
53			REAL asec, bsec, aprim
54	guez	22
55			!----------------------------------------------------------
56
57	guez	266	call assert(size(paprs, 2) == llm + 1, "ozonecm")
58	guez	22
59	guez	315	sint = sin(twopi * (rjour + 15.) / an)
60			cost = cos(twopi * (rjour + 15.) / an)
61	guez	118	field(llm + 1) = 0.
62
63	guez	266	DO i = 1, size(paprs, 1)
64	guez	315	slat = sin(deg_to_rad * rlat(i))
65	guez	118	slat2 = slat * slat
66			gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &
67			* (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &
68			* (- 0.037609 + 0.012248 * sint + 0.00521 * cost + 0.00889 &
69			* slat)) - zo3q3 * ps
70			ppm = 800. - 500. * slat2 - 150. * cost * slat
71			bsec = 2650. + 5000. * slat2
72			a = 4. * bsec*1.5 ppm*1.5 (1. + (bsec / ps)**1.5) &
73			/ (bsec1.5 + ppm1.5)**2
74	guez	316	asec = max(0., (gms - max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &
75			- 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))) &
76			* (1. + (bsec / ps)**1.5))
77	guez	118	aprim = gms - asec / (1. + (bsec / ps)**1.5)
78
79			DO k = 1, llm
80			tozon = aprim / (1. + 3. * (ppm / paprs(i, k))**2) &
81			+ asec / (1. + (bsec / paprs(i, k))**1.5) &
82			+ zo3a3 * paprs(i, k) * paprs(i, k)
83	guez	22	! Convert from Pa to kDU:
84	guez	118	field(k) = tozon / 9.81 / dobson_unit / 1e3
85	guez	22	END DO
86	guez	118
87			forall (k = 1: llm) ozonecm(i, k) = field(k) - field(k + 1)
88	guez	22	END DO
89
90	guez	118	ozonecm = max(ozonecm, 1e-12)
91	guez	22
92			END function ozonecm
93
94			end module ozonecm_m