/[lmdze]/trunk/phylmd/ozonecm.f90
ViewVC logotype

Annotation of /trunk/phylmd/ozonecm.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 36 - (hide annotations)
Thu Dec 2 17:11:04 2010 UTC (13 years, 5 months ago) by guez
Original Path: trunk/libf/phylmd/ozonecm.f90
File size: 3155 byte(s)
Now using the library "NR_util".

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     ! 1976 U.S. Standard Atmosphere, J. Geophys. Res., 81, 4477, (1976).
20    
21     ! Keating, G. M. and D. F. Young, 1985: Interim reference models for the
22     ! middle atmosphere, Handbook for MAP, vol. 16, 205-229.
23    
24     use dimens_m, only: llm
25     USE dimphy, ONLY : klon
26 guez 36 use nr_util, only: assert, pi
27 guez 22 use phyetat0_m, only: rlat
28    
29     REAL, INTENT (IN) :: rjour
30     REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)
31    
32     REAL ozonecm(klon, llm)
33     ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
34     ! between interface "k" and interface "k + 1", in kDU.
35    
36     ! Variables local to the procedure:
37    
38     REAL tozon, pl
39     INTEGER i, k
40    
41     REAL field(klon, llm+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    
45     real, PARAMETER:: ps=101325.
46     REAL, parameter:: an = 360., unit = 2.1415E-5, zo3q3 = 4E-8
47     REAL gms, zslat, zsint, zcost, z, ppm, qpm, a
48     REAL asec, bsec, aprim, zo3a3
49    
50     !----------------------------------------------------------
51    
52     call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")
53    
54     DO k = 1, llm
55     DO i = 1, klon
56     zslat = sin(pi / 180. * rlat(i))
57     zsint = sin(2.*pi*(rjour+15.)/an)
58     zcost = cos(2.*pi*(rjour+15.)/an)
59     z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &
60     + zcost * (-0.001344-0.00346*zslat) &
61     + zslat**2 * (.056222 + zslat**2 &
62     * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))
63     zo3a3 = zo3q3/ps/2.
64     z = z - zo3q3*ps
65     gms = z
66     ppm = 800. - (500.*zslat+150.*zcost)*zslat
67     qpm = 1.74E-5 - (7.5E-6*zslat+1.7E-6*zcost)*zslat
68     bsec = 2650. + 5000.*zslat**2
69     a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))
70     a = a/(bsec**(3./2.)+ppm**(3./2.))**2
71     aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)
72     aprim = amax1(0., aprim)
73     asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))
74     asec = amax1(0.0, asec)
75     aprim = gms - asec/(1.+(bsec/ps)**(3./2.))
76     pl = paprs(i, k)
77     tozon = aprim / (1. + 3. * (ppm / pl)**2) &
78     + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl
79     ! Convert from Pa to kDU:
80     field(i, k) = tozon / 9.81 / unit / 1e3
81     END DO
82     END DO
83    
84     field(:, llm+1) = 0.
85     forall (k = 1: llm) ozonecm(:, k) = field(:, k) - field(:, k+1)
86    
87     END function ozonecm
88    
89     end module ozonecm_m

  ViewVC Help
Powered by ViewVC 1.1.21