Ignore:
Timestamp:
Feb 2, 2017, 7:01:50 PM (7 years ago)
Author:
dcugnet
Message:

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.
Location:
LMDZ5/trunk/libf/phylmd
Files:
3 added
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/clesphys.h

    r2730 r2788  
    8787       LOGICAL :: ok_qch4
    8888       LOGICAL :: ok_conserv_q
     89       LOGICAL :: adjust_tropopause
     90       LOGICAL :: ok_daily_climoz
    8991
    9092       COMMON/clesphys/                                                 &
    … …  
    130132     &     , iflag_rrtm, ok_strato,ok_hines, ok_qch4                    &
    131133     &     , iflag_ice_thermo, ok_gwd_rando, NSW, iflag_albedo          &
    132      &     , ok_chlorophyll,ok_conserv_q, ok_all_xml
     134     &     , ok_chlorophyll,ok_conserv_q, adjust_tropopause             &
     135     &     , ok_daily_climoz, ok_all_xml
    133136     
    134137       save /clesphys/

  • LMDZ5/trunk/libf/phylmd/conf_phys_m.F90

    r2785 r2788  
    222222    LOGICAL,SAVE  :: carbon_cycle_tr_omp
    223223    LOGICAL,SAVE  :: carbon_cycle_cpl_omp
     224    LOGICAL,SAVE  :: adjust_tropopause_omp
     225    LOGICAL,SAVE  :: ok_daily_climoz_omp
    224226
    225227    INTEGER, INTENT(OUT):: read_climoz ! read ozone climatology, OpenMP shared
    … …  
    20352037    !Config Help = ...
    20362038    !
     2039    !
     2040    adjust_tropopause = .FALSE.
     2041    CALL getin('adjust_tropopause', adjust_tropopause_omp)
     2042    !
     2043    !Config Key  = adjust_tropopause
     2044    !Config Desc = Adjust the ozone field from the climoz file by stretching its
     2045    !              tropopause so that it matches the one of LMDZ.
     2046    !Config Def  = .FALSE.
     2047    !Config Help = Ensure tropospheric ozone column conservation.
     2048    !
     2049    !
     2050    ok_daily_climoz = .TRUE.
     2051    CALL getin('ok_daily_climoz', ok_daily_climoz_omp)
     2052    !
     2053    !Config Key  = ok_daily_climoz
     2054    !Config Desc = Interpolate in time the ozone forcings within ce0l.
     2055    !              .TRUE. if backward compatibility is needed.
     2056    !Config Def  = .TRUE.
     2057    !Config Help = .FALSE. ensure much fewer (no calendar dependency)
     2058    !  and lighter monthly climoz files, inetrpolated in time at gcm run time.
    20372059    !
    20382060    ecrit_LES_omp = 1./8.
    … …  
    22912313    callstats = callstats_omp
    22922314    ecrit_LES = ecrit_LES_omp
     2315    adjust_tropopause = adjust_tropopause_omp
     2316    ok_daily_climoz = ok_daily_climoz_omp
    22932317    carbon_cycle_tr = carbon_cycle_tr_omp
    22942318    carbon_cycle_cpl = carbon_cycle_cpl_omp
    … …  
    24852509    write(lunout,*)' aerosol_couple = ', aerosol_couple
    24862510    write(lunout,*)' flag_aerosol = ', flag_aerosol
    2487     write(lunout,*)' flag_aerosol_strat = ', flag_aerosol_strat
     2511    write(lunout,*)' flag_aerosol_strat= ', flag_aerosol_strat
    24882512    write(lunout,*)' new_aod = ', new_aod
    24892513    write(lunout,*)' aer_type = ',aer_type
    … …  
    25682592    write(lunout,*) 'SSO gkwake=',gkwake
    25692593    write(lunout,*) 'SSO gklift=',gklift
     2594    write(lunout,*) 'adjust_tropopause = ', adjust_tropopause
     2595    write(lunout,*) 'ok_daily_climoz = ',ok_daily_climoz
    25702596    write(lunout,*) 'read_climoz = ', read_climoz
    25712597    write(lunout,*) 'carbon_cycle_tr = ', carbon_cycle_tr

  • LMDZ5/trunk/libf/phylmd/limit_read_mod.F90

    r2770 r2788  
    223223          ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
    224224          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
    225           WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
     225          WRITE(abort_message,'(a,2(i3,a))')'limit.nc records number (',nn,') does no'//&
    226226            't match year length (',year_len,')'
    227227          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)

  • LMDZ5/trunk/libf/phylmd/open_climoz_m.F90

    r1907 r2788  
    66contains
    77
    8   subroutine open_climoz(ncid, press_in_edg)
     8  subroutine open_climoz(ncid, press_in_cen, press_in_edg, time_in)
    99
    1010    ! This procedure should be called once per "gcm" run, by a single
    1111    ! thread of each MPI process.
    1212    ! The root MPI process opens "climoz_LMDZ.nc", reads the pressure
    13     ! levels and broadcasts them to the other processes.
     13    ! levels and the times and broadcasts them to the other processes.
    1414
    1515    ! We assume that, in "climoz_LMDZ.nc", the pressure levels are in hPa
    … …  
    2121    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    2222    use mod_phys_lmdz_mpi_transfert, only: bcast_mpi ! broadcast
     23    USE phys_cal_mod, ONLY: calend, year_len, year_cur
    2324
     25    ! OpenMP shares arguments:
    2426    integer, intent(out):: ncid ! of "climoz_LMDZ.nc", OpenMP shared
    2527
    26     real, pointer:: press_in_edg(:)
    27     ! edges of pressure intervals for ozone climatology, in Pa, in strictly
    28     ! ascending order, OpenMP shared
     28    ! pressure levels for ozone climatology, in Pa, strictly ascending order
     29    real, pointer:: press_in_cen(:) !--- at cells centers
     30    real, pointer:: press_in_edg(:) !--- at the interfaces (pressure intervals)
     31
     32    real, pointer:: time_in(:)
     33    ! times of ozone climatology records, in days since Jan. 1st 0h
    2934
    3035    ! Variables local to the procedure:
    3136
    32     real, pointer:: plev(:)
    33     ! (pressure levels for ozone climatology, converted to Pa, in strictly
    34     ! ascending order)
    35 
    3637    integer varid ! for NetCDF
    3738    integer n_plev ! number of pressure levels in the input data
     39    integer n_time ! number of time records    in the input field
    3840    integer k
     41    CHARACTER(LEN=80)  :: modname
     42    CHARACTER(LEN=320) :: msg
    3943
    4044    !---------------------------------------
    4145
    42     print *, "Call sequence information: open_climoz"
     46    modname="open_climoz"
     47    print *, "Call sequence information: "//TRIM(modname)
    4348
    4449    if (is_mpi_root) then
    … …  
    4651
    4752       call nf95_inq_varid(ncid, "plev", varid)
    48        call nf95_gw_var(ncid, varid, plev)
     53       call nf95_gw_var(ncid, varid, press_in_cen)
    4954       ! Convert from hPa to Pa because "paprs" and "pplay" are in Pa:
    50        plev = plev * 100.
    51        n_plev = size(plev)
     55       press_in_cen = press_in_cen * 100.
     56       n_plev = size(press_in_cen)
     57       CALL NF95_INQ_VARID(ncID, "time", varID)
     58       CALL NF95_GW_VAR(ncid, varid, time_in)
     59       n_time = SIZE(time_in)
     60       IF(ALL([year_len,14]/=n_time)) THEN             !--- Check records number
     61         WRITE(msg,'(a,3(i4,a))')TRIM(modname)//': expected input file could be&
     62         & monthly or daily (14 or ',year_len,' records for year ',year_cur,' a&
     63         &nd calendar "'//TRIM(calend)//'"), but ',n_time,' records were found'
     64         CALL abort_physic(modname, msg, 1)
     65       END IF
     66
    5267    end if
    5368
    5469    call bcast_mpi(n_plev)
    55     if (.not. is_mpi_root) allocate(plev(n_plev))
    56     call bcast_mpi(plev)
    57    
     70    if (.not. is_mpi_root) allocate(press_in_cen(n_plev))
     71    call bcast_mpi(press_in_cen)
     72    call bcast_mpi(n_time)
     73    if (.not. is_mpi_root) allocate(time_in(n_time))
     74    call bcast_mpi(time_in)
     75
    5876    ! Compute edges of pressure intervals:
    5977    allocate(press_in_edg(n_plev + 1))
    … …  
    6179       press_in_edg(1) = 0.
    6280       ! We choose edges halfway in logarithm:
    63        forall (k = 2:n_plev) press_in_edg(k) = sqrt(plev(k - 1) * plev(k))
     81       forall (k = 2:n_plev) press_in_edg(k) = sqrt(press_in_cen(k - 1) * press_in_cen(k))
    6482       press_in_edg(n_plev + 1) = huge(0.)
    6583       ! (infinity, but any value guaranteed to be greater than the
    … …  
    6785    end if
    6886    call bcast_mpi(press_in_edg)
    69     deallocate(plev) ! pointer
     87!    deallocate(press_in_cen) ! pointer
    7088
    7189  end subroutine open_climoz

  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r2752 r2788  
    487487!$OMP THREADPRIVATE(budg_sed_part)
    488488#endif
     489      REAL, ALLOCATABLE, SAVE, DIMENSION(:) :: pr_tropopause
     490!$OMP THREADPRIVATE(pr_tropopause)
    489491
    490492CONTAINS
    … …  
    764766      ALLOCATE (vsed_aer(klon,klev))
    765767#endif
     768      ALLOCATE (pr_tropopause(klon))
    766769
    767770END SUBROUTINE phys_local_var_init
    … …  
    10181021      DEALLOCATE (budg_sed_part)
    10191022#endif
     1023      DEALLOCATE (pr_tropopause)
    10201024
    10211025END SUBROUTINE phys_local_var_end

  • LMDZ5/trunk/libf/phylmd/physiq_mod.F90

    r2784 r2788  
    3535    USE change_srf_frac_mod
    3636    USE surface_data,     ONLY : type_ocean, ok_veget, ok_snow
     37    USE tropopause_m,     ONLY: dyn_tropopause
    3738#ifdef CPP_Dust
    3839    USE phytracr_spl_mod, ONLY: phytracr_spl
    … …  
    191192       beta_prec,  &
    192193       rneb,  &
    193        zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
     194       zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic, &
     195       pr_tropopause
    194196       !
    195197    USE phys_state_var_mod ! Variables sauvegardees de la physique
    … …  
    205207    USE phys_output_ctrlout_mod
    206208    use open_climoz_m, only: open_climoz ! ozone climatology from a file
    207     use regr_pr_av_m, only: regr_pr_av
     209    use regr_pr_time_av_m, only: regr_pr_time_av
    208210    use netcdf95, only: nf95_close
    209211    !IM for NMC files
    … …  
    10061008    !$OMP THREADPRIVATE(first)
    10071009
    1008     integer, save::  read_climoz ! read ozone climatology
     1010    ! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
     1011    ! Note that pressure vectors are in Pa and in stricly ascending order
     1012    INTEGER,SAVE :: read_climoz                ! Read ozone climatology
    10091013    !     (let it keep the default OpenMP shared attribute)
    10101014    !     Allowed values are 0, 1 and 2
    … …  
    10131017    !     2: read two ozone climatologies, the average day and night
    10141018    !     climatology and the daylight climatology
    1015 
    1016     integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
    1017     !     (let it keep the default OpenMP shared attribute)
    1018 
    1019     real, pointer, save:: press_climoz(:)
    1020     ! (let it keep the default OpenMP shared attribute)
    1021     ! edges of pressure intervals for ozone climatologies, in Pa, in strictly
    1022     ! ascending order
    1023 
    1024     integer ro3i
    1025     !     required time index in NetCDF file for the ozone fields, between 1
    1026     !     and 360
     1019    INTEGER,SAVE :: ncid_climoz                ! NetCDF file identifier
     1020    REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
     1021    REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
     1022    REAL, POINTER, SAVE :: time_climoz(:)      ! Time vector
     1023    CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
     1024                                  = ["tro3         ","tro3_daylight"]
     1025    ! vars_climoz(1:read_climoz): variables names in climoz file.
    10271026
    10281027    INTEGER ierr
    … …  
    16711670
    16721671       !$omp single
    1673        IF (read_climoz >= 1) THEN
    1674           CALL open_climoz(ncid_climoz, press_climoz)
    1675        ENDIF
     1672       IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz,   &
     1673                                         press_edg_climoz, time_climoz)
    16761674       !$omp end single
    16771675       !
    … …  
    19781976          wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
    19791977       ELSE
    1980           ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
    1981 
    1982           IF (ro3i == 361) ro3i = 360
    1983           ! (This should never occur, except perhaps because of roundup
    1984           ! error. See documentation.)
    1985 
    1986           IF (read_climoz == 1) THEN
    1987              CALL regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
    1988                   press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1978          IF(adjust_tropopause) THEN
     1979            WRITE(*,*)' >> Adjusting O3 field to match gcm tropopause.'
     1980            CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot,  &
     1981                 pr_tropopause)
     1982            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
     1983                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,   &
     1984                 time_climoz,  latitude_deg,  press_cen_climoz, pr_tropopause)
    19891985          ELSE
    1990              ! read_climoz == 2
    1991              CALL regr_pr_av(ncid_climoz, (/"tro3         ", &
    1992                   "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
    1993                   paprs=paprs, v3=wo)
    1994           ENDIF
     1986            WRITE(*,*)' >> Interpolating O3 field directly on gcm levels.'
     1987            CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz),    &
     1988                 days_elapsed+jh_cur-jh_1jan, press_edg_climoz, paprs, wo,  &
     1989                 time_climoz)
     1990          END IF
    19951991          ! Convert from mole fraction of ozone to column density of ozone in a
    19961992          ! cell, in kDU:
    19971993          FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
    19981994               * zmasse / dobson_u / 1e3
    1999           ! (By regridding ozone values for LMDZ only once per day, we
     1995          ! (By regridding ozone values for LMDZ only once a day, we
    20001996          ! have already neglected the variation of pressure in one
    20011997          ! day. So do not recompute "wo" at each time step even if
    … …  
    44654461             CALL nf95_close(ncid_climoz)
    44664462          ENDIF
    4467           DEALLOCATE(press_climoz) ! pointer
     4463          DEALLOCATE(press_edg_climoz) ! pointer
     4464          DEALLOCATE(press_cen_climoz) ! pointer
    44684465       ENDIF
    44694466       !$OMP END MASTER

  • LMDZ5/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90

    r2440 r2788  
    4141
    4242    use mod_grid_phy_lmdz, ONLY : nbp_lat
    43     use regr1_conserv_m, only: regr1_conserv
    44     use regr3_lint_m, only: regr3_lint
     43    use regr_conserv_m, only: regr_conserv
     44    use regr_lint_m, only: regr_lint
    4545    use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, handle_err, &
    4646         nf95_put_var, nf95_gw_var
    … …  
    8484    ! Last dimension is month number.)
    8585
    86     real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, 360)
     86    real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, ndays)
    8787    ! (regridded ozone parameter)
    8888    ! ("o3_par_out(j, l, day)" is at latitude "rlatu(j)", pressure
    … …  
    162162    latitude = latitude / 180. * pi
    163163    n_lat = size(latitude)
    164     ! We need to supply the latitudes to "regr1_conserv" in
     164    ! We need to supply the latitudes to "regr_conserv" in
    165165    ! ascending order, so invert order if necessary:
    166166    desc_lat = latitude(1) > latitude(n_lat)
    … …  
    209209       ! We average with respect to sine of latitude, which is
    210210       ! equivalent to weighting by cosine of latitude:
    211        call regr1_conserv(o3_par_in, xs = sin(lat_in_edg), &
     211       call regr_conserv(1, o3_par_in, xs = sin(lat_in_edg), &
    212212            xt = (/-1., sin((/boundslat_reg(nbp_lat-1:1:-1,south)/)), 1./), &
    213213            vt = v_regr_lat(nbp_lat:1:-1, :, 1:12))
    … …  
    221221
    222222       ! Regrid in time by linear interpolation:
    223        o3_par_out = regr3_lint(v_regr_lat, tmidmonth, tmidday)
     223       call regr_lint(3, v_regr_lat, tmidmonth, tmidday, o3_par_out)
    224224
    225225       ! Write to file:

  • LMDZ5/trunk/libf/phylmd/regr_pr_comb_coefoz_m.F90

    r2346 r2788  
    7676    use dimphy, only: klon
    7777    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    78     use regr_pr_av_m, only: regr_pr_av
     78    use regr_pr_time_av_m, only: regr_pr_time_av
    7979    use regr_pr_int_m, only: regr_pr_int
    8080    use press_coefoz_m, only: press_in_edg, plev
    … …  
    128128    !$omp end master
    129129
    130     call regr_pr_av(ncid, (/"a2       ", "a4       ", "a6       ", &
    131          "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), julien, &
     130    call regr_pr_time_av(ncid, (/"a2       ", "a4       ", "a6       ", &
     131         "P_net_Mob", "r_Mob    ", "temp_Mob ", "R_Het    "/), REAL(julien-1), &
    132132         press_in_edg, paprs, coefoz)
    133133    a2 = coefoz(:, :, 1)

  • LMDZ5/trunk/libf/phylmd/regr_pr_int_m.F90

    r2346 r2788  
    2828    use netcdf, only: nf90_get_var
    2929    use assert_m, only: assert
    30     use regr1_lint_m, only: regr1_lint
     30    use regr_lint_m, only: regr_lint
    3131    use mod_phys_lmdz_mpi_data, only: is_mpi_root
    3232    use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat, nbp_lev
    … …  
    9797    ! Regrid in pressure at each horizontal position:
    9898    do i = 1, klon
    99        v3(i, nbp_lev:1:-1) = regr1_lint(v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1))
     99       call regr_lint(1, v2(i, :), (/0., plev/), pplay(i, nbp_lev:1:-1), &
     100                         v3(i, nbp_lev:1:-1))
    100101       ! (invert order of indices because "pplay" is in descending order)
    101102    end do

  • LMDZ5/trunk/libf/phylmd/regr_pr_o3_m.F90

    r2440 r2788  
    2828    use netcdf, only:  nf90_nowrite, nf90_get_var
    2929    use assert_m, only: assert
    30     use regr1_conserv_m, only: regr1_conserv
     30    use regr_conserv_m, only: regr_conserv
    3131    use press_coefoz_m, only: press_in_edg
    3232    use time_phylmdz_mod, only: day_ref
    … …  
    7575    ! Poles:
    7676    do j = 1, nbp_lat, nbp_lat-1
    77        call regr1_conserv(r_mob(j, :), press_in_edg, &
     77       call regr_conserv(1, r_mob(j, :), press_in_edg, &
    7878            p3d(1, j, nbp_lev + 1:1:-1), o3_mob_regr(1, j, nbp_lev:1:-1))
    7979       ! (invert order of indices because "p3d" is in descending order)
    … …  
    8383    do j = 2, nbp_lat-1
    8484       do i = 1, nbp_lon
    85           call regr1_conserv(r_mob(j, :), press_in_edg, &
     85          call regr_conserv(1, r_mob(j, :), press_in_edg, &
    8686               p3d(i, j, nbp_lev + 1:1:-1), o3_mob_regr(i, j, nbp_lev:1:-1))
    8787          ! (invert order of indices because "p3d" is in descending order)
Note: See TracChangeset for help on using the changeset viewer.