source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/ozonecm_m.F90 @ 1200

Last change on this file since 1200 was 1188, checked in by lguez, 15 years ago

Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6
Fortran 90 interface.
Corrected a few comments in "output.def".

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