/[lmdze]/trunk/libf/phylmd/Mobidic/regr_pr_comb_coefoz.f90
ViewVC logotype

Diff of /trunk/libf/phylmd/Mobidic/regr_pr_comb_coefoz.f90

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

revision 6 by guez, Wed Feb 27 13:16:39 2008 UTC revision 7 by guez, Mon Mar 31 12:24:17 2008 UTC
# Line 1  Line 1 
1  module read_coefoz_m  module regr_pr_comb_coefoz_m
2    
3    ! This module is clean: no C preprocessor directive, no include line.    ! This module is clean: no C preprocessor directive, no include line.
4    
# Line 7  module read_coefoz_m Line 7  module read_coefoz_m
7    
8    implicit none    implicit none
9    
10    real, save:: c_Mob(klon, llm, 12)    real, save:: c_Mob(klon, llm)
11    ! (sum of Mobidic terms in the net mass production rate of ozone    ! (sum of Mobidic terms in the net mass production rate of ozone
12    ! by chemistry, per unit mass of air, in s-1)    ! by chemistry, per unit mass of air, in s-1)
13    ! (On the "physics" grid.    ! (On the "physics" grid.
14    ! Third dimension is the number of the month in the year.    ! "c_Mob(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
   ! "c_Mob(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
15    ! middle of layer "k".)    ! middle of layer "k".)
16    
17    real, save:: a2(klon, llm, 12)    real, save:: a2(klon, llm)
18    ! (derivative of mass production rate of ozone per unit mass of    ! (derivative of mass production rate of ozone per unit mass of
19    ! air with respect to ozone mass fraction, in s-1)    ! air with respect to ozone mass fraction, in s-1)
20    ! (On the "physics" grid.    ! (On the "physics" grid.
21    ! Third dimension is the number of the month in the year.    ! "a2(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
   ! "a2(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
22    ! middle of layer "k".)    ! middle of layer "k".)
23    
24    real, save:: a4_mass(klon, llm, 12)    real, save:: a4_mass(klon, llm)
25    ! (derivative of mass production rate of ozone per unit mass of    ! (derivative of mass production rate of ozone per unit mass of
26    ! air with respect to temperature, in s-1 K-1)    ! air with respect to temperature, in s-1 K-1)
27    ! (On the "physics" grid.    ! (On the "physics" grid.
28    ! Third dimension is the number of the month in the year.    ! "a4_mass(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
   ! "a4_mass(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
29    ! middle of layer "k".)    ! middle of layer "k".)
30    
31    real, save:: a6_mass(klon, llm, 12)    real, save:: a6_mass(klon, llm)
32    ! (derivative of mass production rate of ozone per unit mass of    ! (derivative of mass production rate of ozone per unit mass of
33    ! air with respect to mass column-density of ozone above, in m2 s-1 kg-1)    ! air with respect to mass column-density of ozone above, in m2 s-1 kg-1)
34    ! (On the "physics" grid.    ! (On the "physics" grid.
35    ! Third dimension is the number of the month in the year.    ! "a6_mass(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
   ! "a6_mass(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
36    ! middle of layer "k".)    ! middle of layer "k".)
37    
38    real, save:: r_het_interm(klon, llm, 12)    real, save:: r_het_interm(klon, llm)
39    ! (net mass production rate by heterogeneous chemistry, per unit    ! (net mass production rate by heterogeneous chemistry, per unit
40    ! mass of ozone, corrected for chlorine content and latitude, but    ! mass of ozone, corrected for chlorine content and latitude, but
41    ! not for temperature and sun direction, in s-1)    ! not for temperature and sun direction, in s-1)
42    ! (On the "physics" grid.    ! (On the "physics" grid.
43    ! Third dimension is the number of the month in the year.    ! "r_het_interm(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
   ! "r_het_interm(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
44    ! middle of layer "k".)    ! middle of layer "k".)
45    
46    private klon, llm    private klon, llm
47    
48  contains  contains
49    
50    subroutine read_coefoz    subroutine regr_pr_comb_coefoz(julien)
51    
52      ! This subroutine reads from a file all eight parameters for ozone      ! "regr_pr_comb_coefoz" stands for "regrid pressure combine
53      ! chemistry, at all months.      ! coefficients ozone".
     ! The parameters are packed to the "physics" grid.  
     ! Finally, the eight parameters are combined to define the five  
     ! module variables.  
