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

Diff of /trunk/phylmd/ozonecm.f

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

revision 117 by guez, Tue May 13 17:23:16 2014 UTC revision 118 by guez, Thu Dec 18 17:30:24 2014 UTC
# Line 27  contains Line 27  contains
27      use phyetat0_m, only: rlat      use phyetat0_m, only: rlat
28    
29      REAL, INTENT (IN) :: rjour      REAL, INTENT (IN) :: rjour
30    
31      REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)      REAL, INTENT (IN) :: paprs(:, :) ! (klon, llm+1)
32        ! pression pour chaque inter-couche, en Pa
33    
34      REAL ozonecm(klon, llm)      REAL ozonecm(klon, llm)
35      ! "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
# Line 35  contains Line 37  contains
37    
38      ! Local:      ! Local:
39    
40      REAL tozon, pl      REAL tozon ! equivalent pressure of ozone above interface "k", in Pa
41      INTEGER i, k      INTEGER i, k
42    
43      REAL field(klon, llm+1)      REAL field(llm+1)
44      ! "field(:, k)" is the column-density of ozone between interface      ! "field(k)" is the column-density of ozone between interface
45      ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.      ! "k" and the top of the atmosphere (interface "llm + 1"), in kDU.
46    
47      real, PARAMETER:: ps=101325.      real, PARAMETER:: ps = 101325.
48      REAL, parameter:: an = 360., unit = 2.1415E-5, zo3q3 = 4E-8      REAL, parameter:: an = 360., zo3q3 = 4E-8
49      REAL gms, zslat, zsint, zcost, z, ppm, qpm, a      real, parameter:: zo3a3 = zo3q3 / ps / 2.
50      REAL asec, bsec, aprim, zo3a3      REAL, parameter:: dobson_unit = 2.1415E-5 ! in kg m-2
51        REAL gms, slat, slat2, sint, cost, ppm, a
52        REAL asec, bsec, aprim
53    
54      !----------------------------------------------------------      !----------------------------------------------------------
55    
56      call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")      call assert(shape(paprs) == (/klon, llm + 1/), "ozonecm")
57    
58      DO k = 1, llm      sint = sin(2 * pi * (rjour + 15.) / an)
59         DO i = 1, klon      cost = cos(2 * pi * (rjour + 15.) / an)
60            zslat = sin(pi / 180. * rlat(i))      field(llm + 1) = 0.
61            zsint = sin(2.*pi*(rjour+15.)/an)  
62            zcost = cos(2.*pi*(rjour+15.)/an)      DO i = 1, klon
63            z = 0.0531 + zsint * (-0.001595+0.009443*zslat) &         slat = sin(pi / 180. * rlat(i))
64                 + zcost * (-0.001344-0.00346*zslat) &         slat2 = slat * slat
65                 + zslat**2 * (.056222 + zslat**2 &         gms = 0.0531 + sint * (- 0.001595 + 0.009443 * slat) + cost &
66                 * (-.037609+.012248*zsint+.00521*zcost+.008890*zslat))              * (- 0.001344 - 0.00346 * slat) + slat2 * (0.056222 + slat2 &
67            zo3a3 = zo3q3/ps/2.              * (- 0.037609 + 0.012248 * sint + 0.00521 * cost + 0.00889 &
68            z = z - zo3q3*ps              * slat)) - zo3q3 * ps
69            gms = z         ppm = 800. - 500. * slat2 - 150. * cost * slat
70            ppm = 800. - (500.*zslat+150.*zcost)*zslat         bsec = 2650. + 5000. * slat2
71            qpm = 1.74E-5 - (7.5E-6*zslat+1.7E-6*zcost)*zslat         a = 4. * bsec**1.5 * ppm**1.5 * (1. + (bsec / ps)**1.5) &
72            bsec = 2650. + 5000.*zslat**2              / (bsec**1.5 + ppm**1.5)**2
73            a = 4.0*(bsec)**(3./2.)*(ppm)**(3./2.)*(1.0+(bsec/ps)**(3./2.))         aprim = max(0., (2.666666 * (1.74E-5 - 7.5E-6 * slat2 &
74            a = a/(bsec**(3./2.)+ppm**(3./2.))**2              - 1.7E-6 * cost * slat) * ppm - a * gms) / (1. - a))
75            aprim = (2.666666*qpm*ppm-a*gms)/(1.0-a)         asec = max(0., (gms - aprim) * (1. + (bsec / ps)**1.5))
76            aprim = amax1(0., aprim)         aprim = gms - asec / (1. + (bsec / ps)**1.5)
77            asec = (gms-aprim)*(1.0+(bsec/ps)**(3./2.))  
78            asec = amax1(0.0, asec)         DO k = 1, llm
79            aprim = gms - asec/(1.+(bsec/ps)**(3./2.))            tozon = aprim / (1. + 3. * (ppm / paprs(i, k))**2) &
80            pl = paprs(i, k)                 + asec / (1. + (bsec / paprs(i, k))**1.5) &
81            tozon = aprim / (1. + 3. * (ppm / pl)**2) &                 + zo3a3 * paprs(i, k) * paprs(i, k)
                + asec / (1. + (bsec / pl)**(3./2.)) + zo3a3 * pl * pl  
82            ! Convert from Pa to kDU:            ! Convert from Pa to kDU:
83            field(i, k) = tozon / 9.81 / unit / 1e3            field(k) = tozon / 9.81 / dobson_unit / 1e3
84         END DO         END DO
85    
86           forall (k = 1: llm) ozonecm(i, k) = field(k) - field(k + 1)
87      END DO      END DO
88    
89      field(:, llm+1) = 0.      ozonecm = max(ozonecm, 1e-12)
     forall (k = 1: llm) ozonecm(:, k) = field(:, k) - field(:, k+1)  
90    
91    END function ozonecm    END function ozonecm
92    
Legend:
Removed from v.117  
changed lines
  Added in v.118
  ViewVC Help
Powered by ViewVC 1.1.21