source: LMDZ6/branches/Amaury_dev/libf/phylmd/ozonecm_m.F90 @ 5409

Last change on this file since 5409 was 5119, checked in by abarral, 12 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.5 KB
RevLine 
[524]1! $Header$
[1188]2module ozonecm_m
[641]3
[1188]4  IMPLICIT NONE
[524]5
[5119]6CONTAINS
[1188]7
[1955]8  function ozonecm(rlat, paprs,read_climoz, rjour)
[1188]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
[5116]23    USE dimphy, ONLY: klon, klev
[5117]24    USE lmdz_assert, ONLY: assert
[1188]25
[1215]26    REAL, INTENT (IN) :: rlat(:) ! (klon)
27    REAL, INTENT (IN) :: paprs(:, :) ! (klon,klev+1)
[1220]28    REAL, INTENT (IN) :: rjour
[1955]29    INTEGER, INTENT (IN) :: read_climoz
[1215]30
31    REAL ozonecm(klon,klev)
32    ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
[1188]33    ! between interface "k" and interface "k + 1", in kDU.
34
35    ! Variables local to the procedure:
36
37    REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
[5117]38    REAL pi, pl
[1188]39    INTEGER i, k
40
41    REAL field(klon,klev+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
[5117]45    REAL, PARAMETER:: ps=101325.
[1188]46    REAL, parameter:: an = 360., zo3q3 = 4E-8
47    REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
[1955]48    REAL gms, zslat,zslat2, zsint, zcost, z, ppm, qpm, a
[1188]49    REAL asec, bsec, aprim, zo3a3
50
51    !----------------------------------------------------------
52
[5103]53    CALL assert((/size(rlat), size(paprs, 1)/) == klon, "ozonecm klon")
54    CALL assert(size(paprs, 2) == klev + 1, "ozonecm klev")
[1215]55
[1188]56    pi = 4. * atan(1.)
57    DO k = 1, klev
58       DO i = 1, klon
59          zslat = sin(pi / 180. * rlat(i))
[1955]60          zslat2=zslat*zslat
61          IF (read_climoz==-1) zslat=0. ! Imposing hemispheric symetry
[1215]62          zsint = sin(2 * pi * (rjour + 15.) / an)
63          zcost = cos(2 * pi * (rjour + 15.) / an)
[1188]64          z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
65               + zcost * (-0.001344-0.00346*zslat) &
[1955]66               + zslat2 * (.056222 + zslat2 &
[1188]67               * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))
68          zo3a3 = zo3q3/ps/2.
69          z = z - zo3q3*ps
70          gms = z
[1955]71          ppm = 800. - 500.*zslat2 - 150.*zcost*zslat
72          qpm = 1.74E-5 - 7.5E-6*zslat2 - 1.7E-6*zcost*zslat
73          bsec = 2650. + 5000.*zslat2
[1188]74          a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
75          a = a/(bsec**(3./2.)+ppm**(3./2.))**2
76          aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
77          aprim = amax1(0., aprim)
78          asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
79          asec = amax1(0.0, asec)
80          aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
81          pl = paprs(i, k)
82          tozon = aprim / (1. + 3. * (ppm / pl)**2) &
83               + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl
84          ! Convert from Pa to kDU:
85          field(i, k) = tozon / 9.81 / dobson_unit / 1e3
86       END DO
87    END DO
88
89    field(:,klev+1) = 0.
[1215]90    forall (k = 1: klev) ozonecm(:,k) = field(:,k) - field(:,k+1)
[1231]91    ozonecm = max(ozonecm, 1e-12)
[1188]92
[5103]93!    PRINT*,'ozonecm Version2'
[1955]94
[1215]95  END function ozonecm
[1188]96
[5119]97END MODULE ozonecm_m
Note: See TracBrowser for help on using the repository browser.