source: LMDZ4/trunk/libf/phylmd/ozonecm.F @ 757

Last change on this file since 757 was 699, checked in by Laurent Fairhead, 19 years ago

En fait le bug sur l'ozone n'en etait pas un! on revient a la situation precedente MPL/JLD
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.8 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE ozonecm(rjour, rlat, paprs, o3)
5      IMPLICIT none
6C
7C The ozone climatology is based on an analytic formula which fits the
8C Krueger and Mintzner (1976) profile, as well as the variations with
9C altitude and latitude of the maximum ozone concentrations and the total
10C column ozone concentration of Keating and Young (1986). The analytic
11C formula have been established by J-F Royer (CRNM, Meteo France), who
12C also provided us the code.
13C
14C A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the
15C 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).
16C
17C Keating, G. M. and D. F. Young, 1985: Interim reference models for the
18C middle atmosphere, Handbook for MAP, vol. 16, 205-229.
19C
20
21#include "dimensions.h"
22#include "dimphy.h"
23#include "clesphys.h"
24#include "YOMCST.h"
25      REAL rlat(klon), paprs(klon,klev+1)
26      REAL o3(klon,klev)   ! ozone concentration in kg/kg
27      REAL tozon, rjour, pi, pl
28      INTEGER i, k
29C----------------------------------------------------------
30      REAL field(klon,klev+1)
31      REAL ps
32      PARAMETER (ps=101325.0)
33      REAL an, unit, zo3q3
34      SAVE an, unit, zo3q3
35      REAL mu,gms, zslat, zsint, zcost, z, ppm, qpm, a
36      REAL asec, bsec, aprim, zo3a3
37C----------------------------------------------------------
38c         data an /365.25/   (meteo)
39      DATA an /360.00/
40      DATA unit /2.1415e-05/
41      DATA zo3q3 /4.0e-08/
42
43      pi = 4.0 * ATAN(1.0)
44      DO k = 1, klev
45      DO i = 1, klon
46      zslat = SIN(pi/180.*rlat(i))
47      zsint = SIN(2.*pi*(rjour+15.)/an)
48      zcost = COS(2.*pi*(rjour+15.)/an)
49      z = 0.0531+zsint*(-0.001595+0.009443*zslat) +
50     .  zcost*(-0.001344-0.00346*zslat) +
51     .  zslat**2*(.056222+zslat**2*(-.037609
52     . +.012248*zsint+.00521*zcost+.008890*zslat))
53      zo3a3 = zo3q3/ps/2.
54      z = z-zo3q3*ps
55      gms = z
56      mu = ABS(sin(pi/180.*rlat(i)))
57      ppm = 800.-(500.*zslat+150.*zcost)*zslat
58      qpm = 1.74e-5-(7.5e-6*zslat+1.7e-6*zcost)*zslat
59      bsec = 2650.+5000.*zslat**2
60      a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
61      a = a/(bsec**(3./2.)+ppm**(3./2.))**2
62      aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
63      aprim = amax1(0.,aprim)
64      asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
65      asec = AMAX1(0.0,asec)
66      aprim = gms-asec/(1.+(bsec/ps)**(3./2.))
67      pl = paprs(i,k)
68      tozon = aprim/(1.+3.*(ppm/pl)**2)+asec/(1.+(bsec/pl)**(3./2.))
69     .  + zo3a3*pl*pl
70      tozon = tozon / 9.81  ! en kg/m**2
71      tozon = tozon / unit  ! en kg/m**2  > u dobson (10e-5 m)
72      tozon = tozon / 1000. ! en cm
73      field(i,k) = tozon
74      ENDDO
75      ENDDO
76c
77      DO i = 1, klon
78         field(i,klev+1) = 0.0
79      ENDDO
80      DO k = 1, klev
81        DO i = 1, klon
82          o3(i,k) = field(i,k) - field(i,k+1)
83        ENDDO
84      ENDDO
85c
86      RETURN
87      END
Note: See TracBrowser for help on using the repository browser.