54    
55      use netcdf95, only: nf95_open, nf95_close, nf90_nowrite      ! This subroutine :
56      use o3_Mob_ph_m, only: o3_Mob_ph      ! -- reads from a file all eight parameters for ozone chemistry,
57        !    at the current day ;
58        ! -- regrids the parameters in pressure to the LMDZ vertical grid ;
59        ! -- packs the parameters to the "physics" horizontal grid ;
60        ! -- combines the eight parameters to define the five module variables.
61    
62        ! We assume that, in "coefoz_LMDZ.nc", the pressure levels are in hPa
63        ! and strictly increasing.
64    
65        use netcdf95, only: nf95_open, nf95_close, nf95_get_coord
66        use netcdf, only: nf90_nowrite
67        use regr_pr_coefoz, only: regr_pr_av_coefoz, regr_pr_int_coefoz
68      use phyetat0_m, only: rlat      use phyetat0_m, only: rlat
69    
70        integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
71    
72      ! Variables local to the procedure:      ! Variables local to the procedure:
73      integer ncid ! for NetCDF      integer ncid ! for NetCDF
74    
75      real a6(klon, llm, 12)      real, pointer:: plev(:)
76      ! (derivative of P_net_Mob with respect to column-density of ozone      ! (pressure level of input data, converted to Pa, in strictly
77        ! increasing order)
78    
79        integer n_plev ! number of pressure levels in the input data
80    
81        real, allocatable:: press_in_edg(:)
82        ! (edges of pressure intervals for input data, in Pa, in strictly
83        ! increasing order)
84    
85        real coefoz(klon, llm)
86        ! (temporary storage for an ozone coefficient)
87        ! (On the "physics" grid.
88        ! "coefoz(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
89        ! middle of layer "k".)
90    
91        real a6(klon, llm)
92        ! (derivative of "P_net_Mob" with respect to column-density of ozone
93      ! above, in cm2 s-1)      ! above, in cm2 s-1)
94      ! (On the "physics" grid.      ! (On the "physics" grid.
95      ! Third dimension is the number of the month in the year.      ! "a6(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
     ! "a6(i, k, month)" is at longitude "rlon(i)", latitude "rlat(i)",  
96      ! middle of layer "k".)      ! middle of layer "k".)
97    
98      real, parameter:: amu = 1.6605402e-27 ! atomic mass unit, in kg      real, parameter:: amu = 1.6605402e-27 ! atomic mass unit, in kg
# Line 80  contains Line 100  contains
100      real, parameter:: Clx = 3.8e-9      real, parameter:: Clx = 3.8e-9
101      ! (total chlorine content in the upper stratosphere)      ! (total chlorine content in the upper stratosphere)
102    
103      integer k, month      integer k
104    
105      !------------------------------------      !------------------------------------
106    
# Line 88  contains Line 108  contains
108    
109      call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)      call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)
110    
111      a2 = o3_Mob_ph(ncid, "a2")      call nf95_get_coord(ncid, "plev", plev)
112      a4_mass = o3_Mob_ph(ncid, "a4") * 48. / 29.      ! Convert from hPa to Pa because "regr_pr_av" and "regr_pr_int" require so:
113      a6 = o3_Mob_ph(ncid, "a6")      plev = plev * 100.
114        n_plev = size(plev)
115    
116        ! Compute edges of pressure intervals:
117        allocate(press_in_edg(n_plev + 1))
118        press_in_edg(1) = 0.
119        ! We choose edges halfway in logarithm:
120        forall (k = 2:n_plev) press_in_edg(k) = sqrt(plev(k - 1) * plev(k))
121        press_in_edg(n_plev + 1) = huge(0.)
122        ! (infinity, but any value guaranteed to be greater than the
123        ! surface pressure would do)
124    
125        call regr_pr_av_coefoz(ncid, "a2", julien, press_in_edg, a2)
126    
127        call regr_pr_av_coefoz(ncid, "a4", julien, press_in_edg, a4_mass)
128        a4_mass = a4_mass * 48. / 29.
129    
130        call regr_pr_av_coefoz(ncid, "a6", julien, press_in_edg, a6)
131    
132      ! Compute "a6_mass" avoiding underflow, do not divide by 1e4      ! Compute "a6_mass" avoiding underflow, do not divide by 1e4
133      ! before dividing by molecular mass:      ! before dividing by molecular mass:
134      a6_mass = a6 / (1e4 * 29. * amu)      a6_mass = a6 / (1e4 * 29. * amu)
135      ! (factor 1e4: conversion from cm2 to m2)      ! (factor 1e4: conversion from cm2 to m2)
136    
137      c_Mob = 48. / 29. * (o3_Mob_ph(ncid, "P_net_Mob") &      ! Combine coefficients to get "c_Mob":
138           - a2 * o3_Mob_ph(ncid, "r_Mob") - a6 * o3_Mob_ph(ncid, "Sigma_Mob")) &      ! (We use as few local variables as possible, in order to spare
139           - a4_mass * o3_Mob_ph(ncid, "temp_Mob")      ! main memory.)
140    
141        call regr_pr_av_coefoz(ncid, "P_net_Mob", julien, press_in_edg, c_Mob)
142    
143        call regr_pr_av_coefoz(ncid, "r_Mob", julien, press_in_edg, coefoz)
144        c_mob = c_mob - a2 * coeofoz
145    
146        call regr_pr_int_coefoz(ncid, "Sigma_Mob", julien, plev, top_value=0., &
147             coefoz)
148        c_mob = (c_mob - a6 * coefoz) * 48. / 29.
149    
150        call regr_pr_av_coefoz(ncid, "temp_Mob", julien, press_in_edg, coefoz)
151        c_mob = c_mob - a4_mass * coefoz
152    
153      r_het_interm = o3_Mob_ph(ncid, "R_Het") * (Clx / 3.8e-9)**2      call regr_pr_av_coefoz(ncid, "R_Het", julien, press_in_edg, r_het_interm)
154      ! Heterogeneous chemistry is only at high latitudes:      ! Heterogeneous chemistry is only at high latitudes:
155      forall (k = 1: llm, month = 1: 12)      forall (k = 1: llm)
156         where (abs(rlat) <= 45.) r_het_interm(:, k, month) = 0.         where (abs(rlat) <= 45.) r_het_interm(:, k) = 0.
157      end forall      end forall
158        where (r_het_interm  /= 0.) r_het_interm = r_het_interm * (Clx / 3.8e-9)**2
159    
160        deallocate(plev) ! pointer
161      call nf95_close(ncid)      call nf95_close(ncid)
162    
163    end subroutine read_coefoz    end subroutine regr_pr_comb_coefoz
164    
165  end module read_coefoz_m  end module regr_pr_comb_coefoz_m

Legend:
Removed from v.6  
changed lines
  Added in v.7

  ViewVC Help
Powered by ViewVC 1.1.21