source: LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90 @ 4799

Last change on this file since 4799 was 2788, checked in by dcugnet, 8 years ago

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

  • records 2-13 are obtained with records 1-12 of "climoz.nc".
  • records 1 and 14 are obtained respectively with:
    • record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
    • record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

  • LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
  • file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

  • File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 5.5 KB
Line 
1! $Id$
2module regr_pr_comb_coefoz_m
3
4  implicit none
5
6  ! The five module variables declared here are on the partial
7  ! "physics" grid.
8  ! The value of each variable for index "(i, k)" is at longitude
9  ! "rlon(i)", latitude "rlat(i)" and middle of layer "k".
10
11  real, allocatable, save:: c_Mob(:, :)
12  ! (sum of Mobidic terms in the net mass production rate of ozone
13  ! by chemistry, per unit mass of air, in s-1)
14
15  real, allocatable, save:: a2(:, :)
16  ! (derivative of mass production rate of ozone per unit mass of
17  ! air with respect to ozone mass fraction, in s-1)
18
19  real, allocatable, save:: a4_mass(:, :)
20  ! (derivative of mass production rate of ozone per unit mass of
21  ! air with respect to temperature, in s-1 K-1)
22
23  real, allocatable, save:: a6_mass(:, :)
24  ! (derivative of mass production rate of ozone per unit mass of
25  ! air with respect to mass column-density of ozone above, in m2 s-1 kg-1)
26
27  real, allocatable, save:: r_het_interm(:, :)
28  ! (net mass production rate by heterogeneous chemistry, per unit
29  ! mass of ozone, corrected for chlorine content and latitude, but
30  ! not for temperature and sun direction, in s-1)
31
32  !$omp threadprivate(c_Mob, a2, a4_mass, a6_mass, r_het_interm)
33
34contains
35
36  subroutine alloc_coefoz
37
38    ! This procedure is called once per run.
39    ! It allocates module variables.
40
41    use dimphy, only: klon
42    use mod_grid_phy_lmdz, only: nbp_lev
43
44    ! Variables local to the procedure:
45
46    !---------------------------------------
47
48    !$omp master
49    print *, "Call sequence information: alloc_coefoz"
50    !$omp end master
51    allocate(c_Mob(klon, nbp_lev), a2(klon, nbp_lev), a4_mass(klon, nbp_lev))
52    allocate(a6_mass(klon, nbp_lev), r_het_interm(klon, nbp_lev))
53
54  end subroutine alloc_coefoz
55
56  !*******************************************************
57
58  subroutine regr_pr_comb_coefoz(julien, rlat, paprs, pplay)
59
60    ! "regr_pr_comb_coefoz" stands for "regrid pressure combine
61    ! coefficients ozone".
62
63    ! In this subroutine:
64    ! -- the master thread of the root process reads from a file all
65    !    eight coefficients for ozone chemistry, at the current day;
66    ! -- the coefficients are packed to the "physics" horizontal grid
67    !    and scattered to all threads of all processes;
68    ! -- in all the threads of all the processes, the coefficients are
69    !    regridded in pressure to the LMDZ vertical grid;
70    ! -- in all the threads of all the processes, the eight
71    !    coefficients are combined to define the five module variables.
72
73    use netcdf95, only: nf95_open, nf95_close
74    use netcdf, only: nf90_nowrite
75    use assert_m, only: assert
76    use dimphy, only: klon
77    use mod_phys_lmdz_mpi_data, only: is_mpi_root
78    use regr_pr_time_av_m, only: regr_pr_time_av
79    use regr_pr_int_m, only: regr_pr_int
80    use press_coefoz_m, only: press_in_edg, plev
81    use mod_grid_phy_lmdz, only: nbp_lev
82
83    integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
84
85    REAL, intent(in):: rlat(:)
86    ! (latitude on the partial "physics" grid, in degrees)
87
88    real, intent(in):: paprs(:, :) ! (klon, nbp_lev + 1)
89    ! (pression pour chaque inter-couche, en Pa)
90
91    real, intent(in):: pplay(:, :) ! (klon, nbp_lev)
92    ! (pression pour le mileu de chaque couche, en Pa)
93
94    ! Variables local to the procedure:
95
96    integer ncid ! for NetCDF
97
98    real coefoz(klon, nbp_lev, 7)
99    ! (temporary storage for 7 ozone coefficients)
100    ! (On the partial "physics" grid.
101    ! "coefoz(i, k, :)" is at longitude "rlon(i)", latitude "rlat(i)",
102    ! middle of layer "k".)
103
104    real a6(klon, nbp_lev)
105    ! (derivative of "P_net_Mob" with respect to column-density of ozone
106    ! above, in cm2 s-1)
107    ! (On the partial "physics" grid.
108    ! "a6(i, k)" is at longitude "rlon(i)", latitude "rlat(i)",
109    ! middle of layer "k".)
110
111    real, parameter:: amu = 1.6605402e-27 ! atomic mass unit, in kg
112
113    real, parameter:: Clx = 3.8e-9
114    ! (total chlorine content in the upper stratosphere)
115
116    integer k
117
118    !------------------------------------
119
120    !!print *, "Call sequence information: regr_pr_comb_coefoz"
121    call assert((/size(rlat), size(paprs, 1), size(pplay, 1)/) == klon, &
122         "regr_pr_comb_coefoz klon")
123    call assert((/size(paprs, 2) - 1, size(pplay, 2)/) == nbp_lev, &
124         "regr_pr_comb_coefoz nbp_lev")
125
126    !$omp master
127    if (is_mpi_root) call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)
128    !$omp end master
129
130    call regr_pr_time_av(ncid, (/"a2       ", "a4       ", "a6       ", &
131         "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), REAL(julien-1), &
132         press_in_edg, paprs, coefoz)
133    a2 = coefoz(:, :, 1)
134    a4_mass = coefoz(:, :, 2) * 48. / 29.
135
136    ! Compute "a6_mass" avoiding underflow, do not divide by 1e4
137    ! before dividing by molecular mass:
138    a6_mass = coefoz(:, :, 3) / (1e4 * 29. * amu)
139    ! (factor 1e4: conversion from cm2 to m2)
140
141    ! We can overwrite "coefoz(:, :, 1)", which was saved to "a2":
142    call regr_pr_int(ncid, "Sigma_Mob", julien, plev, pplay, top_value=0., &
143         v3=coefoz(:, :, 1))
144
145    ! Combine coefficients to get "c_Mob":
146    c_mob = (coefoz(:, :, 4) - a2 * coefoz(:, :, 5) &
147         - coefoz(:, :, 3) * coefoz(:, :, 1)) * 48. / 29. &
148         - a4_mass * coefoz(:, :, 6)
149
150    r_het_interm = coefoz(:, :, 7)
151    ! Heterogeneous chemistry is only at high latitudes:
152    forall (k = 1: nbp_lev)
153       where (abs(rlat) <= 45.) r_het_interm(:, k) = 0.
154    end forall
155    r_het_interm = r_het_interm * (Clx / 3.8e-9)**2
156
157    !$omp master
158    if (is_mpi_root) call nf95_close(ncid)
159    !$omp end master
160
161  end subroutine regr_pr_comb_coefoz
162
163end module regr_pr_comb_coefoz_m
Note: See TracBrowser for help on using the repository browser.