Changeset 1263 for LMDZ4


Ignore:
Timestamp:
Nov 17, 2009, 2:00:14 PM (15 years ago)
Author:
lguez
Message:

1) Reactivated ability to read ozone (that was deactivated because of
dependency on version of IOIPSL). Added ability to read a pressure
coordinate in Pa in "regr_lat_time_climoz".

2) Added the ability to read a second ozone climatology, corresponding to
daylight ozone:

-- "read_climoz" is now an integer variable, instead of a logical
variable.

-- Added argument "read_climoz" to "phys_state_var_init",
"phys_output_open" and "regr_lat_time_climoz".

-- Created new variable "ozone_daylight" for "hist*.nc" output files.

-- Added a third dimension to variable "wo" in module
"phys_state_var_mod" and variable "POZON" in "radlwsw": index 1 for
average day-night ozone, index 2 for daylight ozone.

-- Added a fourth dimension to variables "o3_in", "o3_regr_lat" and
"o3_out" in "regr_lat_time_climoz": index 1 for average day-night
ozone, index 2 for daylight ozone.

-- In "physiq", moved call to "conf_phys" before call to
"phys_state_var_init". Thus, "conf_phys" is now inside the block "if
(first)" instead of "IF (debut)". There were definitions of "bl95_b0"
and "bl95_b1" that were useless because the variables were overwritten
by "conf_phys". Removed those definitions.

-- In "radlwsw", we pass the average day-night ozone to "LW_LMDAR4"
and the daylight ozone, if we have it, to "SW_LMDAR4" or
"SW_AEROAR4". If we do not have a specific field for daylight ozone
then "SW_LMDAR4" or "SW_AEROAR4" just get the average day-night ozone.

-- "regr_lat_time_climoz" now manages latitudes where the input ozone
field is missing at all levels (polar night).

-- Encapsulated "radlwsw" in a module.

3) Modifications to make sequential and parallel versions of
"create_etat0_limit" almost identical:

-- In "dyn3dpar/create_etat0_limit.F". No need to call
"phys_state_var_init", removed "use phys_state_var_mod" statement. No
need for "clesphys.h", removed "include" statement.

-- In "dyn3dpar/etat0_netcdf.F". Added argument "tau_ratqs" in call to
"conf_phys" (this bug was already corrected in "dyn3d"). Moved call to
"inifilr" after call to "infotrac_init" (as in "dyn3d").

4) Other peripheral modifications:

-- Added procedures "nf95_get_att" and "nf95_def_var_scalar" in
NetCDF95 interface. Overloaded "nf95_put_var" with three more
procedures: "nf95_put_var_FourByteReal", "nf95_put_var_FourByteInt",
"nf95_put_var_1D_FourByteInt".

-- Overloaded "regr1_step_av" with one more procedure:
"regr14_step_av". Overloaded "regr3_lint" with one more procedure:
"regr34_lint".

-- Corrected call to "Init_Phys_lmdz" in "dyn3d/create_etat0_limit.F":
the last argument should be an array, not a scalar.

-- Encapsulated "conf_phys" in a module.

-- Splitted module "regr_pr" into "regr_pr_av_m" and "regr_pr_int_m".

5) Tests:

This revision was compared to revision 1259, with optimization options
"debug" and "dev", parallelization options "none", "mpi", "omp" and
"mpi_omp", 1 and 2 MPI processes, 1 and 2 OpenMP threads, with the
compiler "FORTRAN90/SX Version 2.0 for SX-8". Both programs
"create_etat0_limit" and "gcm" were tested. In all cases,
parallelization does not change the results. With "read_climoz = 0" in
the ".def" files, the results of revision 1259 and of this revision
are the same.

