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

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

Correction du bug sur l'ozone (unite, calcul). On peut retrouver le
bug, qui est présent dans les simulations IPCC 2005, en positionnant
le flag bug_ozone à true dans physiq.def.
MPL/LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.0 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          IF (.not. bug_ozone) then
84c         convert o3 into kg/kg         
85            o3(i,k)=MAX(o3(i,k),1.0e-12)*RG/46.6968
86     .               /(paprs(i,k)-paprs(i,k+1))
87          ENDIF
88        ENDDO
89      ENDDO
90c
91      RETURN
92      END
Note: See TracBrowser for help on using the repository browser.