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

Diff of /trunk/phylmd/ozonecm.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/phylmd/ozonecm.f revision 118 by guez, Thu Dec 18 17:30:24 2014 UTC trunk/phylmd/ozonecm.f90 revision 328 by guez, Thu Jun 13 14:40:06 2019 UTC
# Line 16  contains Line 16  contains
16      ! also provided us the code.      ! also provided us the code.
17    
18      ! A. J. Krueger and R. A. Minzner, A Mid-Latitude Ozone Model for the      ! 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).      ! 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      ! Keating 1985 k1032
     ! middle atmosphere, Handbook for MAP, vol. 16, 205-229.  
22    
23      use dimens_m, only: llm      use nr_util, only: assert, twopi, deg_to_rad
24      USE dimphy, ONLY : klon  
25      use nr_util, only: assert, pi      use dimensions, only: llm
26      use phyetat0_m, only: rlat      use phyetat0_m, only: rlat
27    
28      REAL, INTENT (IN) :: rjour      REAL, INTENT (IN) :: rjour
29    
30      REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)      REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm + 1)
31      ! pression pour chaque inter-couche, en Pa      ! pression pour chaque inter-couche, en Pa
32    
33      REAL ozonecm(klon, llm)      REAL ozonecm(size(paprs, 1), llm) ! (klon, llm)
34      ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is      ! "ozonecm(j, k)" is the column-density of ozone in cell "(j, k)", that is
35      ! between interface "k" and interface "k + 1", in kDU.      ! between interface "k" and interface "k + 1", in kDU.
36    
# Line 40  contains Line 39  contains
39      REAL tozon ! equivalent pressure of ozone above interface "k", in Pa      REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
40      INTEGER i, k      INTEGER i, k
41    
42      REAL field(llm+1)      REAL field(llm + 1)
43      ! "field(k)" is the column-density of ozone between interface      ! "field(k)" is the column-density of ozone between interface
44      ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.      ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
45    
# Line 53  contains Line 52  contains
52    
53      !----------------------------------------------------------      !----------------------------------------------------------
54    
55      call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")      call assert(size(paprs, 2) == llm + 1, "ozonecm")
56    
57      sint = sin(2 * pi * (rjour + 15.) / an)      sint = sin(twopi * (rjour + 15.) / an)
58      cost = cos(2 * pi * (rjour + 15.) / an)      cost = cos(twopi * (rjour + 15.) / an)
59      field(llm + 1) = 0.      field(llm + 1) = 0.
60    
61      DO i = 1, klon      DO i = 1, size(paprs, 1)
62         slat = sin(pi / 180. * rlat(i))         slat = sin(deg_to_rad * rlat(i))
63         slat2 = slat * slat         slat2 = slat * slat
64         gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &         gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &
65              * (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &              * (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &
# Line 70  contains Line 69  contains
69         bsec = 2650. + 5000. * slat2         bsec = 2650. + 5000. * slat2
70         a = 4. * bsec**1.5 * ppm**1.5 * (1. + (bsec / ps)**1.5) &         a = 4. * bsec**1.5 * ppm**1.5 * (1. + (bsec / ps)**1.5) &
71              / (bsec**1.5 + ppm**1.5)**2              / (bsec**1.5 + ppm**1.5)**2
72         aprim = max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &         asec = max(0., (gms - max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &
73              - 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))              - 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))) &
74         asec = max(0., (gms - aprim) * (1. + (bsec / ps)**1.5))              * (1. + (bsec / ps)**1.5))
75         aprim = gms - asec / (1. + (bsec / ps)**1.5)         aprim = gms - asec / (1. + (bsec / ps)**1.5)
76    
77         DO k = 1, llm         DO k = 1, llm
Legend:
Removed from v.118  
changed lines
  Added in v.328
  ViewVC Help
Powered by ViewVC 1.1.21