Location:
LMDZ4/branches/LMDZ4-dev/libf
Files:
2 added
17 edited
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/netcdf95.F90

    r1157 r1263  
    2525  ! criticisms for some (not all) procedures.
    2626
     27  ! "nf95_get_att" is more secure than "nf90_get_att" because it
     28  ! checks that the "values" argument is long enough and removes the
     29  ! null terminator, if any.
     30
    2731  ! This module replaces some of the official NetCDF procedures.
    2832  ! This module also provides the procedures "handle_err" and "nf95_gw_var".
    … …  
    3539  use nf95_gw_var_m
    3640  use nf95_put_att_m
     41  use nf95_get_att_m
    3742  use simple
    3843  use handle_err_m

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_def_var_m.F90

    r1157 r1263  
    11! $Id$
    22module nf95_def_var_m
     3
     4  ! The generic procedure name "nf90_def_var" applies to
     5  ! "nf90_def_var_Scalar" but we cannot apply the generic procedure name
     6  ! "nf95_def_var" to "nf95_def_var_scalar" because of the additional
     7  ! optional argument.
     8  ! "nf95_def_var_scalar" cannot be distinguished from "nf95_def_var_oneDim".
    39
    410  implicit none
    … …  
    915
    1016  private
    11   public nf95_def_var
     17  public nf95_def_var, nf95_def_var_scalar
    1218
    1319contains
     20
     21  subroutine nf95_def_var_scalar(ncid, name, xtype, varid, ncerr)
     22
     23    use netcdf, only: nf90_def_var
     24    use handle_err_m, only: handle_err
     25
     26    integer,               intent( in) :: ncid
     27    character (len = *),   intent( in) :: name
     28    integer,               intent( in) :: xtype
     29    integer,               intent(out) :: varid
     30    integer, intent(out), optional:: ncerr
     31
     32    ! Variable local to the procedure:
     33    integer ncerr_not_opt
     34
     35    !-------------------
     36
     37    ncerr_not_opt = nf90_def_var(ncid, name, xtype, varid)
     38    if (present(ncerr)) then
     39       ncerr = ncerr_not_opt
     40    else
     41       call handle_err("nf95_def_var_scalar " // name, ncerr_not_opt, ncid)
     42    end if
     43
     44  end subroutine nf95_def_var_scalar
     45
     46  !***********************
    1447
    1548  subroutine nf95_def_var_oneDim(ncid, name, xtype, dimids, varid, ncerr)

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/nf95_put_var_m.F90

    r1157 r1263  
    55
    66  interface nf95_put_var
    7      module procedure nf95_put_var_1D_FourByteReal, &
     7     module procedure nf95_put_var_FourByteReal, nf95_put_var_FourByteInt, &
     8          nf95_put_var_1D_FourByteReal, nf95_put_var_1D_FourByteInt, &
    89          nf95_put_var_2D_FourByteReal, nf95_put_var_3D_FourByteReal, &
    910          nf95_put_var_4D_FourByteReal
    … …  
    1920contains
    2021
     22  subroutine nf95_put_var_FourByteReal(ncid, varid, values, start, ncerr)
     23
     24    use netcdf, only: nf90_put_var
     25    use handle_err_m, only: handle_err
     26
     27    integer, intent( in) :: ncid, varid
     28    real, intent( in) :: values
     29    integer, dimension(:), optional, intent( in) :: start
     30    integer, intent(out), optional:: ncerr
     31
     32    ! Variable local to the procedure:
     33    integer ncerr_not_opt
     34
     35    !-------------------
     36
     37    ncerr_not_opt = nf90_put_var(ncid, varid, values, start)
     38    if (present(ncerr)) then
     39       ncerr = ncerr_not_opt
     40    else
     41       call handle_err("nf95_put_var_FourByteReal", ncerr_not_opt, ncid, &
     42            varid)
     43    end if
     44
     45  end subroutine nf95_put_var_FourByteReal
     46
     47  !***********************
     48
     49  subroutine nf95_put_var_FourByteInt(ncid, varid, values, start, ncerr)
     50
     51    use netcdf, only: nf90_put_var
     52    use handle_err_m, only: handle_err
     53
     54    integer, intent( in) :: ncid, varid
     55    integer, intent( in) :: values
     56    integer, dimension(:), optional, intent( in) :: start
     57    integer, intent(out), optional:: ncerr
     58
     59    ! Variable local to the procedure:
     60    integer ncerr_not_opt
     61
     62    !-------------------
     63
     64    ncerr_not_opt = nf90_put_var(ncid, varid, values, start)
     65    if (present(ncerr)) then
     66       ncerr = ncerr_not_opt
     67    else
     68       call handle_err("nf95_put_var_FourByteInt", ncerr_not_opt, ncid, &
     69            varid)
     70    end if
     71
     72  end subroutine nf95_put_var_FourByteInt
     73
     74  !***********************
     75
    2176  subroutine nf95_put_var_1D_FourByteReal(ncid, varid, values, start, count, &
    2277       stride, map, ncerr)
    … …  
    45100
    46101  end subroutine nf95_put_var_1D_FourByteReal
     102
     103  !***********************
     104
     105  subroutine nf95_put_var_1D_FourByteInt(ncid, varid, values, start, count, &
     106       stride, map, ncerr)
     107
     108    use netcdf, only: nf90_put_var
     109    use handle_err_m, only: handle_err
     110
     111    integer,                         intent(in) :: ncid, varid
     112    integer, intent(in) :: values(:)
     113    integer, dimension(:), optional, intent(in) :: start, count, stride, map
     114    integer, intent(out), optional:: ncerr
     115
     116    ! Variable local to the procedure:
     117    integer ncerr_not_opt
     118
     119    !-------------------
     120
     121    ncerr_not_opt = nf90_put_var(ncid, varid, values, start, count, stride, &
     122         map)
     123    if (present(ncerr)) then
     124       ncerr = ncerr_not_opt
     125    else
     126       call handle_err("nf95_put_var_1D_FourByteInt", ncerr_not_opt, ncid, &
     127            varid)
     128    end if
     129
     130  end subroutine nf95_put_var_1D_FourByteInt
    47131
    48132  !***********************

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90

    r1157 r1263  
    1717     ! The difference between the procedures is the rank of the first argument.
    1818
    19      module procedure regr11_step_av, regr12_step_av, regr13_step_av
     19     module procedure regr11_step_av, regr12_step_av, regr13_step_av, &
     20          regr14_step_av
    2021  end interface
    2122
    … …  
    203204  end function regr13_step_av
    204205
     206  !********************************************
     207
     208  function regr14_step_av(vs, xs, xt) result(vt)
     209
     210    ! "vs" has rank 4.
     211
     212    use assert_eq_m, only: assert_eq
     213    use assert_m, only: assert
     214    use interpolation, only: locate
     215
     216    real, intent(in):: vs(:, :, :, :) ! values of steps on the source grid
     217    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
     218
     219    real, intent(in):: xs(:)
     220    ! (edges of steps on the source grid, in strictly increasing order)
     221
     222    real, intent(in):: xt(:)
     223    ! (edges of cells of the target grid, in strictly increasing order)
     224
     225    real vt(size(xt) - 1, size(vs, 2), size(vs, 3), size(vs, 4))
     226    ! (average values on the target grid)
     227    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
     228
     229    ! Variables local to the procedure:
     230    integer is, it, ns, nt
     231    real left_edge
     232
     233    !---------------------------------------------
     234
     235    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_step_av ns")
     236    nt = size(xt) - 1
     237
     238    ! Quick check on sort order:
     239    call assert(xs(1) < xs(2), "regr14_step_av xs bad order")
     240    call assert(xt(1) < xt(2), "regr14_step_av xt bad order")
     241
     242    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
     243         "regr14_step_av extrapolation")
     244
     245    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
     246    do it = 1, nt
     247       ! 1 <= is <= ns
     248       ! xs(is) <= xt(it) < xs(is + 1)
     249       ! Compute "vt(it, :, :, :)":
     250       left_edge = xt(it)
     251       vt(it, :, :, :) = 0.
     252       do while (xs(is + 1) < xt(it + 1))
     253          ! 1 <= is <= ns - 1
     254          vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - left_edge) &
     255               * vs(is, :, :, :)
     256          is = is + 1
     257          left_edge = xs(is)
     258       end do
     259       ! 1 <= is <= ns
     260       vt(it, :, :, :) = (vt(it, :, :, :) + (xt(it + 1) - left_edge) &
     261            * vs(is, :, :, :)) / (xt(it + 1) - xt(it))
     262       if (xs(is + 1) == xt(it + 1)) is = is + 1
     263       ! 1 <= is <= ns .or. it == nt
     264    end do
     265
     266  end function regr14_step_av
     267
    205268end module regr1_step_av_m

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/regr3_lint_m.F90

    r1157 r1263  
    1111     ! input array.
    1212     ! The difference betwwen the procedures is the rank of the first argument.
    13      module procedure regr33_lint
     13     module procedure regr33_lint, regr34_lint
    1414  end interface
    1515
    … …  
    5757  end function regr33_lint
    5858
     59  !*********************************************************
     60
     61  function regr34_lint(vs, xs, xt) result(vt)
     62
     63    ! "vs" has rank 4.
     64
     65    use assert_eq_m, only: assert_eq
     66    use interpolation, only: hunt
     67
     68    real, intent(in):: vs(:, :, :, :)
     69    ! (values of the function at source points "xs")
     70
     71    real, intent(in):: xs(:)
     72    ! (abscissas of points in source grid, in strictly monotonic order)
     73
     74    real, intent(in):: xt(:)
     75    ! (abscissas of points in target grid)
     76
     77    real vt(size(vs, 1), size(vs, 2), size(xt), size(vs, 4))
     78    ! (values of the function on the target grid)
     79
     80    ! Variables local to the procedure:
     81    integer is, it, ns
     82    integer is_b ! "is" bound between 1 and "ns - 1"
     83
     84    !--------------------------------------
     85
     86    ns = assert_eq(size(vs, 3), size(xs), "regr34_lint ns")
     87
     88    is = -1 ! go immediately to bisection on first call to "hunt"
     89
     90    do it = 1, size(xt)
     91       call hunt(xs, xt(it), is)
     92       is_b = min(max(is, 1), ns - 1)
     93       vt(:, :, it, :) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b, :) &
     94            + (xt(it) - xs(is_b)) * vs(:, :, is_b+1, :)) &
     95            / (xs(is_b+1) - xs(is_b))
     96    end do
     97
     98  end function regr34_lint
     99
    59100end module regr3_lint_m

  • LMDZ4/branches/LMDZ4-dev/libf/bibio/simple.F90

    r1157 r1263  
    118118    ! unknown in the calling procedure, before the call.
    119119    ! Here we use a better solution: a pointer argument array.
    120     ! This procedure associates and defines '"dimids" if it is present.
     120    ! This procedure associates and defines "dimids" if it is present.
    121121
    122122    use netcdf, only: nf90_inquire_variable, nf90_max_var_dims

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/create_etat0_limit.F

    r1222 r1263  
    4747      END IF
    4848
    49       CALL Init_Phys_lmdz(iim,jjp1,llm,1,(jjm-1)*iim+2)
     49      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     50      PRINT *,'---> klon=',klon
    5051      call InitComgeomphy
    5152

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3d/etat0_netcdf.F

    r1245 r1263  
    1515      USE filtreg_mod
    1616      use regr_lat_time_climoz_m, only: regr_lat_time_climoz
     17      use conf_phys_m, only: conf_phys
    1718#endif
    1819!#endif of #ifdef CPP_EARTH
    … …  
    143144      REAL      :: solarlong0
    144145      real :: seuil_inversion
    145       logical  read_climoz ! read ozone climatology
     146
     147      integer  read_climoz ! read ozone climatology
     148C     Allowed values are 0, 1 and 2
     149C     0: do not read an ozone climatology
     150C     1: read a single ozone climatology that will be used day and night
     151C     2: read two ozone climatologies, the average day and night
     152C     climatology and the daylight climatology
    146153
    147154      !
    … …  
    179186     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
    180187
    181 
    182188! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
    183189      co2_ppm0 = co2_ppm
    … …  
    194200
    195201      CALL inifilr()
    196       CALL phys_state_var_init()
     202      CALL phys_state_var_init(read_climoz)
    197203      !
    198204      latfi(1) = ASIN(1.0)
    … …  
    419425      !
    420426
    421       if (read_climoz) call regr_lat_time_climoz ! ozone climatology
     427!     Ozone climatology:
     428      if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz)
    422429
    423430      varname = 'tsol'

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/create_etat0_limit.F

    r1222 r1263  
    99       USE mod_phys_lmdz_para
    1010       USE mod_const_mpi
    11        USE phys_state_var_mod     
    1211       USE infotrac
    1312#ifdef CPP_IOIPSL
    … …  
    3837#include "indicesol.h"
    3938#include  "control.h"
    40 #include "clesphys.h"
    4139      REAL :: masque(iip1,jjp1)
    4240!      REAL :: pctsrf(iim*(jjm-1)+2, nbsrf)
    … …  
    6260     &                 for 1 process and 1 task')
    6361      ENDIF
    64       CALL phys_state_var_init
    6562      call InitComgeomphy
    66      
     63
    6764#ifdef CPP_IOIPSL
    6865      call ioconf_calendar('360d')

  • LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/etat0_netcdf.F

    r1227 r1263  
    1515      USE filtreg_mod
    1616      use regr_lat_time_climoz_m, only: regr_lat_time_climoz
     17      use conf_phys_m, only: conf_phys
    1718#endif
    1819!#endif of #ifdef CPP_EARTH
    … …  
    132133      REAL                 :: bl95_b0, bl95_b1
    133134      real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
     135      real                 :: tau_ratqs
    134136      integer              :: iflag_cldcon
    135137      integer              :: iflag_ratqs
    … …  
    142144      REAL      :: solarlong0
    143145      real :: seuil_inversion
    144       logical  read_climoz ! read ozone climatology
     146
     147      integer  read_climoz ! read ozone climatology
     148C     Allowed values are 0, 1 and 2
     149C     0: do not read an ozone climatology
     150C     1: read a single ozone climatology that will be used day and night
     151C     2: read two ozone climatologies, the average day and night
     152C     climatology and the daylight climatology
    145153
    146154      !
    … …  
    170178     &                 fact_cldcon, facttemps,ok_newmicro,iflag_radia,  &
    171179     &                 iflag_cldcon,                                    &
    172      &                 iflag_ratqs,ratqsbas,ratqshaut,                  &
     180     &                 iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,        &
    173181     &                 ok_ade, ok_aie, aerosol_couple,                  &
    174182     &                 flag_aerosol, new_aod,                           &
    … …  
    187195      CALL inigeom()
    188196
    189       CALL inifilr()
    190197! Initialisation pour traceurs
    191198      call infotrac_init
    192199      ALLOCATE(q3d(iip1, jjp1, llm, nqtot))
    193 !     CALL phys_state_var_init()
     200
     201      CALL inifilr()
     202      CALL phys_state_var_init(read_climoz)
    194203      !
    195204      latfi(1) = ASIN(1.0)
    … …  
    416425      !
    417426
    418       if (read_climoz) call regr_lat_time_climoz ! ozone climatology
     427!     Ozone climatology:
     428      if (read_climoz >= 1) call regr_lat_time_climoz(read_climoz)
    419429
    420430      varname = 'tsol'

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90

    r1259 r1263  
    55!
    66!
     7module conf_phys_m
     8
     9   implicit none
     10
     11contains
    712
    813  subroutine conf_phys(ok_journe, ok_mensuel, ok_instan, ok_hf, &
    … …  
    2227   USE surface_data
    2328   USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl
    24 
    25    implicit none
    2629
    2730 include "conema3.h"
    … …  
    152155  LOGICAL,SAVE      :: carbon_cycle_cpl_omp
    153156
    154   logical, intent(out):: read_climoz ! read ozone climatology, OpenMP shared
    155 !
     157  integer, intent(out):: read_climoz ! read ozone climatology, OpenMP shared
     158  ! Allowed values are 0, 1 and 2
     159  ! 0: do not read an ozone climatology
     160  ! 1: read a single ozone climatology that will be used day and night
     161  ! 2: read two ozone climatologies, the average day and night
     162  ! climatology and the daylight climatology
    156163
    157164!$OMP MASTER
    … …  
    13111318  call getin('ecrit_LES', ecrit_LES_omp)
    13121319!
    1313   read_climoz = .false. ! default value
     1320  read_climoz = 0 ! default value
    13141321  call getin('read_climoz', read_climoz)
    13151322
    … …  
    16011608  end subroutine conf_phys
    16021609
     1610end module conf_phys_m
    16031611!
    16041612!#################################################################

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_mod.F90

    r1258 r1263  
    390390  type(ctrl_out),save :: o_rhum         = ctrl_out((/ 2, 10, 10, 10, 10 /),'rhum')
    391391  type(ctrl_out),save :: o_ozone        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone')
     392  type(ctrl_out),save :: o_ozone_light        = ctrl_out((/ 2, 10, 10, 10, 10 /),'ozone_daylight')
    392393  type(ctrl_out),save :: o_upwd         = ctrl_out((/ 2, 10, 10, 10, 10 /),'upwd')
    393394  type(ctrl_out),save :: o_dtphy        = ctrl_out((/ 2, 10, 10, 10, 1 /),'dtphy')
    … …  
    483484 
    484485  SUBROUTINE phys_output_open(jjmp1,nlevSTD,clevSTD,nbteta, &
    485                               ctetaSTD,dtime, ok_veget, &
    486                               type_ocean, iflag_pbl,ok_mensuel,ok_journe, &
    487                               ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)   
     486       ctetaSTD,dtime, ok_veget, &
     487       type_ocean, iflag_pbl,ok_mensuel,ok_journe, &
     488       ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, read_climoz)
    488489
    489490  USE iophy
    … …  
    506507  logical                               :: ok_mensuel, ok_journe, ok_hf, ok_instan
    507508  logical                               :: ok_LES,ok_ade,ok_aie
     509
     510  integer, intent(in)::  read_climoz ! read ozone climatology
     511  !     Allowed values are 0, 1 and 2
     512  !     0: do not read an ozone climatology
     513  !     1: read a single ozone climatology that will be used day and night
     514  !     2: read two ozone climatologies, the average day and night
     515  !     climatology and the daylight climatology
     516
    508517  real                                  :: dtime
    509518  integer                               :: idayref
    … …  
    979988 CALL histdef3d(iff,o_rhum%flag,o_rhum%name, "Relative humidity", "-")
    980989 CALL histdef3d(iff,o_ozone%flag,o_ozone%name, "Ozone mole fraction", "-")
     990 if (read_climoz == 2) &
     991      CALL histdef3d(iff,o_ozone_light%flag,o_ozone_light%name, &
     992      "Daylight ozone mole fraction", "-")
    981993 CALL histdef3d(iff,o_dtphy%flag,o_dtphy%name, "Physics dT", "K/s")
    982994 CALL histdef3d(iff,o_dqphy%flag,o_dqphy%name, "Physics dQ", "(kg/kg)/s")

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h

    r1258 r1263  
    10041004      IF (o_ozone%flag(iff)<=lev_files(iff)) THEN
    10051005         CALL histwrite_phy(nid_files(iff), o_ozone%name, itau_w,
    1006      $        wo * dobson_u * 1e3 / zmasse / rmo3 * rmd)
     1006     $        wo(:, :, 1) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
     1007      ENDIF
     1008
     1009      IF (o_ozone_light%flag(iff)<=lev_files(iff) .and.
     1010     $     read_climoz == 2) THEN
     1011         CALL histwrite_phy(nid_files(iff), o_ozone_light%name, itau_w,
     1012     $        wo(:, :, 2) * dobson_u * 1e3 / zmasse / rmo3 * rmd)
    10071013      ENDIF
    10081014

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90

    r1227 r1263  
    203203!$OMP THREADPRIVATE(albsol1,albsol2)
    204204
    205       REAL, ALLOCATABLE, SAVE:: wo(:,:)
     205      REAL, ALLOCATABLE, SAVE:: wo(:, :, :)
    206206      ! column-density of ozone in a layer, in kilo-Dobsons
     207      ! Third dimension has size 1 or 2.
     208      ! "wo(:, :, 1)" is for the average day-night field,
     209      ! "wo(:, :, 2)" is for daylight time.
    207210      !$OMP THREADPRIVATE(wo)
    208211
    … …  
    284287
    285288!======================================================================
    286 SUBROUTINE phys_state_var_init
     289SUBROUTINE phys_state_var_init(read_climoz)
    287290use dimphy
    288291use aero_mod
    289292IMPLICIT NONE
     293
     294integer, intent(in)::  read_climoz
     295! read ozone climatology
     296! Allowed values are 0, 1 and 2
     297! 0: do not read an ozone climatology
     298! 1: read a single ozone climatology that will be used day and night
     299! 2: read two ozone climatologies, the average day and night
     300! climatology and the daylight climatology
     301
    290302#include "indicesol.h"
    291303#include "control.h"
    … …  
    365377      ALLOCATE(paire_ter(klon))
    366378      ALLOCATE(albsol1(klon), albsol2(klon))
    367       ALLOCATE(wo(klon,klev))
     379
     380      if (read_climoz <= 1) then
     381         ALLOCATE(wo(klon,klev, 1))
     382      else
     383         ! read_climoz == 2
     384         ALLOCATE(wo(klon,klev, 2))
     385      end if
     386     
    368387      ALLOCATE(clwcon0(klon,klev),rnebcon0(klon,klev))
    369388      ALLOCATE(heat(klon,klev), heat0(klon,klev))

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F

    r1259 r1263  
    1313
    1414      USE ioipsl, only: histbeg, histvert, histdef, histend, histsync,
    15      $     histwrite, ju2ymds, ymds2ju
    16 !JG     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
     15     $     histwrite, ju2ymds, ymds2ju, ioget_year_len
    1716      USE comgeomphy
    1817      USE phys_cal_mod
    … …  
    3332      USE phys_output_mod
    3433      use open_climoz_m, only: open_climoz ! ozone climatology from a file
    35       use regr_pr, only: regr_pr_av
     34      use regr_pr_av_m, only: regr_pr_av
    3635      use netcdf95, only: nf95_close
    3736      use mod_phys_lmdz_mpi_data, only: is_mpi_root
    3837      USE aero_mod
    3938      use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
     39      use conf_phys_m, only: conf_phys
     40      use radlwsw_m, only: radlwsw
    4041
    4142      IMPLICIT none
    … …  
    10931094      integer iunit
    10941095
    1095       logical, save::  read_climoz ! read ozone climatology from a file
    1096       integer, save:: ncid_climoz ! NetCDF file containing ozone climatology
     1096      integer, save::  read_climoz ! read ozone climatology
     1097C     Allowed values are 0, 1 and 2
     1098C     0: do not read an ozone climatology
     1099C     1: read a single ozone climatology that will be used day and night
     1100C     2: read two ozone climatologies, the average day and night
     1101C     climatology and the daylight climatology
     1102
     1103      integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
    10971104
    10981105      real, pointer, save:: press_climoz(:)
    1099 !     edges of pressure intervals for ozone climatology, in Pa, in strictly
     1106!     edges of pressure intervals for ozone climatologies, in Pa, in strictly
    11001107!     ascending order
    11011108
    1102       integer,save:: co3i = 0 ! time index in NetCDF file of current ozone field
     1109      integer, save:: co3i = 0
     1110!     time index in NetCDF file of current ozone fields
    11031111c$OMP THREADPRIVATE(co3i)
    11041112
    11051113      integer ro3i
    1106 !     required time index in NetCDF file for the ozone field, between 1 and 360
     1114!     required time index in NetCDF file for the ozone fields, between 1
     1115!     and 360
    11071116
    11081117#include "YOMCST.h"
    … …  
    11651174      print*, 'Allocation des variables locales et sauvegardees'
    11661175      call phys_local_var_init
    1167       call phys_state_var_init
     1176c     appel a la lecture du run.def physique
     1177      call conf_phys(ok_journe, ok_mensuel,
     1178     .     ok_instan, ok_hf,
     1179     .     ok_LES,
     1180     .     solarlong0,seuil_inversion,
     1181     .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
     1182     .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
     1183     .     ok_ade, ok_aie, aerosol_couple,
     1184     .     flag_aerosol, new_aod,
     1185     .     bl95_b0, bl95_b1,
     1186     .     iflag_thermals,nsplit_thermals,tau_thermals,
     1187     .     iflag_thermals_ed,iflag_thermals_optflux,
     1188c     nv flags pour la convection et les poches froides
     1189     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz)
     1190      call phys_state_var_init(read_climoz)
    11681191      print*, '================================================='
    11691192
    … …  
    11811204        first=.false.
    11821205
    1183       endif  ! fisrt
     1206      endif  ! first
    11841207
    11851208       modname = 'physiq'
    … …  
    12091232c
    12101233       IF (debut) THEN
    1211 C
    12121234!rv
    12131235cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
    … …  
    12201242         rain_con(:)=0.
    12211243         snow_con(:)=0.
    1222          bl95_b0=0.
    1223          bl95_b1=0.
    12241244         topswai(:)=0.
    12251245         topswad(:)=0.
    … …  
    12471267         IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
    12481268c
    1249 c appel a la lecture du run.def physique
    1250 c
    1251          call conf_phys(ok_journe, ok_mensuel,
    1252      .                  ok_instan, ok_hf,
    1253      .                  ok_LES,
    1254      .                  solarlong0,seuil_inversion,
    1255      .                  fact_cldcon, facttemps,ok_newmicro,iflag_radia,
    1256      .         iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
    1257      .                  ok_ade, ok_aie, aerosol_couple,
    1258      .                  flag_aerosol, new_aod,
    1259      .                  bl95_b0, bl95_b1,
    1260      .                  iflag_thermals,nsplit_thermals,tau_thermals,
    1261      .                  iflag_thermals_ed,iflag_thermals_optflux,
    1262 cnv flags pour la convection et les poches froides
    1263      .                   iflag_coupl,iflag_clos,iflag_wake, read_climoz)
    1264 
    12651269      print*,'iflag_coupl,iflag_clos,iflag_wake',
    12661270     .   iflag_coupl,iflag_clos,iflag_wake
    … …  
    14531457     &                        ctetaSTD,dtime,ok_veget,
    14541458     &                        type_ocean,iflag_pbl,ok_mensuel,ok_journe,
    1455      &                        ok_hf,ok_instan,ok_LES,ok_ade,ok_aie)
     1459     &              ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, read_climoz)
    14561460c$OMP END MASTER
    14571461c$OMP BARRIER
    … …  
    15441548
    15451549C$omp single
    1546       if (read_climoz) then
     1550      if (read_climoz >= 1) then
    15471551         call open_climoz(ncid_climoz, press_climoz)
    15481552      END IF
    … …  
    17111715c Prescrire l'ozone et calculer l'albedo sur l'ocean.
    17121716c
    1713       if (read_climoz) then
     1717      if (read_climoz >= 1) then
    17141718C        Ozone from a file
    17151719!        Update required ozone index:
    1716 !JG : ioget_year_len n'existe pas dans les versions officiels d'ioipsl
    1717 !JG         ro3i = int((days_elapsed + jh_cur - jh_1jan)
    1718 !JG     $        / ioget_year_len(year_cur) * 360.) + 1
    1719          CALL abort_gcm(modname,
    1720      $     'JG : read_climoz temporary desactivated',1)
    1721 
     1720         ro3i = int((days_elapsed + jh_cur - jh_1jan)
     1721     $        / ioget_year_len(year_cur) * 360.) + 1
    17221722         if (ro3i == 361) ro3i = 360
    17231723C        (This should never occur, except perhaps because of roundup
    … …  
    17251725         if (ro3i /= co3i) then
    17261726C           Update ozone field:
    1727             call regr_pr_av(ncid_climoz, "tro3", julien=ro3i,
    1728      &           press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1727            if (read_climoz == 1) then
     1728               call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i,
     1729     $              press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1730            else
     1731C              read_climoz == 2
     1732               call regr_pr_av(ncid_climoz,
     1733     $              (/"tro3         ", "tro3_daylight"/),
     1734     $              julien=ro3i, press_in_edg=press_climoz, paprs=paprs,
     1735     $              v3=wo)
     1736            end if
    17291737!           Convert from mole fraction of ozone to column density of ozone in a
    17301738!           cell, in kDU:
    1731             wo = wo * rmo3 / rmd * zmasse / dobson_u / 1e3
     1739            forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l)
     1740     $           * rmo3 / rmd * zmasse / dobson_u / 1e3
    17321741C           (By regridding ozone values for LMDZ only once every 360th of
    17331742C           year, we have already neglected the variation of pressure in one
    … …  
    17381747      elseif (MOD(itap-1,lmt_pas) == 0) THEN
    17391748C        Once per day, update ozone from Royer:
    1740          wo = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
     1749         wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
    17411750      ENDIF
    17421751c
    … …  
    28362845     $                          u,
    28372846     $                          v,
    2838      $                          wo,
     2847     $                          wo(:, :, 1),
    28392848     $                          q_seri,
    28402849     $                          zxtsol,
    … …  
    29192928     e        (kdlon,kflev,dist, rmu0, fract, solaire,
    29202929     e        paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
    2921      e        wo,
     2930     e        wo(:, :, 1),
    29222931     e        cldfra, cldemi, cldtau,
    29232932     s        heat,heat0,cool,cool0,radsol,albpla,
    … …  
    35293538!         close(97)
    35303539C$OMP MASTER
    3531          if (read_climoz) then
     3540         if (read_climoz >= 1) then
    35323541            if (is_mpi_root) then
    35333542               call nf95_close(ncid_climoz)

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F90

    r1246 r1263  
     1module radlwsw_m
     2
     3  IMPLICIT NONE
     4
     5contains
     6
    17SUBROUTINE radlwsw( &
    28   dist, rmu0, fract, &
    … …  
    2430
    2531  USE DIMPHY
    26 
    27   IMPLICIT NONE
     32  use assert_m, only: assert
     33
    2834  !======================================================================
    2935  ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
    … …  
    103109  REAL,    INTENT(in)  :: t(KLON,KLEV), q(KLON,KLEV)
    104110
    105   REAL, INTENT(in)::wo(KLON,KLEV)
     111  REAL, INTENT(in):: wo(:, :, :) ! dimension(KLON,KLEV, 1 or 2)
    106112  ! column-density of ozone in a layer, in kilo-Dobsons
     113  ! "wo(:, :, 1)" is for the average day-night field,
     114  ! "wo(:, :, 2)" is for daylight time.
    107115
    108116  LOGICAL, INTENT(in)  :: ok_ade, ok_aie                                 ! switches whether to use aerosol direct (indirect) effects or not
    … …  
    157165  REAL(KIND=8) PTAVE(kdlon,kflev)
    158166  REAL(KIND=8) PWV(kdlon,kflev), PQS(kdlon,kflev)
    159   real(kind=8) POZON(kdlon,kflev) ! mass fraction of ozone
     167
     168  real(kind=8) POZON(kdlon, kflev, size(wo, 3)) ! mass fraction of ozone
     169  ! "POZON(:, :, 1)" is for the average day-night field,
     170  ! "POZON(:, :, 2)" is for daylight time.
     171
    160172  REAL(KIND=8) PAER(kdlon,kflev,5)
    161173  REAL(KIND=8) PCLDLD(kdlon,kflev)
    … …  
    187199  real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
    188200
     201  call assert(size(wo, 1) == klon, size(wo, 2) == klev, "radlwsw wo")
    189202  ! initialisation
    190203  tauaero(:,:,:,:)=0.
    … …  
    245258        PWV(i,k) = MAX (q(iof+i,k), 1.0e-12)
    246259        PQS(i,k) = PWV(i,k)
    247         POZON(i,k) = wo(iof+i, k) * RG * dobson_u * 1e3 &
     260        POZON(i,k, :) = wo(iof+i, k, :) * RG * dobson_u * 1e3 &
    248261             / (paprs(iof+i, k) - paprs(iof+i, k+1))
    249262        PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k)
    … …  
    298311    IF (iflag_rrtm == 0) THEN
    299312       ! Old radiation scheme, used for AR4 runs
     313       ! average day-night ozone for longwave
    300314       CALL LW_LMDAR4(&
    301315            PPMB, PDP,&
    302316            PPSOL,PDT0,PEMIS,&
    303             PTL, PTAVE, PWV, POZON, PAER,&
     317            PTL, PTAVE, PWV, POZON(:, :, 1), PAER,&
    304318            PCLDLD,PCLDLU,&
    305319            PVIEW,&
    … …  
    309323            ZFLUP, ZFLDN, ZFLUP0,ZFLDN0)
    310324
    311 
     325       ! daylight ozone, if we have it, for short wave
    312326       IF (.NOT. new_aod) THEN
    313327          ! use old version
    … …  
    315329               PPMB, PDP, &
    316330               PPSOL, PALBD, PALBP,&
    317                PTAVE, PWV, PQS, POZON, PAER,&
     331               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
    318332               PCLDSW, PTAU, POMEGA, PCG,&
    319333               zheat, zheat0,&
    … …  
    330344               PPMB, PDP,&
    331345               PPSOL, PALBD, PALBP,&
    332                PTAVE, PWV, PQS, POZON, PAER,&
     346               PTAVE, PWV, PQS, POZON(:, :, size(wo, 3)), PAER,&
    333347               PCLDSW, PTAU, POMEGA, PCG,&
    334348               zheat, zheat0,&
    … …  
    438452 ENDDO ! j = 1, nb_gr
    439453
    440 
    441454END SUBROUTINE radlwsw
    442455
    443 
     456end module radlwsw_m

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_lat_time_climoz_m.F90

    r1239 r1263  
    1111contains
    1212
    13   subroutine regr_lat_time_climoz
     13  subroutine regr_lat_time_climoz(read_climoz)
    1414
    1515    ! "regr_lat_time_climoz" stands for "regrid latitude time
    … …  
    2323
    2424    ! If the input field has missing values, they must be signaled by
    25     ! the "missing_value" attribute. At each latitude and each date, the
    26     ! missing values are replaced by the lowest valid value above missing
    27     ! values.
     25    ! the "missing_value" attribute.
    2826
    2927    ! We assume that the input field is a step function of latitude
    … …  
    3432    ! The values of "rlatu" are taken to be the centers of intervals.
    3533
    36     ! Regridding in time is by linear interpolation.
    37     ! Monthly values are processed to get daily values, on the basis
    38     ! of a 360-day calendar.
     34    ! We assume that in the input file:
     35
     36    ! -- Latitude is in degrees.
     37
     38    ! -- Latitude and pressure are strictly monotonic (as all NetCDF
     39    ! coordinate variables should be).
     40
     41    ! -- The time coordinate is in ascending order (even though we do
     42    ! not use its values).
    3943    ! The input file may contain either values for 12 months or values
    4044    ! for 14 months.
    4145    ! If there are 14 months then we assume that we have (in that order):
    4246    ! December, January, February, ..., November, December, January
    43     ! We use the first December value to interpolate values between January
    44     ! 1st and mid-January.
     47
     48    ! -- Missing values are contiguous, at the bottom of
     49    ! the vertical domain and at the latitudinal boundaries.
     50
     51    ! If values are all missing at a given latitude and date, then we
     52    ! replace those missing values by values at the closest latitude,
     53    ! equatorward, with valid values.
     54    ! Then, at each latitude and each date, the missing values are replaced
     55    ! by the lowest valid value above missing values.
     56
     57    ! Regridding in time is by linear interpolation.
     58    ! Monthly values are processed to get daily values, on the basis
     59    ! of a 360-day calendar.
     60    ! If there are 14 months, we use the first December value to
     61    ! interpolate values between January 1st and mid-January.
    4562    ! We use the last January value to interpolate values between
    4663    ! mid-December and end of December.
    47     ! If there are only 12 months in the input file than we assume
     64    ! If there are only 12 months in the input file then we assume
    4865    ! periodicity for interpolation at the beginning and at the end of the
    4966    ! year.
    5067
    51     ! We assume that in the input file:
    52     ! -- latitude is in degrees;
    53     ! -- pressure is in hPa (even though we do not use pressure values
    54     ! here, we write the unit of pressure in the NetCDF header, and we
    55     ! will use the assumption later, when we regrid in pressure).
    56     ! -- latitude and pressure are strictly monotonic (as all NetCDF
    57     ! coordinate variables should be);
    58     ! -- time increases (even though we do not use values of the input
    59     ! time coordinate);
    60     ! -- missing values are vertically contiguous, at the bottom of
    61     ! the vertical domain;
    62     ! -- there is no missing value at the top level of the atmosphere.
    63 
    6468    use regr1_step_av_m, only: regr1_step_av
    6569    use regr3_lint_m, only: regr3_lint
    66     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
    67          nf95_put_var, nf95_gw_var
    68     use netcdf, only: nf90_nowrite, nf90_get_att, nf90_noerr
     70    use netcdf95, only: handle_err, nf95_close, nf95_get_att, nf95_gw_var, &
     71         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, nf95_open, &
     72         nf95_put_var
     73    use netcdf, only: nf90_get_att, nf90_get_var, nf90_noerr, nf90_nowrite
     74    use assert_m, only: assert
     75
     76    integer, intent(in):: read_climoz ! read ozone climatology
     77    ! Allowed values are 1 and 2
     78    ! 1: read a single ozone climatology that will be used day and night
     79    ! 2: read two ozone climatologies, the average day and night
     80    ! climatology and the daylight climatology
    6981
    7082    ! Variables local to the procedure:
    … …  
    7991    ! (for "pi")
    8092
    81     integer ncid_in, ncid_out ! NetCDF IDs for input and output files
    8293    integer n_plev ! number of pressure levels in the input data
    83     integer n_lat! number of latitudes in the input data
     94    integer n_lat ! number of latitudes in the input data
     95    integer n_month ! number of months in the input data
    8496
    8597    real, pointer:: latitude(:)
    … …  
    91103
    92104    real, pointer:: plev(:)
    93     ! pressure levels of input data, sorted in strictly ascending order
     105    ! pressure levels of input data, sorted in strictly ascending
     106    ! order, converted to hPa
    94107
    95108    logical desc_lat ! latitude in descending order in the input file
    96109    logical desc_plev ! pressure levels in descending order in the input file
    97110
    98     real, pointer:: o3_in(:, :, :) ! (n_lat, n_plev, 12 or 0:13)
    99     ! (ozone climatology from the input file)
    100     ! ("o3_in(j, k, l)" is at latitude "latitude(j)" and pressure
    101     ! level "plev(k)". "l" is between 1 and 12 or between 0 and 13)
     111    real, allocatable:: o3_in(:, :, :, :)
     112    ! (n_lat, n_plev, n_month, read_climoz)
     113    ! ozone climatologies from the input file
     114    ! "o3_in(j, k, :, :)" is at latitude "latitude(j)" and pressure
     115    ! level "plev(k)".
     116    ! Third dimension is month index, first value may be December or January.
     117    ! "o3_in(:, :, :, 1)" is for the day- night average, "o3_in(:, :, :, 2)"
     118    ! is for daylight.
    102119
    103120    real missing_value
    104121
    105     real, allocatable:: o3_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13)
    106     ! (mean of "o3_in" over a latitude interval of LMDZ)
    107     ! (First dimension is latitude interval.
    108     ! The latitude interval for "o3_regr_lat(j,:, :)" contains "rlatu(j)".
     122    real, allocatable:: o3_regr_lat(:, :, :, :)
     123    ! (jjm + 1, n_plev, 0:13, read_climoz)
     124    ! mean of "o3_in" over a latitude interval of LMDZ
     125    ! First dimension is latitude interval.
     126    ! The latitude interval for "o3_regr_lat(j,:, :, :)" contains "rlatu(j)".
    109127    ! If "j" is between 2 and "jjm" then the interval is:
    110128    ! [rlatv(j), rlatv(j-1)]
    … …  
    114132    ! [- pi / 2, rlatv(jjm)]
    115133    ! respectively.
    116     ! "o3_regr_lat(:, k, :)" is for pressure level "plev(k)".
    117     ! Last dimension is month number.)
    118 
    119     real, allocatable:: o3_out(:, :, :) ! (jjm + 1, n_plev, 360)
    120     ! (regridded ozone climatology)
    121     ! ("o3_out(j, k, l)" is at latitude "rlatu(j)", pressure
     134    ! "o3_regr_lat(:, k, :, :)" is for pressure level "plev(k)".
     135    ! Third dimension is month number, 1 for January.
     136    ! "o3_regr_lat(:, :, :, 1)" is average day and night,
     137    ! "o3_regr_lat(:, :, :, 2)" is for daylight.
     138
     139    real, allocatable:: o3_out(:, :, :, :)
     140    ! (jjm + 1, n_plev, 360, read_climoz)
     141    ! regridded ozone climatology
     142    ! "o3_out(j, k, l, :)" is at latitude "rlatu(j)", pressure
    122143    ! level "plev(k)" and date "January 1st 0h" + "tmidday(l)", in a
    123     ! 360-day calendar.)
    124 
    125     integer j, k, l
     144    ! 360-day calendar.
     145    ! "o3_out(:, :, :, 1)" is average day and night,
     146    ! "o3_out(:, :, :, 2)" is for daylight.
     147
     148    integer j, k, l,m
    126149
    127150    ! For NetCDF:
    128     integer varid_in, varid_out, varid_plev, varid_time, varid, ncerr
     151    integer ncid_in, ncid_out ! IDs for input and output files
     152    integer varid_plev, varid_time, varid, ncerr, dimid
     153    character(len=80) press_unit ! pressure unit
     154
     155    integer varid_in(read_climoz), varid_out(read_climoz)
     156    ! index 1 is for average ozone day and night, index 2 is for
     157    ! daylight ozone.
    129158
    130159    real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * l, l = 0, 13)/)
    … …  
    141170
    142171    print *, "Call sequence information: regr_lat_time_climoz"
     172    call assert(read_climoz == 1 .or. read_climoz == 2, "regr_lat_time_climoz")
    143173
    144174    call nf95_open("climoz.nc", nf90_nowrite, ncid_in)
    … …  
    172202    desc_plev = plev(1) > plev(n_plev)
    173203    if (desc_plev) plev = plev(n_plev:1:-1)
     204    call nf95_get_att(ncid_in, varid, "units", press_unit)
     205    if (press_unit == "Pa") then
     206       ! Convert to hPa:
     207       plev = plev / 100.
     208    elseif (press_unit /= "hPa") then
     209       print *, "regr_lat_time_climoz: the only recognized units are Pa " &
     210            // "and hPa."
     211       stop 1
     212    end if
    174213
    175214    ! Create the output file and get the variable IDs:
    … …  
    183222    deallocate(plev) ! pointer
    184223
    185     call nf95_inq_varid(ncid_in, "tro3", varid_in)
    186     call nf95_gw_var(ncid_in, varid_in, o3_in)
    187     if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :)
    188     if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :)
    189 
    190     ncerr = nf90_get_att(ncid_in, varid_in, "missing_value", missing_value)
    191     if (ncerr == nf90_noerr) then
    192        do l = 1, size(o3_in, 3)
    193           do j = 1, n_lat
    194              ! Find missing values, starting from top of atmosphere
    195              ! and going down.
    196              ! We assume that the highest level contains no missing value.
    197              k = 2
    198              do
    199                 if (o3_in(j, k, l) == missing_value .or. k == n_plev) exit
    200                 k = k + 1
     224    ! Get the  number of months:
     225    call nf95_inq_dimid(ncid_in, "time", dimid)
     226    call nf95_inquire_dimension(ncid_in, dimid, len=n_month)
     227
     228    allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
     229
     230    call nf95_inq_varid(ncid_in, "tro3", varid_in(1))
     231    ncerr = nf90_get_var(ncid_in, varid_in(1), o3_in(:, :, :, 1))
     232    call handle_err("regr_lat_time_climoz nf90_get_var tro3", ncerr, ncid_in)
     233
     234    if (read_climoz == 2) then
     235       call nf95_inq_varid(ncid_in, "tro3_daylight", varid_in(2))
     236       ncerr = nf90_get_var(ncid_in, varid_in(2), o3_in(:, :, :, 2))
     237       call handle_err("regr_lat_time_climoz nf90_get_var tro3_daylight", &
     238            ncerr, ncid_in, varid_in(2))
     239    end if
     240
     241    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :, :)
     242    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :, :)
     243
     244    do m = 1, read_climoz
     245       ncerr = nf90_get_att(ncid_in, varid_in(m), "missing_value", &
     246            missing_value)
     247       if (ncerr == nf90_noerr) then
     248          do l = 1, n_month
     249             ! Take care of latitudes where values are all missing:
     250
     251             ! Next to the south pole:
     252             j = 1
     253             do while (o3_in(j, 1, l, m) == missing_value)
     254                j = j + 1
    201255             end do
    202              ! Replace missing values with the valid value at the
    203              ! lowest level above missing values:
    204              if (o3_in(j, k, l) == missing_value) &
    205                   o3_in(j, k:n_plev, l) = o3_in(j, k-1, l)
     256             if (j > 1) o3_in(:j-1, :, l, m) = &
     257                  spread(o3_in(j, :, l, m), dim=1, ncopies=j-1)
     258             
     259             ! Next to the north pole:
     260             j = n_lat
     261             do while (o3_in(j, 1, l, m) == missing_value)
     262                j = j - 1
     263             end do
     264             if (j < n_lat) o3_in(j+1:, :, l, m) = &
     265                  spread(o3_in(j, :, l, m), dim=1, ncopies=n_lat-j)
     266
     267             ! Take care of missing values at high pressure:
     268             do j = 1, n_lat
     269                ! Find missing values, starting from top of atmosphere
     270                ! and going down.
     271                ! We have already taken care of latitudes full of
     272                ! missing values so the highest level has a valid value.
     273                k = 2
     274                do while  (o3_in(j, k, l, m) /= missing_value .and. k < n_plev)
     275                   k = k + 1
     276                end do
     277                ! Replace missing values with the valid value at the
     278                ! lowest level above missing values:
     279                if (o3_in(j, k, l, m) == missing_value) &
     280                     o3_in(j, k:n_plev, l, m) = o3_in(j, k-1, l, m)
     281             end do
    206282          end do
    207        end do
    208     else
    209        print *, "regr_lat_time_climoz: no missing value attribute"
    210     end if
    211    
     283       else
     284          print *, "regr_lat_time_climoz: field ", m, &
     285               ", no missing value attribute"
     286       end if
     287    end do
     288
    212289    call nf95_close(ncid_in)
    213290
    214     allocate(o3_regr_lat(jjm + 1, n_plev, 0:13))
    215     allocate(o3_out(jjm + 1, n_plev, 360))
     291    allocate(o3_regr_lat(jjm + 1, n_plev, 0:13, read_climoz))
     292    allocate(o3_out(jjm + 1, n_plev, 360, read_climoz))
    216293
    217294    ! Regrid in latitude:
    218295    ! We average with respect to sine of latitude, which is
    219296    ! equivalent to weighting by cosine of latitude:
    220     if (size(o3_in, 3) == 12) then
    221        print *, "Found 12 months in ozone climatology, assuming periodicity..."
    222        o3_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_in, &
     297    if (n_month == 12) then
     298       print *, &
     299            "Found 12 months in ozone climatologies, assuming periodicity..."
     300       o3_regr_lat(jjm+1:1:-1, :, 1:12, :) = regr1_step_av(o3_in, &
    223301            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
    224302       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
    … …  
    227305       ! Duplicate January and December values, in preparation of time
    228306       ! interpolation:
    229        o3_regr_lat(:, :, 0) = o3_regr_lat(:, :, 12)
    230        o3_regr_lat(:, :, 13) = o3_regr_lat(:, :, 1)
     307       o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
     308       o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
    231309    else
    232        print *, "Using 14 months in ozone climatology..."
    233        o3_regr_lat(jjm+1:1:-1, :, :) = regr1_step_av(o3_in, &
     310       print *, "Using 14 months in ozone climatologies..."
     311       o3_regr_lat(jjm+1:1:-1, :, :, :) = regr1_step_av(o3_in, &
    234312            xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
    235313       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
    236314       ! in descending order)
    237315    end if
    238    
    239     deallocate(o3_in) ! pointer
    240316
    241317    ! Regrid in time by linear interpolation:
    … …  
    243319
    244320    ! Write to file:
    245     call nf95_put_var(ncid_out, varid_out, o3_out(jjm+1:1:-1, :, :))
    246     ! (The order of "rlatu" is inverted in the output file)
     321    do m = 1, read_climoz
     322       call nf95_put_var(ncid_out, varid_out(m), o3_out(jjm+1:1:-1, :, :, m))
     323       ! (The order of "rlatu" is inverted in the output file)
     324    end do
    247325
    248326    call nf95_close(ncid_out)
    … …  
    263341
    264342    integer, intent(in):: ncid_in, n_plev
    265     integer, intent(out):: ncid_out, varid_out, varid_plev, varid_time
     343    integer, intent(out):: ncid_out, varid_plev, varid_time
     344
     345    integer, intent(out):: varid_out(:) ! dim(1 or 2)
     346    ! "varid_out(1)" is for average ozone day and night,
     347    ! "varid_out(2)" is for daylight ozone.
    266348
    267349    ! Variables local to the procedure:
    … …  
    307389    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
    308390
    309     ! Define the primary variable:
     391    ! Define the primary variables:
    310392
    311393    call nf95_def_var(ncid_out, "tro3", nf90_float, &
    312          (/dimid_rlatu, dimid_plev, dimid_time/), varid_out)
    313     call nf95_put_att(ncid_out, varid_out, "long_name", "ozone mole fraction")
    314     call nf95_put_att(ncid_out, varid_out, "standard_name", &
     394         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(1))
     395    call nf95_put_att(ncid_out, varid_out(1), "long_name", &
     396         "ozone mole fraction")
     397    call nf95_put_att(ncid_out, varid_out(1), "standard_name", &
    315398         "mole_fraction_of_ozone_in_air")
     399
     400    if (size(varid_out) == 2) then
     401       call nf95_def_var(ncid_out, "tro3_daylight", nf90_float, &
     402            (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(2))
     403       call nf95_put_att(ncid_out, varid_out(2), "long_name", &
     404            "ozone mole fraction in daylight")
     405    end if
    316406
    317407    ! Global attributes:

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/regr_pr_int_m.F90

    r1239 r1263  
    11! $Id$
    2 module regr_pr
     2module regr_pr_int_m
    33
    44  ! Author: Lionel GUEZ
    5 
    6   ! In both procedures of this module:
    7   ! -- the root process reads a 2D latitude-pressure field from a
    8   !    NetCDF file, at a given day.
    9   ! -- the field is packed to the LMDZ horizontal "physics"
    10   !    grid and scattered to all threads of all processes;
    11   ! -- in all the threads of all the processes, the field is regridded in
    12   !    pressure to the LMDZ vertical grid.
    13   ! We assume that, in the input file, the field has 3 dimensions:
    14   ! latitude, pressure, julian day.
    15   ! We assume that latitudes are in ascending order in the input file.
    16 
    175
    186  implicit none
    … …  
    208contains
    219
    22   subroutine regr_pr_av(ncid, name, julien, press_in_edg, paprs, v3)
    23 
    24     ! "regr_pr_av" stands for "regrid pressure averaging".
    25     ! The target vertical LMDZ grid is the grid of layer boundaries.
    26     ! Regridding in pressure is done by averaging a step function of pressure.
    27 
    28     use dimphy, only: klon
    29     use netcdf95, only: nf95_inq_varid, handle_err
    30     use netcdf, only: nf90_get_var
    31     use assert_m, only: assert
    32     use regr1_step_av_m, only: regr1_step_av
    33     use mod_phys_lmdz_mpi_data, only: is_mpi_root
    34 
    35     use mod_phys_lmdz_transfert_para, only: scatter2d
    36     ! (pack to the LMDZ horizontal "physics" grid and scatter)
    37 
    38     integer, intent(in):: ncid ! NetCDF ID of the file
    39     character(len=*), intent(in):: name ! of the NetCDF variable
    40     integer, intent(in):: julien ! jour julien, 1 <= julien <= 360
    41 
    42     real, intent(in):: press_in_edg(:)
    43     ! edges of pressure intervals for input data, in Pa, in strictly
    44     ! ascending order
    45 
    46     real, intent(in):: paprs(:, :) ! (klon, llm + 1)
    47     ! (pression pour chaque inter-couche, en Pa)
    48 
    49     real, intent(out):: v3(:, :) ! (klon, llm)
    50     ! (regridded field on the partial "physics" grid)
    51     ! ("v3(i, k)" is at longitude "xlon(i)", latitude
    52     ! "xlat(i)", in pressure interval "[paprs(i, l+1), paprs(i, l)]".)
    53 
    54     ! Variables local to the procedure:
    55 
    56     include "dimensions.h"
    57     integer varid, ncerr ! for NetCDF
    58 
    59     real  v1(iim, jjm + 1, size(press_in_edg) - 1)
    60     ! (input field at day "julien", on the global "dynamics" horizontal grid)
    61     ! (First dimension is for longitude.
    62     ! The value is the same for all longitudes.
    63     ! "v1(:, j, k)" is at latitude "rlatu(j)" and for
    64     ! pressure interval "[press_in_edg(k), press_in_edg(k+1)]".)
    65 
    66     real v2(klon, size(press_in_edg) - 1)
    67     ! (field scattered to the partial "physics" horizontal grid)
    68     ! ("v2(i, k)" is at longitude "xlon(i)", latitude "xlat(i)" and
    69     ! for pressure interval "[press_in_edg(k), press_in_edg(k+1)]".)
    70 
    71     integer i
    72 
    73     !--------------------------------------------
    74 
    75     call assert(shape(v3) == (/klon, llm/), "regr_pr_av v3")
    76     call assert(shape(paprs) == (/klon, llm+1/), "regr_pr_av paprs")
    77 
    78     !$omp master
    79     if (is_mpi_root) then
    80        call nf95_inq_varid(ncid, name, varid)
    81 
    82        ! Get data at the right day from the input file:
    83        ncerr = nf90_get_var(ncid, varid, v1(1, :, :), start=(/1, 1, julien/))
    84        call handle_err("regr_pr_av nf90_get_var " // name, ncerr, ncid)
    85        ! Latitudes are in ascending order in the input file while
    86        ! "rlatu" is in descending order so we need to invert order:
    87        v1(1, :, :) = v1(1, jjm+1:1:-1, :)
    88 
    89        ! Duplicate on all longitudes:
    90        v1(2:, :, :) = spread(v1(1, :, :), dim=1, ncopies=iim-1)
    91     end if
    92     !$omp end master
    93 
    94     call scatter2d(v1, v2)
    95 
    96     ! Regrid in pressure at each horizontal position:
    97     do i = 1, klon
    98        v3(i, llm:1:-1) = regr1_step_av(v2(i, :), press_in_edg, &
    99             paprs(i, llm+1:1:-1))
    100        ! (invert order of indices because "paprs" is in descending order)
    101     end do
    102 
    103   end subroutine regr_pr_av
    104 
    105   !***************************************************************
    106 
    10710  subroutine regr_pr_int(ncid, name, julien, plev, pplay, top_value, v3)
    10811
    10912    ! "regr_pr_int" stands for "regrid pressure interpolation".
     13    ! In this procedure:
     14    ! -- the root process reads a 2D latitude-pressure field from a
     15    !    NetCDF file, at a given day.
     16    ! -- the field is packed to the LMDZ horizontal "physics"
     17    !    grid and scattered to all threads of all processes;
     18    ! -- in all the threads of all the processes, the field is regridded in
     19    !    pressure to the LMDZ vertical grid.
     20    ! We assume that, in the input file, the field has 3 dimensions:
     21    ! latitude, pressure, julian day.
     22    ! We assume that latitudes are in ascending order in the input file.
    11023    ! The target vertical LMDZ grid is the grid of mid-layers.
    11124    ! Regridding is by linear interpolation.
    … …  
    191104  end subroutine regr_pr_int
    192105
    193 end module regr_pr
     106end module regr_pr_int_m
Note: See TracChangeset for help on using the changeset viewer.