ecearth3/sources/tm5mp/proj/cb05/boundary.F90

#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if

#include "tm5.inc"

!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:     BOUNDARY
!
! !DESCRIPTION: Set up boundary conditions:
!
!      O3, O3S and CH4 are relaxed to climatology at stratospheric layers
!       
!      For O3 we use either
!      CMIP6 ozone 3D mixing ratios 1850-2100 (Hegglin et al., in prep) or
!      Multi Sensor Reanalysis v2
!      Source MSR : www.temis.nl; Allaart, van der A, manuscript in preparation, 2009)
!      Source MSR2: FIXME
!       
!      For CH4 we use Climatology HALOE CH4, October 1991 to August 2002 by:
!      Technical note: A stratospheric climatology for O3, H2O, CH4, NOx, HCl
!        and HF derived from HALOE measurements by J.-U. Grooss and J. M.
!        Russell III, Atmospheric Chemistry and Physics, 5, 2797-2807, 2005
!        SRef-ID: 1680-7324/acp/2005-5-2797
!       
!      HNO3 is prescribed at level lm using HNO3:03 ratios from ODIN (fallback on UARS). 
!
!      CO is nudged using the ODIN ratio of CO/O3 @ 3 levels (JEW, 2014).
!
!      Note about ODIN datasets: a slow relaxation between monthly mean values
!      is incorporated to improve on the large variability which exists between
!      consecutive months (JEW, 2014).
!\\
!\\
! !INTERFACE: 
!
MODULE BOUNDARY
  !
  ! !USES:
  !
  USE GO,            ONLY : gol, goPr, goErr, goLabel
  USE TM5_DISTGRID,  ONLY : dgrid, Get_DistGrid, gather, scatter_i_band
  USE DIMS,          ONLY : nregions, idate
  USE chem_param,    ONLY : ntracet
  USE global_types,  ONLY : d23_data, d3_data
  USE Grid,          ONLY : TLevelInfo

#ifdef with_budgets
  USE budget_global, ONLY : nbudg, nbud_vg, nzon_vg, budg_dat
#endif

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: BOUNDARY_INIT, BOUNDARY_DONE, BOUNDARY_APPLY
  !
  ! !PUBLIC DATA MEMBERS:
  !
  LOGICAL,        PUBLIC :: use_o3du
  TYPE(d23_data), PUBLIC :: o3du(nregions) ! optionally used in photolysis only
  LOGICAL,        PUBLIC :: LCMIP6_CO2
  REAL,           PUBLIC :: co2_glob
  !
  ! !PRIVATE DATA MEMBERS:
  !
  ! Flags for O3, CO and CH4 stratospheric nudging
  LOGICAL            :: LCMIP6_O3, LCMIP6_CH4
  LOGICAL            :: o3_piclim
  REAL, ALLOCATABLE  :: weight_cmip62tm5(:)
  INTEGER, ALLOCATABLE :: jlow_cmip62tm5(:)
  LOGICAL            :: use_MSR
  LOGICAL            :: use_HALOE
  LOGICAL            :: use_ODIN   
  LOGICAL            :: emis_ch4_single, emis_ch4_fix3d 
  LOGICAL            :: o3_fixyear, ch4_fixyear, co2_fixyear
  INTEGER            :: o3_year, ch4_year, co2_year

  character(len=256) :: o3vmr_dirname
  character(len=256) :: o3du_dirname
  character(len=256) :: ch4vmr_dirname
  character(len=256) :: covmr_dirname
  character(len=256) :: cmip6_ch4_dirname_strat  
  character(len=256) :: cmip6_co2_dirname

  ! Volume Mixing Ratio and Ratio, and values in DU
  TYPE(d23_data)     :: o3vmrpm(nregions)
  TYPE(d23_data)     :: o3vmr(nregions)
  TYPE(d23_data)     :: o3vmrnm(nregions)
  TYPE(d23_data)     :: ch4vmrpm(nregions)
  TYPE(d23_data)     :: ch4vmr(nregions)
  TYPE(d23_data)     :: ch4vmrnm(nregions)
  TYPE(d23_data)     :: o3ratpm(nregions)
  TYPE(d23_data)     :: o3rat(nregions)
  TYPE(d23_data)     :: o3ratnm(nregions)
  TYPE(d23_data)     :: ch4ratpm(nregions)
  TYPE(d23_data)     :: ch4rat(nregions)
  TYPE(d23_data)     :: ch4ratnm(nregions)

  TYPE(d23_data)     :: odin_hno3_o3(nregions)    ! hno3/o3 ratio at 4 pressure levels
  TYPE(d23_data)     :: odin_hno3_o3_pm(nregions) ! hno3/o3 ratio at 4 pressure levels (previous month)
  TYPE(d23_data)     :: odin_hno3_o3_nm(nregions) ! hno3/o3 ratio at 4 pressure levels (next month)
  TYPE(d23_data)     :: odin_co_o3(nregions)      ! co/o3 ratio at 4 pressure levels
  TYPE(d23_data)     :: odin_co_o3_pm(nregions)   ! co/o3 ratio at 4 pressure levels (previous month)
  TYPE(d23_data)     :: odin_co_o3_nm(nregions)   ! co/o3 ratio at 4 pressure levels (next month)

  TYPE(TLevelInfo)   :: LeviX_msr, LeviX_haloe, LeviX_cmip6_o3
    
  real, allocatable  :: wrk2d_ML(:,:) ! 2D work array, model levels
  real, allocatable  :: wrk2d_4L(:,:) ! 2D work array,     4 levels
    
  real, allocatable  :: wrk3dpm_ML(:,:,:), wrk3d_ML(:,:,:), wrk3dnm_ML(:,:,:)! 3D work array, model levels
  real, allocatable  :: wrk3dratio(:,:,:)
  real, allocatable  :: wrk3dratiopm(:,:,:)  ! for smoothing CO and HNO3 values at monthly scale  
  real, allocatable  :: wrk3drationm(:,:,:)  ! for smoothing CO and HNO3 values at monthly scale

#ifdef with_budgets
!  REAL, DIMENSION(nregions)                ::  sum_stratosphere
  REAL,     PUBLIC   ::  budstratg(nbudg, nbud_vg, ntracet)
#endif

  integer ::  itim_init, itim_appl ! Timers id

  CHARACTER(len=*), PARAMETER  ::  mname = 'boundary'
  !
  ! !REVISION HISTORY: 
  !   15 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !   03 Mar 2014 - J. E. Williams - added nudging for CO from ODIN measurements
  !   Dec 2016 - Mar 2017 - T. van Noije - added options for using CMIP6 data sets 
  !                                        for O3, CH4, and CO2
  !
  ! !REMARKS:
  !   (1) All reference to sum_stratosphere have been commented, since the
  !       variable is neither printed nor saved nor used to compute something
  !       else.
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BOUNDARY_INIT
  !
  ! !DESCRIPTION: Act as both INIT and MONTHLY_UPDATE methods.
  ! 
  !               As INIT           : read settings from rc file & allocate data.
  !               As MONTHLY_UPDATE : read data in.
  !               
  !               Called at run start and at start of every month (from
  !                sources_sinks/sources_sinks_init & ss_monthly_update).
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE BOUNDARY_INIT( first, status )
    !
    ! !USES:
    !
    USE GO,            ONLY : TrcFile, Init, Done, ReadRc
    USE GO,            ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
    USE dims,          ONLY : im, jm, lm, lme, okdebug, idate, iglbsfc
    USE dims,          ONLY : nregions, nlat180, nlon360
    USE dims,          ONLY : newyr
    USE global_data,   ONLY : rcfile, inputdir
    USE partools,      ONLY : isRoot, par_broadcast
    USE MDF,           ONLY : MDF_Open, MDF_HDF4, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close

    USE binas,         ONLY : p0
    USE Grid,          ONLY : FillGrid, Fill3D, FillLevels
    USE Grid,          ONLY : TLevelInfo
    USE Grid,          ONLY : TllGridInfo, Init
    USE MeteoData,     ONLY : global_lli, lli_z, levi
    !
    ! !INPUT PARAMETERS:
    !
    LOGICAL, INTENT(in)    :: first   ! is a new run
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)   :: status
    !
    ! !REVISION HISTORY: 
    !   13 Dec 2010 - P. Le Sager - bug fix in reading o3du : use correct
    !                  number of levels to read data, and allow vertical
    !                  regridding of O3du.
    !   23 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !                             - fixed : now run longer than a month have correct boundary
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*),      PARAMETER ::  rname = mname//'/Boundary_Init'
    INTEGER,               PARAMETER ::  avg_field = 1
    INTEGER, DIMENSION(2), PARAMETER :: msr2_valid = (/1979, 2015/) ! 2013-15 are GOME2 data

    INTEGER,               PARAMETER :: nlon_cmip6_o3 = 144
    INTEGER,               PARAMETER :: nlat_cmip6_o3 =  96
    INTEGER,               PARAMETER :: nlev_cmip6_o3 =  66
    REAL, ALLOCATABLE                :: a_cmip6_o3(:),b_cmip6_o3(:)
    REAL, ALLOCATABLE                :: lat_cmip6_o3(:)

    CHARACTER(len=9)     :: rgtype

    TYPE(TrcFile)        ::  rcF
    character(len=256)   ::  o3vmr_fnamepm, o3vmr_fname, o3vmr_fnamenm, o3du_fname
    character(len=256)   ::  ch4vmr_fnamepm, ch4vmr_fname, ch4vmr_fnamenm
    CHARACTER(len=256)   ::  filemon
    CHARACTER(len=4)     ::  ecstring

    INTEGER              ::  nmw, i1, i2, j2, j1, hid, varid
    INTEGER              ::  i,j,lmr, inp_levs, k
    INTEGER              ::  valid_year_pm, valid_year, valid_year_nm
    REAL*4               ::  field3d_r4(nlat180,4,12)
    INTEGER              ::  nm, pm
    TYPE(TLevelInfo)     ::  src_lev

    ! some arrays on 1x1 resolution
    REAL, ALLOCATABLE    ::  field3d_pm(:,:,:), field3d(:,:,:), field3d_nm(:,:,:), field3d_h(:,:,:)
    REAL, ALLOCATABLE    ::  sp_z(:,:)

    ! CMIP6 ozone fields
    REAL, ALLOCATABLE    ::  field4d_cmip6_h(:,:,:,:)
    REAL, ALLOCATABLE    ::  field3d_cmip6_pm(:,:), field3d_cmip6(:,:), field3d_cmip6_nm(:,:)
    REAL                 ::  tm5_lat
    INTEGER              ::  j_cmip6, jlow_cmip6
    INTEGER              ::  l ! TESTING ONLY

    ! Scale factors for stratospheric CH4 from HALOE
    real                 ::  ch4_scale, ch4_scale_pm, ch4_scale_nm
    real                 ::  ch4_ref
    integer              ::  iyear, target_year
    ! CMIP6 global annual mean mixing ratios for CH4 and CO2
    real*4, allocatable  ::  ch4_hist(:,:), ch4_sce(:,:), co2_hist(:,:), co2_sce(:,:)
    character(len=512)   ::  ch4_hist_fname, ch4_sce_fname, co2_hist_fname, co2_sce_fname
    !
    !  Scaling on HALOE climatology from 2000 
    !
    REAL,Dimension(37)   ::  RCP6_CH4_STRAT ! covers 1999-2035

    data RCP6_CH4_STRAT /1.00, 1.00, 1.0043444078, 1.0099155518, 1.0159814566, 1.021545823,  1.0269678609, 1.032430564, &
                       1.037717051,  1.0432000868, 1.0489868922, 1.0547802041, 1.0604597888, 1.0660828758, &
                       1.071659984,  1.0772118797, 1.0827461809, 1.0882515012, 1.0937098125, 1.0991152047, &
                       1.1044713649, 1.1097834167, 1.1150356906, 1.1202446152, 1.1254644653, 1.1307196942, &
                       1.1360154799, 1.1413492741, 1.146701693,  1.1520482561, 1.1573782007, 1.1626917436, &
                       1.1679937104, 1.1733492199, 1.1788363765, 1.1844880647, 1.1903212557/

    ! --- begin ----------------------------------------

    CALL goLabel(rname)

    !--------------------------------------------------
    ! ** TRUE INIT : only once
    !--------------------------------------------------    
    IF ( FIRST ) THEN

       ! read settings from rcfile:
       CALL Init( rcF, rcfile, status )
       IF_NOTOK_RETURN(status=1)

       CALL ReadRc( rcF, 'input.conc.o3.cmip6', LCMIP6_O3, status, default=.FALSE. )
       IF_ERROR_RETURN(status=1)

       IF (LCMIP6_O3) THEN
          write(gol,*) 'Stratospheric O3 based on CMIP6 mixing ratios'; call goPr

          CALL ReadRc( rcF, 'input.conc.o3.cmip6.dir', o3vmr_dirname, status )
          IF_NOTOK_RETURN(status=1)

          CALL ReadRc( rcF, 'input.conc.o3.cmip6.piclim', o3_piclim, status, default=.FALSE. )
          IF_ERROR_RETURN(status=1)

          IF (o3_piclim) THEN
             write(gol,*) 'Stratospheric O3 nudged to pre-industrial climatology from CMIP6'; call goPr
          ENDIF
       ENDIF
       
       CALL ReadRc( rcF, 'input.climat.MSR', use_MSR, status, default=.FALSE. )
       IF_ERROR_RETURN(status=1)

       IF (use_MSR) THEN
          CALL ReadRc( rcF, 'input.climat.o3vmr', o3vmr_dirname, status )
          IF_NOTOK_RETURN(status=1)

          CALL ReadRc( rcF, 'input.climat.use_o3du', use_o3du, status, default=.FALSE. )
          IF_ERROR_RETURN(status=1)

          IF (use_o3du) THEN
             CALL ReadRc( rcF, 'input.climat.o3du',  o3du_dirname, status  )
             IF_NOTOK_RETURN(status=1)
          ENDIF
       ENDIF

       IF (LCMIP6_O3 .and. use_MSR) THEN
          write (gol,'("ERROR - CMIP6 and MSR O3 cannot both be selected")') ; call goErr
          status=1; TRACEBACK; return 
       ENDIF

       ! It is now possible to nudge stratospheric O3 mixing ratios
       ! to the fields for a fixed year,
       ! both for the CMIP6 and MSR data sets.
       CALL ReadRc( rcF, 'input.o3.fixyear', o3_fixyear, status, default = .FALSE. )
       IF_ERROR_RETURN(status=1)

       IF (o3_fixyear) THEN
          CALL ReadRc( rcF, 'input.o3.year', o3_year, status )
          IF_NOTOK_RETURN(status=1)
          
          IF (LCMIP6_O3) THEN
             IF (.not.o3_piclim) THEN
                write(gol,*) 'Stratospheric O3 nudged to fixed year: ', o3_year; call goPr
             ENDIF
          ELSE IF (use_MSR) THEN
             ! o3_piclim not defined
             write(gol,*) 'Stratospheric O3 nudged to fixed year: ', o3_year; call goPr
          ENDIF
       ENDIF

       CALL ReadRc( rcF, 'input.climat.HALOE', use_HALOE, status, default=.FALSE. )
       IF_ERROR_RETURN(status=1)

       IF (use_HALOE) THEN
          CALL ReadRc( rcF, 'input.climat.ch4vmr', ch4vmr_dirname, status )
          IF_NOTOK_RETURN(status=1)
       ENDIF

       CALL ReadRc( rcF, 'input.conc.ch4.cmip6', LCMIP6_CH4, status, default=.TRUE. )
       IF_ERROR_RETURN(status=1)

       IF (LCMIP6_CH4) THEN

          write(gol,*) 'Stratospheric CH4 scaled based on CMIP6 surface concentration'; call goPr

          CALL ReadRc( rcF, 'input.conc.ch4.cmip6.dir.year', cmip6_ch4_dirname_strat, status )
          IF_NOTOK_RETURN(status=1)

          IF (.not.use_HALOE) THEN
             write (gol,'("ERROR - The HALOE stratospheric CH4 climatology should be used with CMIP6 CH4")') ; call goErr
             status=1; TRACEBACK; return
          ENDIF

       ELSE

          CALL ReadRc( rcF, 'input.emis.ch4.single', emis_ch4_single, status )
          IF_NOTOK_RETURN(status=1)
 
          IF (emis_ch4_single) THEN          
             call ReadRc( rcF, 'input.emis.ch4.fix3d', emis_ch4_fix3d, status, default=.true. )
             IF_NOTOK_RETURN(status=1)

             IF ( use_HALOE .and. emis_ch4_fix3d ) THEN
               write (gol,'("ERROR - Do not use the HALOE stratospheric CH4 climatology")') ; call goErr
               write (gol,'("ERROR - when a single mixing ratio is prescribed in the whole atmosphere")') ; call goErr
               status=1; TRACEBACK; return
             ENDIF
          ENDIF

       ENDIF

       CALL ReadRc( rcF, 'input.ch4.fixyear', ch4_fixyear, status, default = .FALSE. )
       IF_ERROR_RETURN(status=1)

       IF (ch4_fixyear) THEN
          CALL ReadRc( rcF, 'input.ch4.year', ch4_year, status )
          IF_NOTOK_RETURN(status=1)

          write(gol,*) 'Base year for stratospheric CH4 scaling: ', ch4_year; call goPr 
       ENDIF 

       CALL ReadRc( rcF, 'input.climat.ODIN', use_ODIN, status, default=.FALSE. )
       IF_ERROR_RETURN(status=1)       

       IF (use_ODIN) THEN
          CALL ReadRc( rcF, 'input.climat.covmr', covmr_dirname, status )
          IF_NOTOK_RETURN(status=1)
       ENDIF

       CALL ReadRc( rcF, 'input.conc.co2.cmip6', LCMIP6_CO2, status, default=.TRUE. )
       IF_ERROR_RETURN(status=1)

       IF (LCMIP6_CO2) THEN
          write(gol,*) 'Global annual mean CO2 mixing ratios in pH calculation based on CMIP6'; call goPr

          CALL ReadRc( rcF, 'input.conc.co2.cmip6.dir', cmip6_co2_dirname, status )
          IF_NOTOK_RETURN(status=1)

          CALL ReadRc( rcF, 'input.co2.fixyear', co2_fixyear, status, default = .FALSE. )
          IF_ERROR_RETURN(status=1)

          IF (co2_fixyear) THEN
             CALL ReadRc( rcF, 'input.co2.year', co2_year, status )
             IF_NOTOK_RETURN(status=1)

             write(gol,*) 'CO2 mixing ratio fixed to year: ', co2_year; call goPr
          ENDIF 
       ELSE    
          ! It would be better to change this warning into an error
          write (gol,'("WARNING - CO2 mixing ratio in pH calculation fixed to global annual mean for year 2000")') ; call goPr
          write (gol,'("WARNING - It is recommended to use the CMIP6 time series")') ; call goPr
       ENDIF

       CALL Done( rcF, status )
       IF_NOTOK_RETURN(status=1)

       IF (isRoot) THEN

          IF (LCMIP6_O3) THEN

             allocate(a_cmip6_o3(0:nlev_cmip6_o3))
             allocate(b_cmip6_o3(0:nlev_cmip6_o3))
             allocate(lat_cmip6_o3(1:nlat_cmip6_o3))
             allocate(jlow_cmip62tm5(1:jm(1)))
             allocate(weight_cmip62tm5(1:jm(1)))

             ! CMIP6 fields are provided on a pressure grid spanning from 1000 hPa to 0.0001 hPa
             ! half-level pressure coefficient a (Pa)
             a_cmip6_o3(:) = 100. * &
                  (/1037.5,   962.5,   887.5,   825.0,   790.0,  765.0,  725.0,  675.0,  625.0,  550.0, &
                     475.0,   425.0,   375.0,   325.0,   292.5,  267.5,  225.0,  185.0,  160.0,  140.0, &
                     122.5,   107.5,    95.0,    85.0,    75.0,   65.0,   55.0,   45.0,   37.5,   32.5, &
                      27.5,    22.5,    17.5,    12.5,     8.5,    6.0,    4.5,    3.5,    2.5,   1.75, &
                      1.25,    0.85,     0.6,    0.45,    0.35,   0.25,  0.175,  0.125,  0.085,   0.06, &
                     0.045,   0.035,   0.025,  0.0175,  0.0125, 0.0085,  0.006, 0.0045, 0.0035, 0.0025, &
                   0.00175, 0.00125,  0.0009, 0.00065,  0.0004, 0.0002, 0.0000 /)
             
             ! half-level pressure coefficient b
             b_cmip6_o3(:) = 0.

             CALL Init(LeviX_cmip6_o3, nlev_cmip6_o3, a_cmip6_o3, b_cmip6_o3, status)
             IF_NOTOK_RETURN(status=1)

             ! Calculate weights for interpolation
             ! from CMIP6 to model latitudes.
             ! Note that in IFS the Cartesian coordinate
             ! sin(lat) is used to linearly interpolate
             ! to the model's Gaussian grid
             ! (see e.g. cmip6_piaer_mxr_interp.F90).
             ! For the regular grid of TM5,
             ! it is more appropriate to use
             ! the latitude coordinate directly.
             ! The difference is negligible anyway.
             ! The alternative would be to impose
             ! local mass conservation,
             ! but that is not necessary here.
             lat_cmip6_o3(:) = &
                   (/-90., -88.10526, -86.21053, -84.31579, -82.42105, -80.52631, -78.63158, &
                -76.73684, -74.8421 , -72.94736, -71.05264, -69.1579, -67.26316, &
                -65.36842, -63.47368, -61.57895, -59.68421, -57.78947, -55.89474, -54., &
                -52.10526, -50.21053, -48.31579, -46.42105, -44.52632, -42.63158, &
                -40.73684, -38.84211, -36.94737, -35.05263, -33.15789, -31.26316, &
                -29.36842, -27.47368, -25.57895, -23.68421, -21.78947, -19.89474, -18., &
                -16.10526, -14.21053, -12.31579, -10.42105, -8.526316, -6.631579, &
                -4.736842, -2.842105, -0.9473684, 0.9473684, 2.842105, 4.736842, &
                 6.631579, 8.526316, 10.42105, 12.31579, 14.21053, 16.10526, 18., 19.89474, &
                 21.78947, 23.68421, 25.57895, 27.47368, 29.36842, 31.26316, 33.15789, &
                 35.05263, 36.94737, 38.84211, 40.73684, 42.63158, 44.52632, 46.42105, &
                 48.31579, 50.21053, 52.10526, 54., 55.89474, 57.78947, 59.68421, 61.57895, &
                 63.47368, 65.36842, 67.26316, 69.1579, 71.05264, 72.94736, 74.8421, &
                 76.73684, 78.63158, 80.52631, 82.42105, 84.31579, 86.21053, 88.10526, 90./)

             jlow_cmip6 = 1
             DO j=1,jm(1)
                tm5_lat=lli_z(1)%lat_deg(j)     ! lat in degrees
                ! Since CMIP6 mid-point latitudes run from -90 to +90,
                ! the target latitudes of the model always fall within this range
                DO j_cmip6=jlow_cmip6,nlat_cmip6_o3
                   IF ( (tm5_lat .ge. lat_cmip6_o3(j_cmip6)) .and. &
                        (tm5_lat .lt. lat_cmip6_o3(j_cmip6+1)) ) THEN
                      jlow_cmip6 = j_cmip6
                      EXIT
                   ENDIF
                ENDDO
                jlow_cmip62tm5(j)=jlow_cmip6
                weight_cmip62tm5(j) = (tm5_lat - lat_cmip6_o3(jlow_cmip6)) / &
                                      (lat_cmip6_o3(jlow_cmip6+1) - lat_cmip6_o3(jlow_cmip6))
             ENDDO

             deallocate(a_cmip6_o3)
             deallocate(b_cmip6_o3) 
             deallocate(lat_cmip6_o3)

          ENDIF

          IF  (use_MSR) THEN
             ! Define vertical grid of MSR input files (MSR2 only on 60 levels)
             ecstring='ec60'

             CALL Init(LeviX_msr, ecstring, status)
             IF_NOTOK_RETURN(status=1)
          ENDIF

          IF (use_HALOE) THEN
             ! Define vertical grid of HALOE input files
             select case (lme)
               case(60)
                 ! offline TM5 resolution
                 ecstring='ec60'
               case(34, 91)
                 ! standard EC-Earth3 resolution
                 ecstring='ec91'
               case(40)
                 ecstring='ec40'
               case(62)
                 ecstring='ec62'
               case default
                 write (gol,'("ERROR - HALOE input files not available for this vertical resolution")') ; call goErr
                 status=1; TRACEBACK; return
             end select

             CALL Init(LeviX_haloe, ecstring, status)
             IF_NOTOK_RETURN(status=1)
          ENDIF

       ENDIF

       ! Allocate 
       CALL Get_DistGrid( dgrid(1), J_STRT=j1, J_STOP=j2 ) ! grid size

       lmr = lm(1)

       IF (use_o3du) THEN
          ALLOCATE( o3du(1)%d23(j1:j2, lmr) )
          o3du(1)%d23  = 0.0
       END IF

       allocate( o3ratpm      (1)%d23( j1:j2, lmr) )       
       allocate( o3rat        (1)%d23( j1:j2, lmr) )
       allocate( o3ratnm      (1)%d23( j1:j2, lmr) ) 
       allocate( ch4ratpm     (1)%d23( j1:j2, lmr) )
       allocate( ch4rat       (1)%d23( j1:j2, lmr) )
       allocate( ch4ratnm     (1)%d23( j1:j2, lmr) )
       allocate( odin_hno3_o3 (1)%d23( j1:j2,   4) )
       allocate( odin_hno3_o3_pm(1)%d23( j1:j2, 4) )
       allocate( odin_hno3_o3_nm(1)%d23( j1:j2, 4) )
       allocate( o3vmrpm      (1)%d23( j1:j2, lmr) )     
       allocate( o3vmr        (1)%d23( j1:j2, lmr) )
       allocate( o3vmrnm      (1)%d23( j1:j2, lmr) )
       allocate( ch4vmrpm     (1)%d23( j1:j2, lmr) )
       allocate( ch4vmr       (1)%d23( j1:j2, lmr) ) ! only if Haloe
       allocate( ch4vmrnm     (1)%d23( j1:j2, lmr) )
       allocate( odin_co_o3   (1)%d23( j1:j2,   4) )
       allocate( odin_co_o3_pm(1)%d23( j1:j2,   4) )
       allocate( odin_co_o3_nm(1)%d23( j1:j2,   4) )

       ! work arrays
       IF (isRoot) THEN
          allocate( wrk3dpm_ML  (1, jm(1), lmr) )
          allocate( wrk3d_ML    (1, jm(1), lmr) )
          allocate( wrk3dnm_ML  (1, jm(1), lmr) )
          allocate( wrk3dratio(1, jm(1),   4) )
          allocate( wrk3dratiopm(1, jm(1),   4) )         
          allocate( wrk3drationm(1, jm(1),   4) )        
          allocate( wrk2d_ML    (   jm(1), lmr) )

          allocate( wrk2d_4L  (   jm(1),   4) )
       ELSE
          allocate( wrk3dpm_ML  (1,1,1) )
          allocate( wrk3d_ML    (1,1,1) )
          allocate( wrk3dnm_ML  (1,1,1) )
          allocate( wrk3dratio(1,1,1) )
          allocate( wrk3dratiopm(1,1,1) )
          allocate( wrk3drationm(1,1,1) )

          ALLOCATE( wrk2d_ML(1,1) )
          ALLOCATE( wrk2d_4L(1,1) )
       END IF

       ! safety
       o3ratpm(1)%d23 = 0.0         ; o3rat(1)%d23 = 0.0        ; o3ratnm(1)%d23 = 0.0
       ch4ratpm(1)%d23 = 0.0        ; ch4rat(1)%d23 = 0.0       ; ch4ratnm(1)%d23 = 0.0
       odin_hno3_o3_pm(1)%d23 = 0.0 ; odin_hno3_o3(1)%d23 = 0.0 ; odin_hno3_o3_nm(1)%d23 = 0.0       
       odin_co_o3_pm(1)%d23 = 0.0   ; odin_co_o3  (1)%d23 = 0.0 ; odin_co_o3_nm(1)%d23 = 0.0

       ! Budgets init
       !--------------------------------------------------    
#ifdef with_budgets
       ! sum_stratosphere(1) = 0.0
       budstratg(:,:,:)         = 0.0
#endif

       ! Timers
       call GO_Timer_Def( itim_init, 'boundary init', status )
       IF_NOTOK_RETURN(status=1)
       call GO_Timer_Def( itim_appl, 'boundary appl', status )
       IF_NOTOK_RETURN(status=1)

    ELSE

       ! start timing
       call GO_Timer_Start( itim_init, status )
       IF_NOTOK_RETURN(status=1)

       !--------------------------------------------------
       ! ** OZONE CLIMATOLOGIES
       !--------------------------------------------------

       !--------------------------------------------------
       ! ** o3 vmr
       !--------------------------------------------------
       ROOT : IF (isRoot) THEN

          IF (LCMIP6_O3) THEN
             ! CMIP6 ozone files
             ! containing float vmro3(time, plev, lat, lon) 
             ! with 14 months

             allocate( field4d_cmip6_h(nlon_cmip6_o3, nlat_cmip6_o3, nlev_cmip6_o3,3) ) 
             allocate( field3d_cmip6_pm  (nlat_cmip6_o3, nlev_cmip6_o3) )
             allocate( field3d_cmip6     (nlat_cmip6_o3, nlev_cmip6_o3) )
             allocate( field3d_cmip6_nm  (nlat_cmip6_o3, nlev_cmip6_o3) )
             allocate( field3d_pm        (1, jm(1), nlev_cmip6_o3) )
             allocate( field3d           (1, jm(1), nlev_cmip6_o3) )
             allocate( field3d_nm        (1, jm(1), nlev_cmip6_o3) )

             if (o3_piclim) then
                ! pre-industrial control using 1850 O3 climatology
                write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_clim_1850.nc'
             else 
                if (.not.o3_fixyear) then
                   target_year = MIN(2100,MAX(idate(1),1850))
                else
                   target_year = MIN(2100,MAX(o3_year,1850))
                endif 

                ! To prevent problems with the last entry,
                ! containing the field for January 2015,
                ! O3 data are read from the scenario files
                ! already when the target_year is 2014.
                if ( target_year .gt. 2013 ) then
                   write (gol,'("ERROR - Nudging of stratospheric O3 mixing ratios")') ; call goPr
                   write (gol,'("ERROR - not implemented yet for CMIP6 scenarios")') ; call goErr
                   status=1; TRACEBACK; return
                   ! include case statement to use different filename for each scenario
                   write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'xxxvmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_',target_year,'.nc'
                else 
                   write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_',target_year,'.nc'
                endif
             endif

             CALL MDF_Open( TRIM(o3vmr_fname), MDF_NETCDF, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Inq_VarID( hid,'vmro3', varid, status )
             IF_NOTOK_RETURN(status=1)

             ! Ozone input files contain 14 month entries,
             ! 1: December previous year
             ! 2-13: January to December current year
             ! 14: January next year
             ! Besides the current month, 
             ! also the previous and next month are needed.
             ! By default, these have indices idate(2), idate(2)+1, and idate(2)+2.
             CALL MDF_Get_Var( hid, varid, field4d_cmip6_h, status, start=(/1,1,1,idate(2)/), count=(/nlon_cmip6_o3,nlat_cmip6_o3,nlev_cmip6_o3,3/) )
             IF_NOTOK_RETURN(status=1)

             ! When a fixed year is used for O3
             ! information from the previous or next year
             ! should not be used and overwritten
             ! with the data for the target year.
             ! This is not necessary 
             ! when the pre-industrial climatology is used.
             IF (o3_fixyear .and. .not.o3_piclim) THEN
                IF (idate(2).eq.1) THEN
                   ! overwrite previous month with field for December of the target year
                   CALL MDF_Get_Var( hid, varid, field4d_cmip6_h(:,:,:,1), status, start=(/1,1,1,13/), count=(/nlon_cmip6_o3,nlat_cmip6_o3,nlev_cmip6_o3,1/) )
                ELSE IF (idate(2).eq.12) THEN
                   ! overwrite next month with field for January of the current year
                   CALL MDF_Get_Var( hid, varid, field4d_cmip6_h(:,:,:,3), status, start=(/1,1,1,2/), count=(/nlon_cmip6_o3,nlat_cmip6_o3,nlev_cmip6_o3,1/) )
                ENDIF
             ENDIF

             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1)

             ! calculate zonal mean
             DO k=1,nlev_cmip6_o3
                DO j=1,nlat_cmip6_o3
                   field3d_cmip6_pm(j,k)= sum(field4d_cmip6_h(:,j,k,1))/float(nlon_cmip6_o3)
                   field3d_cmip6(j,k)   = sum(field4d_cmip6_h(:,j,k,2))/float(nlon_cmip6_o3)
                   field3d_cmip6_nm(j,k)= sum(field4d_cmip6_h(:,j,k,3))/float(nlon_cmip6_o3)
                ENDDO
             ENDDO

             deallocate( field4d_cmip6_h )

             DO j=1,jm(1)
               jlow_cmip6=jlow_cmip62tm5(j)
               field3d_pm(1,j,:) = & 
                    field3d_cmip6_pm(jlow_cmip6,:) + weight_cmip62tm5(j) * &
                  ( field3d_cmip6_pm(jlow_cmip6+1,:) - field3d_cmip6_pm(jlow_cmip6,:) )
               field3d(1,j,:) = &
                    field3d_cmip6(jlow_cmip6,:) + weight_cmip62tm5(j) * &
                  ( field3d_cmip6(jlow_cmip6+1,:) - field3d_cmip6(jlow_cmip6,:) )
               field3d_nm(1,j,:) = &
                    field3d_cmip6_nm(jlow_cmip6,:) + weight_cmip62tm5(j) * &
                  ( field3d_cmip6_nm(jlow_cmip6+1,:) - field3d_cmip6_nm(jlow_cmip6,:) )
             ENDDO

             deallocate (field3d_cmip6_pm)
             deallocate (field3d_cmip6) 
             deallocate (field3d_cmip6_nm)

             ! convert from mol/mol to ppmv,
             ! the unit used in the MSR input files
             field3d_pm = 1.e6 * field3d_pm
             field3d    = 1.e6 * field3d
             field3d_nm = 1.e6 * field3d_nm
             
             ! Distribute to model vertical resolution
             src_lev = leviX_cmip6_o3
             rgtype  = 'mass-aver'

             ALLOCATE( sp_z(1,jm(1)) ) ! dummy surface pressure:
             sp_z = p0                 ! 1e5 Pa

             CALL FillLevels( levi,    'n',     sp_z,   wrk3dpm_ML, &
                              src_lev, field3d_pm, rgtype, status    )

             CALL FillLevels( levi,    'n',     sp_z,   wrk3d_ML, &
                              src_lev, field3d, rgtype, status    )

             CALL FillLevels( levi,    'n',     sp_z,   wrk3dnm_ML, &
                              src_lev, field3d_nm, rgtype, status    )

          ELSE 

            IF (use_MSR) THEN

              inp_levs = LeviX_msr%nlev

              allocate( field3D_pm  (1, nlat180, inp_levs) )
              allocate( field3D     (1, nlat180, inp_levs) )
              allocate( field3D_nm  (1, nlat180, inp_levs) )
              allocate( field3d_h   (1, nlat180, inp_levs) )

              field3d_pm=0.0 ; field3d=0.0 ; field3d_nm=0.0 ; field3d_h=0.0

              nmw=idate(2) ! pick out right month in file

              ! filename (includes year and nb of levels)
              if (.not.o3_fixyear) then
                valid_year = MIN(msr2_valid(2), MAX(idate(1), msr2_valid(1)))
                valid_year_pm = MAX(valid_year-1, msr2_valid(1))
                valid_year_nm = MIN(valid_year+1, msr2_valid(2))
              else
                valid_year = MIN(msr2_valid(2), MAX(o3_year, msr2_valid(1))) 
                valid_year_pm = valid_year
                valid_year_nm = valid_year
              endif
              ! Currently available: 1979-2012
              WRITE(o3vmr_fnamepm,'(a,"msr2_o3vmr",i4,"_",i2,".hdf")')TRIM(o3vmr_dirname),valid_year_pm,inp_levs
              WRITE(o3vmr_fname,'  (a,"msr2_o3vmr",i4,"_",i2,".hdf")')TRIM(o3vmr_dirname),valid_year,   inp_levs
              WRITE(o3vmr_fnamenm,'(a,"msr2_o3vmr",i4,"_",i2,".hdf")')TRIM(o3vmr_dirname),valid_year_nm,inp_levs
              !
              ! January
              !
              if(nmw .eq. 1) then

                CALL MDF_Open( TRIM(o3vmr_fnamepm), MDF_HDF4, MDF_READ, hid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Inq_VarID( hid, 'ozone12', varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Close( hid, status )
                IF_NOTOK_RETURN(status=1)
                ! Flip the data over (upside down), since leviX_msr, associated with
                ! field3d when regridding below (Fill3D), is flipped w/r/t to levi
                ! [it is the case, since ecstring starts with "ec" and not "tm"
                ! --see grid_type_hyb.F90--]
                DO k=1,inp_levs
                   field3d_pm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                CALL MDF_Open( TRIM(o3vmr_fname), MDF_HDF4, MDF_READ, hid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Inq_VarID( hid, 'ozone1', varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)     

                DO k=1,inp_levs
                   field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                CALL MDF_Inq_VarID( hid, 'ozone2', varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Close( hid, status )
                IF_NOTOK_RETURN(status=1)       

                DO k=1,inp_levs
                   field3d_nm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                src_lev = leviX_msr
                rgtype  = 'mass-aver'

                deallocate( field3d_h )
              endif
              !
              ! Here the previous and next monthly fields are located in the same year
              !
              if(nmw .gt. 1 .and. nmw .lt. 12) then
                ! field name (has month: ozone1, ozone2, ..., ozone12)
                IF (nmw.LT.11) THEN
                   nmw=idate(2)-1
                   WRITE(filemon,'(a,i1)')'ozone',nmw
                ELSE
                   nmw=idate(2)-1
                   WRITE(filemon,'(a,i2)')'ozone',nmw
                ENDIF

                CALL MDF_Open( TRIM(o3vmr_fname), MDF_HDF4, MDF_READ, hid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                ! Flip the data over (upside down), since leviX_msr, associated with
                ! field3d when regridding below (Fill3D), is flipped w/r/t to levi
                ! [it is the case, since ecstring starts with "ec" and not "tm"
                ! --see grid_type_hyb.F90--]
                DO k=1,inp_levs
                   field3d_pm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                nmw=idate(2)

                ! field name (has month: ozone1, ozone2, ..., ozone12)
                IF (nmw.LT.10) THEN
                   WRITE(filemon,'(a,i1)')'ozone', nmw
                ELSE
                   WRITE(filemon,'(a,i2)')'ozone', nmw
                ENDIF

                CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)      

                DO k=1,inp_levs
                   field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                ! field name (has month: ozone1, ozone2, ..., ozone12)
                IF (nmw.LT.9) THEN
                   nmw=idate(2)+1
                   WRITE(filemon,'(a,i1)')'ozone',nmw
                ELSE
                   nmw=idate(2)+1
                   WRITE(filemon,'(a,i2)')'ozone',nmw
                ENDIF

                CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Close( hid, status )
                IF_NOTOK_RETURN(status=1)

                DO k=1,inp_levs
                   field3d_nm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                src_lev = leviX_msr
                rgtype  = 'mass-aver'

                deallocate( field3d_h )

                nmw=idate(2)
              endif
              !
              ! December
              !
              if(nmw .eq. 12) then
                nmw=idate(2)-1
                WRITE(filemon,'(a,i2)')'ozone',nmw

                CALL MDF_Open( TRIM(o3vmr_fname), MDF_HDF4, MDF_READ, hid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Inq_VarID( hid,TRIM(filemon), varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                ! Flip the data over (upside down), since leviX_msr, associated with
                ! field3d when regridding below (Fill3D), is flipped w/r/t to levi
                ! [it is the case, since ecstring starts with "ec" and not "tm"
                ! --see grid_type_hyb.F90--]
                DO k=1,inp_levs
                   field3d_pm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO
                nmw=12
                WRITE(filemon,'(a,i2)')'ozone',nmw
                CALL MDF_Inq_VarID( hid,TRIM(filemon) , varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Close( hid, status )
                IF_NOTOK_RETURN(status=1)       

                DO k=1,inp_levs
                   field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                nmw=1
                WRITE(filemon,'(a,i1)')'ozone',nmw
                CALL MDF_Open( TRIM(o3vmr_fnamenm), MDF_HDF4, MDF_READ, hid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
                IF_NOTOK_RETURN(status=1)
                CALL MDF_Close( hid, status )
                IF_NOTOK_RETURN(status=1)       

                DO k=1,inp_levs
                   field3d_nm(1,:,k)=field3d_h(1,:,inp_levs+1-k)
                ENDDO

                src_lev = leviX_msr
                rgtype  = 'mass-aver'

                deallocate( field3d_h )
              endif
             
            ELSE ! use other file already remapped to model levels (e.g. o3vmr2000_tropo25.hdf)

              allocate( field3d(1, nlat180, lm(1) ))
              field3d=0.0

              o3vmr_fname=o3vmr_dirname

              CALL MDF_Open( TRIM(o3vmr_fname), MDF_HDF4, MDF_READ, hid, status )
              IF_NOTOK_RETURN(status=1)
              CALL MDF_Get_Var( hid, idate(2), field3d(1,:,:), status )
              IF_NOTOK_RETURN(status=1)
              CALL MDF_Close( hid, status )
              IF_NOTOK_RETURN(status=1)

             src_lev = levi
             rgtype  = 'area-aver'

            ENDIF

            ! Coarsen/distribute to each 1 resolution

            ALLOCATE( sp_z(1,jm(1)) ) ! dummy surface pressure:
            sp_z = p0                 ! 1e5 Pa

            CALL Fill3D( lli_z(1),       levi,    'n',     sp_z,   wrk3dpm_ML, &
                         lli_z(iglbsfc), src_lev, field3d_pm, rgtype, status    )

            CALL Fill3D( lli_z(1),       levi,    'n',     sp_z,   wrk3d_ML, &
                         lli_z(iglbsfc), src_lev, field3d, rgtype, status    )

            CALL Fill3D( lli_z(1),       levi,    'n',     sp_z,   wrk3dnm_ML, &
                         lli_z(iglbsfc), src_lev, field3d_nm, rgtype, status    )

          ENDIF

          IF_NOTOK_RETURN(status=1)

          deallocate( sp_z )
          deallocate( field3d_pm)
          deallocate( field3d )
          deallocate( field3d_nm)

       END IF ROOT

       ! scatter along latitude direction, then broadcast
       wrk2d_ML = wrk3dpm_ML(1,:,:) ! reshape

       CALL SCATTER_I_BAND( dgrid(1), o3vmrpm(1)%d23, wrk2d_ML, status )
       IF_NOTOK_RETURN(status=1) 

       CALL PAR_BROADCAST( o3vmrpm(1)%d23, status, row=.TRUE. )
       IF_NOTOK_RETURN(status=1)

       wrk2d_ML = wrk3d_ML(1,:,:) ! reshape

       CALL SCATTER_I_BAND( dgrid(1), o3vmr(1)%d23, wrk2d_ML, status )
       IF_NOTOK_RETURN(status=1)

       CALL PAR_BROADCAST( o3vmr(1)%d23, status, row=.TRUE. )
       IF_NOTOK_RETURN(status=1)

       wrk2d_ML = wrk3dnm_ML(1,:,:) ! reshape

       CALL SCATTER_I_BAND( dgrid(1), o3vmrnm(1)%d23, wrk2d_ML, status )
       IF_NOTOK_RETURN(status=1) 

       CALL PAR_BROADCAST( o3vmrnm(1)%d23, status, row=.TRUE. )
       IF_NOTOK_RETURN(status=1)

       !--------------------------------------------------
       ! ** O3 DU
       !--------------------------------------------------
       IF (use_o3du) THEN

          IF(isRoot) THEN

             ! o3du only required for 40 or 62 levels
             inp_levs = lme

             ALLOCATE( field3d  (1, nlat180, inp_levs) )
             ALLOCATE( field3d_h(1, nlat180, inp_levs) )

             nmw=idate(2) ! pick out right month in file

             WRITE(o3du_fname,'(a,"msr2_o3du",i4,"_",i2,".hdf")')TRIM(o3du_dirname),idate(1),inp_levs
             IF (nmw.LT.10) THEN
                WRITE(filemon,'(a,i1)')'ozone', idate(2)
             ELSE
                WRITE(filemon,'(a,i2)')'ozone', idate(2)
             ENDIF

             CALL MDF_Open( TRIM(o3du_fname), MDF_HDF4, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, varid, field3d_h, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1)       

             ! flip data so that they are on the leviX levels
             DO k=1,inp_levs
                field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)
             ENDDO

             ! remap to model levels
             ALLOCATE( sp_z(1,jm(1)) )
             sp_z = p0  ! 1e5 Pa

             ! Should we really use leviX_msr here?
             ! Doesn't seem to match with the use of lme above.
             ! Anyway, o3du is normally not used anymore,
             ! since the vertical resolution has increased.
             CALL Fill3D( lli_z(1),       levi, 'n',       sp_z,       wrk3d_ML, &
                          lli_z(iglbsfc), leviX_msr, field3d, 'area-aver', status       )
             IF_NOTOK_RETURN(status=1)

             DEALLOCATE( sp_z )
             DEALLOCATE( field3d )
             DEALLOCATE( field3d_h )

          END IF

          ! scatter along latitude direction, then broadcast
          wrk2d_ML = wrk3d_ML(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), o3du(1)%d23, wrk2d_ML, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( o3du(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

       ENDIF ! use_o3du

       !--------------------------------------------------
       ! ** CH4 HALOE
       !--------------------------------------------------
       IF (use_HALOE) THEN

          IF(isRoot) THEN

             inp_levs = LeviX_haloe%nlev

             if (LCMIP6_CH4) then
               allocate(ch4_hist(1,2015)) ! years 0-2014
!               WRITE(ch4_hist_fname,'(a,"mole_fraction_of_methane_in_air_input4MIPs_GHGConcentrations_CMIP_UoM-CMIP-1-2-0_gr1-GMNHSH_0000-2014.nc")')TRIM(cmip6_ch4_dirname_strat)
               WRITE(ch4_hist_fname,'(a,"CMIP6_histo_mole_fraction_of_methane_in_air_input4MIPs_gr1-GMNHSH.nc")')TRIM(cmip6_ch4_dirname_strat)
               CALL MDF_Open( TRIM(ch4_hist_fname), MDF_NETCDF, MDF_READ, hid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Inq_VarID( hid, TRIM('mole_fraction_of_methane_in_air'), varid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Get_Var( hid, varid, ch4_hist(1,:), status, start = (/1,1/), count = (/1,2015/) )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Close( hid, status )
               IF_NOTOK_RETURN(status=1)

               ! Following the recommendation by Meinshausen et al. (GMDD, 2016),
               ! a one year delay is assumed between the surface and the stratosphere.
               ! The HALOE climatology is based on measurements 
               ! from October 1991 to August 2002.
               ! HALOE therefore provides the full annual cycle 
               ! for the 10-year period 1992-2001.
               ! Since our scaling is based on annual mean concentrations,
               ! we take this period as our reference.
               ! To account for the delay of one year,
               ! this translates into 1991-2000 at the surface.
               ch4_ref=0.
               do iyear=1991,2000
                 ! First year of ch4_hist is the year 0
                 ch4_ref=ch4_ref+ch4_hist(1,iyear+1)
               enddo
               ch4_ref=ch4_ref/10.

               if ( (  (.not.ch4_fixyear) .and. &
                       ( (idate(1) .gt. 2015)  .or. (idate(1)==2015 .and. idate(2)==12) ) ) &
                    .or. &
                    (ch4_fixyear .and. ch4_year.gt.2014) ) then
                 write (gol,'("ERROR - Scaling of stratospheric CH4 concentrations from HALOE")') ; call goPr
                 write (gol,'("ERROR - not implemented yet for CMIP6 scenarios")') ; call goErr
                 status=1; TRACEBACK; return

                 allocate(ch4_sce(1,86)) ! years 2015-2100
                 ! include case statement to use different filename for each scenario
                 WRITE(ch4_sce_fname,'(a,"xxxx")')TRIM(cmip6_ch4_dirname_strat)
                 CALL MDF_Open( TRIM(ch4_sce_fname), MDF_NETCDF, MDF_READ, hid, status )
                 IF_NOTOK_RETURN(status=1)
                 CALL MDF_Inq_VarID( hid, TRIM('mole_fraction_of_methane_in_air'), varid, status )
                 IF_NOTOK_RETURN(status=1)
                 CALL MDF_Get_Var( hid, varid, ch4_sce(1,:), status, start = (/1,1/), count = (/1,86/) )
                 IF_NOTOK_RETURN(status=1)
                 CALL MDF_Close( hid, status )
                 IF_NOTOK_RETURN(status=1)
               endif

               ! Following CMIP methodology,
               ! surface concentrations before 1850
               ! are set to their 1850 level.
               !
               ! target_year is the year used in the scaling
               ! of surface concentrations.

               ! scale factor for previous month:
               if (.not.ch4_fixyear) then
                 target_year = MIN(2100,MAX(idate(1)-1,1850))
                 if(idate(2) .eq. 1) then
                   target_year = MIN(2100,MAX(idate(1)-2,1850))
                 endif
               else
                 target_year = MIN(2100,MAX(ch4_year,1850)) 
               endif
               write (gol,*) 'Stratospheric CH4 year for previous month: ', target_year; call goPr
               if (target_year .le. 2014) then
                 ch4_scale_pm=ch4_hist(1,target_year+1)/ch4_ref
               else
                 ch4_scale_pm=ch4_sce(1,target_year-2014)/ch4_ref
               endif

               ! scale factor for the current month:
               if (.not.ch4_fixyear) then
                 target_year = MIN(2100,MAX(idate(1)-1,1850))
               else
                 target_year = MIN(2100,MAX(ch4_year,1850))
               endif
               write (gol,*) 'Stratospheric CH4 year for current month: ', target_year; call goPr
               if (target_year .le. 2014) then
                 ch4_scale=ch4_hist(1,target_year+1)/ch4_ref
               else
                 ch4_scale=ch4_sce(1,target_year-2014)/ch4_ref
               endif

               ! scale factor for the next month: 
               if (.not.ch4_fixyear) then
                 target_year = MIN(2100,MAX(idate(1)-1,1850))
                 if(idate(2) .eq. 12) then
                   target_year = MIN(2100,MAX(idate(1),1850))
                 endif
               else
                 target_year = MIN(2100,MAX(ch4_year,1850))
               endif
               write (gol,*) 'Stratospheric CH4 year for next month: ', target_year; call goPr 
               if (target_year .le. 2014) then
                 ch4_scale_nm=ch4_hist(1,target_year+1)/ch4_ref
               else
                 ch4_scale_nm=ch4_sce(1,target_year-2014)/ch4_ref
               endif 

               deallocate(ch4_hist)
               if (allocated(ch4_sce)) deallocate(ch4_sce)

             else 

               ! The pre-CMIP6 calculation of the scale factors is not fully
               ! consistent with the new implementation for CMIP6.
               ! This could be changed, but I haven't done that.
               ! I would propose to use LCMIP6_CH4 by default, 
               ! and remove the scaling based on RCP6 from the code. 
               write (gol,'("WARNING - Scaling of stratospheric CH4 concentrations from HALOE")') ; call goPr
               write (gol,'("WARNING - based on old implementation, which uses RCP6 for years later than 2005.")') ; call goPr

               if (ch4_fixyear) then
                  write (gol,'("ERROR - Fixing CH4 concentrations to a specific year")') ; call goErr
                  write (gol,'("ERROR - only works for CMIP6 data")') ; call goErr
                  status=1; TRACEBACK; return
               endif

               ! The scale factor for the previous month
               ! used to be based on the previous year
               ! irrespective of the month.
               ! This is incorrect and has been fixed:
               !valid_year = MIN(2035,MAX(idate(1)-1,1999)) 
               valid_year = MIN(2035,MAX(idate(1),1999))
               if(idate(2) .eq. 1) then
                 valid_year = MIN(2035,MAX(idate(1)-1,1999))
               endif
               ch4_scale_pm=RCP6_CH4_STRAT(valid_year-1998)

               ! scale factor for the current month
               valid_year = MIN(2035,MAX(idate(1),1999))
               ch4_scale=RCP6_CH4_STRAT(valid_year-1998)

               ! scale factor for the next month
               valid_year = MIN(2035,MAX(idate(1),1999))
               if(idate(2) .eq. 12) then
                 valid_year = MIN(2035,MAX(idate(1)+1,1999))
               endif
               ch4_scale_nm=RCP6_CH4_STRAT(valid_year-1998)

             endif

             allocate( field3d_pm  (1, nlat180, inp_levs) )             
             allocate( field3d     (1, nlat180, inp_levs) )
             allocate( field3d_nm  (1, nlat180, inp_levs) )
             allocate( field3d_h   (1, nlat180, inp_levs) )

             field3d_pm=0.0 ; field3d=0.0 ; field3d_nm=0.0 ; field3d_h=0.0
             !
             ! For each month, data are read, then flipped over (upside down), otherwise Fill3D
             ! doesn't work properly, and scaled up to account for growth from 2000.             
             !
             ! - Previous month
             !
             if(idate(2) .eq. 1) then
                nmw=12
                WRITE(filemon,'(a,i2)')'haloe',nmw
             else
                nmw=idate(2)-1
                if(idate(2) .lt. 11) WRITE(filemon,'(a,i1)')'haloe',nmw
                if(idate(2) .gt. 10) WRITE(filemon,'(a,i2)')'haloe',nmw
             endif
             !
             WRITE(ch4vmr_fname,'(a,"haloe_ch4vmr_",i2,".hdf")')TRIM(ch4vmr_dirname),inp_levs

             CALL MDF_Open( TRIM(ch4vmr_fname), MDF_HDF4, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
             IF_NOTOK_RETURN(status=1)
             
             DO k=1,inp_levs
                field3d_pm(1,:,k)=field3d_h(1,:,inp_levs+1-k)*ch4_scale_pm
             ENDDO

             !
             ! - Current month
             !
             nmw=idate(2)

             IF (nmw.LT.10) THEN
                WRITE(filemon,'(a,i1)')'haloe', nmw
             ELSE
                WRITE(filemon,'(a,i2)')'haloe', nmw
             ENDIF

             CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
             IF_NOTOK_RETURN(status=1)      

             DO k=1,inp_levs
                field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)*ch4_scale
             ENDDO

             !
             ! - Next month
             !             
             nmw=idate(2)+1
             if(idate(2) .eq. 12) then
                nmw=1
             endif

             if(nmw .lt. 10) WRITE(filemon,'(a,i1)')'haloe', nmw
             if(nmw .gt. 9)  WRITE(filemon,'(a,i2)')'haloe', nmw

             CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, varid, field3d_h(1,:,:), status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1)       

             DO k=1,inp_levs
                field3d_nm(1,:,k)=field3d_h(1,:,inp_levs+1-k)*ch4_scale_nm
             ENDDO
             
             !
             ! Coarsen or distribute on region #1 the three datasets
             !
             ALLOCATE( sp_z(1,jm(1)) )
             sp_z = p0  ! 1e5 Pa

             CALL Fill3D( lli_z(1),       levi,  'n',     sp_z,        wrk3dpm_ML, &
                          lli_z(iglbsfc), leviX_haloe, field3d_pm, 'mass-aver', status    )            

             CALL Fill3D( lli_z(1),       levi,  'n',     sp_z,        wrk3d_ML, &
                          lli_z(iglbsfc), leviX_haloe, field3d, 'mass-aver', status    )

             CALL Fill3D( lli_z(1),       levi,  'n',     sp_z,        wrk3dnm_ML, &
                          lli_z(iglbsfc), leviX_haloe, field3d_nm, 'mass-aver', status    ) 

             IF_NOTOK_RETURN(status=1)

             ! clear allocatables
             deallocate( sp_z )
             deallocate( field3d_pm )
             deallocate( field3d )
             deallocate( field3d_nm )
             deallocate( field3d_h )

          END IF ! Read & Regrid on Root

          ! Scatter along latitude direction, then broadcast to other cores in same zonal band

          wrk2d_ML = wrk3dpm_ML(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), ch4vmrpm(1)%d23, wrk2d_ML, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( ch4vmrpm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

          wrk2d_ML = wrk3d_ML(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), ch4vmr(1)%d23, wrk2d_ML, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( ch4vmr(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

          wrk2d_ML = wrk3dnm_ML(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), ch4vmrnm(1)%d23, wrk2d_ML, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( ch4vmrnm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

       ENDIF  ! use_HALO
       
       !--------------------------------------------------
       ! ** ODIN climatology : HNO3/O3 and CO/O3 ratios
       !--------------------------------------------------
       IF (use_ODIN) then

          inp_levs=4

          !--------------------------------------------------
          ! * ODIN climatology : HNO3/O3
          !--------------------------------------------------
          IF ( isRoot ) THEN

             ! open, read ratio
             CALL MDF_Open( TRIM(covmr_dirname)//'/ODIN_Climatology_HNO3_O3_4levels.hdf', MDF_HDF4, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Inq_VarID( hid, 'Ratio', varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid,varid, field3d_r4, status )
             IF_NOTOK_RETURN(status=1)

             pm=idate(2)-1

             ! previous month
             allocate( field3d(1,nlat180,4) )
             if(idate(2) .eq. 1) field3d(1,:,:)=field3d_r4(:,:,12)
             if(idate(2) .gt. 1) field3d(1,:,:)=field3d_r4(:,:,pm)

             ! remap
             DO i=1,inp_levs
                CALL FillGrid( lli_z(1),       'n', wrk3dratiopm(:,:,i),  &
                               lli_z(iglbsfc), 'n', field3d(:,:,i), 'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             ! current month
             field3d(1,:,:)=field3d_r4(:,:,idate(2))

             ! remap
             DO i=1,inp_levs
                CALL FillGrid( lli_z(1),       'n', wrk3dratio(:,:,i),  &
                               lli_z(iglbsfc), 'n', field3d(:,:,i), 'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             nm=idate(2)+1

             if(idate(2) .lt. 12) field3d(1,:,:)=field3d_r4(:,:,nm)
             if(idate(2) .eq. 12) field3d(1,:,:)=field3d_r4(:,:,1)

             ! remap             
             DO i=1,inp_levs
                CALL FillGrid( lli_z(1),       'n', wrk3drationm(:,:,i),  &
                               lli_z(iglbsfc), 'n', field3d(:,:,i), 'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1) 

             DEALLOCATE( field3d )

          END IF  ! root

          ! Scatter along latitude direction, then broadcast to other cores in same zonal band

          wrk2d_4L = wrk3dratiopm(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_hno3_o3_pm(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_hno3_o3_pm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)         

          wrk2d_4L = wrk3dratio(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_hno3_o3(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_hno3_o3(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

          wrk2d_4L = wrk3drationm(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_hno3_o3_nm(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_hno3_o3_nm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)             

          !--------------------------------------------------
          ! * ODIN climatology : CO O3
          !--------------------------------------------------
          IF ( isRoot ) THEN

             ! open, read ratio             
             CALL MDF_Open(TRIM(covmr_dirname)//'ODIN_CO_O3_ratio_4levels.hdf', MDF_HDF4, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Inq_VarID( hid, 'CO_O3_ratio', varid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, varid, field3d_r4, status )
             IF_NOTOK_RETURN(status=1)

             allocate( field3d(1,nlat180,4) ) ! 1D array along latitude

             pm=idate(2)-1

             ! previous month
             if(idate(2) .eq. 1) field3d(1,:,:)=field3d_r4(:,:,12)
             if(idate(2) .gt. 1) field3d(1,:,:)=field3d_r4(:,:,pm)

             ! remap
             DO i=1,inp_levs         
                CALL FillGrid( lli_z(1),       'n', wrk3dratiopm(:,:,i),                     &
                               lli_z(iglbsfc), 'n', field3d(:,:,i),       'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             ! current month
             field3d(1,:,:)=field3d_r4(:,:,idate(2))

             ! remap
             DO i=1,inp_levs
                CALL FillGrid( lli_z(1),       'n', wrk3dratio(:,:,i),                     &
                               lli_z(iglbsfc), 'n', field3d(:,:,i),       'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             nm=idate(2)+1

             ! correct month
             if(idate(2) .lt. 12) field3d(1,:,:)=field3d_r4(:,:,nm)
             if(idate(2) .eq. 12) field3d(1,:,:)=field3d_r4(:,:,1)

             ! remap
             DO i=1,inp_levs         
                CALL FillGrid( lli_z(1),       'n', wrk3drationm(:,:,i),                     &
                               lli_z(iglbsfc), 'n', field3d(:,:,i),       'area-aver', status )
                IF_NOTOK_RETURN(status=1)
             END DO

             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1)           
             deallocate( field3d )                  

          END IF  ! root

          ! scatter along latitude direction, then broadcast

          wrk2d_4L = wrk3dratiopm(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_co_o3_pm(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_co_o3_pm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1) 

          wrk2d_4L = wrk3dratio(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_co_o3(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_co_o3(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1) 

          ! scatter along latitude direction, then broadcast

          wrk2d_4L = wrk3drationm(1,:,:)

          CALL SCATTER_I_BAND( dgrid(1), odin_co_o3_nm(1)%d23(:,:), wrk2d_4L, status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_co_o3_nm(1)%d23, status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)  

       ELSE

          !--------------------------------------------------
          ! ** UARS climatology : HNO3 O3 (Fallback on this one if we are not using ODIN)
          !--------------------------------------------------
          IF ( isRoot ) THEN

             ALLOCATE( field3d(1,nlat180,1) ) ! 1D array along latitude
             field3d = 0.0

             CALL MDF_Open( TRIM(inputdir)//'/boundary/HNO3_top/uars_ratio.hdf', MDF_HDF4, MDF_READ, hid, status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Get_Var( hid, idate(2), field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)
             CALL MDF_Close( hid, status )
             IF_NOTOK_RETURN(status=1)       

             ! remap
             CALL FillGrid( lli_z(1),       'n', wrk3dratio(:,:,1), &
                            lli_z(iglbsfc), 'n', field3d(:,:,1),    'area-aver', status )
             IF_NOTOK_RETURN(status=1)

             DEALLOCATE( field3d )

          END IF  ! root

          ! scatter along latitude direction, then broadcast
          wrk2d_ML(:,1) = wrk3dratio(1,:,1)

          CALL SCATTER_I_BAND( dgrid(1), odin_hno3_o3(1)%d23(:,1), wrk2d_ML(:,1), status )
          IF_NOTOK_RETURN(status=1)

          CALL PAR_BROADCAST( odin_hno3_o3(1)%d23(:,1), status, row=.TRUE. )
          IF_NOTOK_RETURN(status=1)

       ENDIF

       ! ------------------------------------------------
       ! Annual mean CO2 mixing ratios for pH calculation 
       ! used in aqueous chemistry.
       !-------------------------------------------------
       ! To avoid MPI communication, 
       ! reading is done on all processors, 
       ! but only once per year.
       !-------------------------------------------------
       if (LCMIP6_CO2 .and. newyr) then

          IF ( isRoot ) THEN

             ! Following CMIP methodology,
             ! CO2 mixing ratios before 1850
             ! are set to their 1850 level.
             !
             if (.not.co2_fixyear) then
               target_year = MIN(2100,MAX(idate(1),1850))
             else
               target_year = MIN(2100,MAX(co2_year,1850))
             endif
             write (gol,*) 'In pH calculation, using CO2 for year: ', target_year; call goPr

             if (target_year .le. 2014) then
               allocate(co2_hist(1,2015)) ! years 0-2014
!               WRITE(co2_hist_fname,'(a,"mole_fraction_of_carbon_dioxide_in_air_input4MIPs_GHGConcentrations_CMIP_UoM-CMIP-1-2-0_gr1-GMNHSH_0000-2014.nc")')TRIM(cmip6_co2_dirname)
               WRITE(co2_hist_fname,'(a,"CMIP6_histo_mole_fraction_of_carbon_dioxide_in_air_input4MIPs_gr1-GMNHSH.nc")')TRIM(cmip6_co2_dirname)
               CALL MDF_Open( TRIM(co2_hist_fname), MDF_NETCDF, MDF_READ, hid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Inq_VarID( hid, TRIM('mole_fraction_of_carbon_dioxide_in_air'), varid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Get_Var( hid, varid, co2_hist(1,:), status, start = (/1,1/), count = (/1,2015/) )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Close( hid, status )
               IF_NOTOK_RETURN(status=1)
               co2_glob=co2_hist(1,target_year+1)
               deallocate(co2_hist)
             else
               write (gol,'("ERROR - Reading of CO2 mixing ratios")') ; call goPr
               write (gol,'("ERROR - not implemented yet for CMIP6 scenarios")') ; call goErr
               status=1; TRACEBACK; return

               allocate(co2_sce(1,86)) ! years 2015-2100
               ! include case statement to use different filename for each scenario
               WRITE(co2_sce_fname,'(a,"xxxx")')TRIM(cmip6_co2_dirname)
               CALL MDF_Open( TRIM(co2_sce_fname), MDF_NETCDF, MDF_READ, hid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Inq_VarID( hid, TRIM('mole_fraction_of_carbon_dioxide_in_air'), varid, status )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Get_Var( hid, varid, co2_sce(1,:), status, start = (/1,1/), count = (/1,86/) )
               IF_NOTOK_RETURN(status=1)
               CALL MDF_Close( hid, status )
               IF_NOTOK_RETURN(status=1)
               co2_glob=co2_sce(1,target_year-2014)
               deallocate(co2_sce)
             endif

          ENDIF ! root

          CALL PAR_BROADCAST( co2_glob, status )
          IF_NOTOK_RETURN(status=1)

          write (gol,*) 'Annual mean CO2 mixing ratio (ppm) = ', co2_glob; call goPr

       ENDIF

       ! end timing:
       call GO_Timer_End( itim_init, status )
       IF_NOTOK_RETURN(status=1)

    END IF

    CALL goLabel()
    status = 0

  END SUBROUTINE BOUNDARY_INIT
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BOUNDARY_DONE
  !
  ! !DESCRIPTION: deallocate arrays (budgets gathering/reduction is done in
  !               chemistry/chemistry_done)
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE BOUNDARY_DONE( status )
    !
    ! !USES:
    !
    USE partools,      ONLY : isRoot
    USE Grid,          ONLY : Done
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !   15 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/Boundary_Done'

    ! --- begin ----------------------------------------

    CALL goLabel(rname)

    IF (isRoot) THEN
       IF (LCMIP6_O3) THEN
          CALL Done( LeviX_cmip6_o3, status )
          IF_NOTOK_RETURN(status=1)
       ENDIF

       IF (use_MSR) THEN
          CALL Done( LeviX_msr, status )
          IF_NOTOK_RETURN(status=1)
       ENDIF

       IF (use_HALOE) THEN
          CALL Done( LeviX_haloe, status )
          IF_NOTOK_RETURN(status=1)
       ENDIF
    ENDIF

    IF (LCMIP6_O3) THEN
     if (allocated(weight_cmip62tm5)) deallocate (weight_cmip62tm5)
     if (allocated(jlow_cmip62tm5))   deallocate (jlow_cmip62tm5)
    ENDIF

    IF (use_o3du) deallocate( o3du(1)%d23  )
     deallocate( o3vmrpm      (1)%d23  )    
     deallocate( o3vmr        (1)%d23  )
     deallocate( o3vmrnm      (1)%d23  )
     deallocate( ch4vmrpm     (1)%d23  )  
     deallocate( ch4vmr       (1)%d23  )
     deallocate( ch4vmrnm     (1)%d23  )
     deallocate( o3ratpm      (1)%d23  )
     deallocate( o3rat        (1)%d23  )
     deallocate( o3ratnm      (1)%d23  ) 
     deallocate( ch4ratpm     (1)%d23  )   
     deallocate( ch4rat       (1)%d23  )
     deallocate( ch4ratnm     (1)%d23  )

     deallocate( odin_hno3_o3   (1)%d23  )
     deallocate( odin_hno3_o3_pm(1)%d23  )
     deallocate( odin_hno3_o3_nm(1)%d23  )
     deallocate( odin_co_o3     (1)%d23  )
     deallocate( odin_co_o3_pm  (1)%d23  )
     deallocate( odin_co_o3_nm  (1)%d23  )
    
     deallocate( wrk3dpm_ML    )    
     deallocate( wrk3d_ML      )
     deallocate( wrk3dnm_ML    ) 
     deallocate( wrk3dratio    )
     deallocate( wrk3dratiopm  )
     deallocate( wrk3drationm  )
     deallocate( wrk2d_ML      )
     deallocate( wrk2d_4L      )

    CALL goLabel()

    status = 0

  END SUBROUTINE BOUNDARY_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    BOUNDARY_APPLY
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE BOUNDARY_APPLY( region, status )
    !
    ! !USES:
    !
    USE dims,        ONLY : im, jm, lm
    USE dims,        ONLY : ybeg
    USE dims,        ONLY : xref, yref, tref
    USE dims,        ONLY : at, bt
    USE dims,        ONLY : dx, dy
    USE dims,        ONLY : nsrce
    USE dims,        ONLY : okdebug
    USE dims,        ONLY : sec_day,sec_month
    USE toolbox,     ONLY : lvlpress
    USE chem_param,  ONLY : io3, io3s, ich4, ihno3, ico
    USE chem_param,  ONLY : xmair, xmhno3, xmo3, xmch4, xmco, fscale, ra
    USE meteodata,   ONLY : m_dat
    USE global_data, ONLY : region_dat, mass_dat
    USE partools,    ONLY : par_reduce_element
    use GO,          ONLY : GO_Timer_End, GO_Timer_Start
    !
    ! !INPUT PARAMETERS:
    !
    INTEGER, INTENT(in)     ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !   23 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER  ::  rname = mname//'/Boundary_Apply'

    ! Nudging coefficients for O3 and CH4
    REAL, PARAMETER :: gnud_pole  = 2.894e-6   !  1 / (4   days)
    REAL, PARAMETER :: gnud_hilat = 3.858e-6   !  1 / (3   days)
    REAL, PARAMETER :: gnud_trop  = 4.630e-6   !  1 / (2.5 days)
    ! Nudging coefficients for HNO3 and CO (tested by B. Maasakkers, see KNMI internship report)
    REAL, PARAMETER :: gnud_hno3_5hpa   = 2.315e-6   !  1 / (5   days)
    REAL, PARAMETER :: gnud_hno3_10hpa  = 5.587e-7   !  1 / (10  days) 
    REAL, PARAMETER :: gnud_hno3_28hpa  = 1.987e-7   !  1 / (60  days)
    REAL, PARAMETER :: gnud_hno3_43hpa  = 2.572e-7   !  1 / (45  days)

    REAL,DIMENSION(:,:,:),   POINTER :: m, ppm
    REAL,DIMENSION(:,:,:,:), POINTER :: rm

    REAL,DIMENSION(:,:), ALLOCATABLE :: znl_lcl, zonal_glbl  ! zonal (total or average), local and global

    INTEGER :: lmr, nlevels, nudge_lev, daystep
    INTEGER :: i, j, l, ipos, nzone, nzone_v, n, lm_min, nday
    REAL    :: pr, dtime, dyy, xlat, gnud, o3za, ch4za, coza
    REAL    :: rmold, rmold1, rmold2, z, z1, z2, rmobs, rmobs1, rmobs2

    INTEGER :: lsave, i1, i2, j1, j2
    REAL    :: fraction_pm, fraction_cm, fraction_nm, daily_scale
    REAL    :: diffmin, pressure(4), gnud_hno3(4), gnud_co(4)
    REAL    :: nudgeFactor, tm5_o3_hno3_mix_ratio, tm5_o3_co_mix_ratio
    !
    !  no. days from 20th month till end
    !
    REAL, DIMENSION(12) :: nday_em = (/11.0, 8.0, 11.0, 10.0, 11.0, 10.0, 11.0, 11.0, 10.0, 11.0, 10.0, 11.0/)
    
    ! --- Begin ----------------------------------

    ! O3, O3S and CH4 are relaxed to climatology at stratospheric layers
    ! HNO3 is prescribed at level lm using ODIN or UARS HNO3:03 ratios
    ! CO is prescribed at level lm using ODIN CO:03 ratios

    ! start timing:
    call GO_Timer_Start( itim_appl, status )
    IF_NOTOK_RETURN(status=1)

    CALL goLabel(rname)

    CALL Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr=lm(region)
    
    dtime=nsrce/(2*tref(region))  
    dyy=dy/yref(region)  !gridsize at current region

    rm => mass_dat(region)%rm
    m  => m_dat   (region)%data

    nday=sec_month/sec_day
    !
    !   Define the fractions from the previous and next months
    !   For the period 10-20th of each month fraction_cm=1.0
    !
    fraction_pm = 0.0 ; fraction_nm = 0.0 ; fraction_cm = 1.0
    if(idate(3) .lt. 10) then
       fraction_pm = 1.0-((idate(3)+nday_em(idate(2)))/((nday_em(idate(2))+10)))
       fraction_cm = 1.0-fraction_pm
    endif
    if(idate(3) .gt. 20) then
       fraction_nm = 0.5*(real(idate(3))-20.5)/(real(nday)-20.0)
       fraction_cm = 1.0-fraction_nm
    endif

    ! ---------------------------------------------------------------------
    ! Get the NUDGING FACTORS
    ! ---------------------------------------------------------------------
    ! Only for the global domain, for which we are sure a good Zonal Average
    ! can be calculated. The nudging factors are given by the ratio between
    ! the climatological zonal mean and the observed zonal mean.
    
    IF (region==1) THEN

       ! Calculate the ratio of the climatology and the zonal mean model field.
       ! The concentration at a given longitude is relaxed towards the ratio of
       ! climatology and zonal mean model field multiplied by the model value 

       ALLOCATE( zonal_glbl(j1:j2, lm(1)) )
       ALLOCATE( znl_lcl   (j1:j2, lm(1)) )
       ALLOCATE( ppm(i1:i2, j1:j2, lm(1)) )

       ! ------ O3
       ppm = 1e6*fscale(io3)*rm(i1:i2,j1:j2,:,io3)/m(i1:i2,j1:j2,:) ! into ppmv       
       znl_lcl=sum(ppm,dim=1) ! total zonal O3 in ppmv (local: i.e. only tile is accounted for)
       
       CALL PAR_REDUCE_ELEMENT(znl_lcl, 'SUM', zonal_glbl, status, all=.true., row=.true.)
       IF_NOTOK_RETURN(status=1)
       
       ! zonal average
       zonal_glbl = zonal_glbl/im(1)
       
       DO l=1,lm(1)
       DO j=j1,j2

          !ratio clim/obs   
          IF ( zonal_glbl(j,l) == 0.0 ) THEN
             o3ratpm(1)%d23(j,l) = 0.0
             o3rat  (1)%d23(j,l) = 0.0
             o3ratnm(1)%d23(j,l) = 0.0
          ELSE
             o3ratpm(1)%d23(j,l) = o3vmrpm(1)%d23(j,l)/zonal_glbl(j,l)
             o3rat  (1)%d23(j,l) = o3vmr(1)  %d23(j,l)/zonal_glbl(j,l)
             o3ratnm(1)%d23(j,l) = o3vmrnm(1)%d23(j,l)/zonal_glbl(j,l)
          END IF
       ENDDO
       ENDDO

       ! ------ CH4
       IF (use_HALOE) THEN  ! else CH4rat remains 0 everywhere

          ppm=1e6*fscale(ich4)*rm(i1:i2,j1:j2,:,ich4)/m(i1:i2,j1:j2,:)
          znl_lcl=sum(ppm,dim=1) ! total zonal (local: i.e. tile only)

          CALL PAR_REDUCE_ELEMENT(znl_lcl, 'SUM', zonal_glbl, status, all=.true., row=.true.)
          IF_NOTOK_RETURN(status=1)

          ! zonal average
          zonal_glbl=zonal_glbl/im(1)
          
          DO l=1,lm(1)
             DO j=j1, j2
                ! climatology/model ratio
                IF ( zonal_glbl(j,l) == 0.0 ) THEN
                   CH4ratpm(1)%d23(j,l) = 0.0                   
                   CH4rat  (1)%d23(j,l) = 0.0
                   CH4ratnm(1)%d23(j,l) = 0.0
                ELSE
                   CH4ratpm(1)%d23(j,l) = ch4vmrpm(1)%d23(j,l)/zonal_glbl(j,l)
                   CH4rat  (1)%d23(j,l) = ch4vmr  (1)%d23(j,l)/zonal_glbl(j,l)
                   CH4ratnm(1)%d23(j,l) = ch4vmrnm(1)%d23(j,l)/zonal_glbl(j,l)
                END IF
             ENDDO
          ENDDO

       END IF

       DEALLOCATE( ppm, znl_lcl, zonal_glbl)
       
    ENDIF

    ! ---------------------------------------------------------------------
    ! Apply to O3 and CH4
    ! ---------------------------------------------------------------------

    DO j=j1,j2

          xlat = ybeg(region) + (j+0.5)*dyy

          IF (ABS(xlat) <= 30.0) THEN
             lm_min = LVLPRESS(region,  4500., 98400.)
             gnud = gnud_trop
          ELSEIF(ABS(xlat)<66.0 ) THEN
             lm_min = LVLPRESS(region,  9500., 98400.) 
             gnud = gnud_hilat
          ELSE
             lm_min = LVLPRESS(region, 12000., 98400.)
             gnud = gnud_pole
          ENDIF

          DO l=lm_min,lm(region)
             DO i=i1,i2
                !
                ! CMK: strictly spoken,  the procedure om multi PEs is somewhat different. 
                ! Now they are independently nudged towards the same ZA!
                ! In the long run, this will be the same! (hopefully) (FIXME : vraag MK wat betekend dat)
                !
#ifndef without_o3_nudging

                if(idate(3) .lt. 10) &
                     daily_scale = ((o3ratpm(region)%d23(j,l)*fraction_pm) + (o3rat(region)%d23(j,l)*fraction_cm))
                if(idate(3) .gt. 9 .and. idate(3) .lt. 21) &
                     daily_scale = o3rat(region)%d23(j,l)
                if(idate(3) .gt. 20) &
                     daily_scale = ((o3ratnm(region)%d23(j,l)*fraction_nm) + (o3rat(region)%d23(j,l)*fraction_cm)) 

                rmobs = rm(i,j,l,io3)*daily_scale  ! correct za
                rmold = rm(i,j,l,io3)                         ! uncorrected value
                z = (rmold+dtime*gnud*rmobs)/(1.+dtime*gnud)  ! new concentration
                rm(i,j,l,io3) = z

                rmobs1 = rm(i,j,l,io3s)*daily_scale  ! correct za
                rmold1 = rm(i,j,l,io3s)                         ! uncorrected value
                z1 = (rmold1+dtime*gnud*rmobs1)/(1.+dtime*gnud) ! new concentration
                rm(i,j,l,io3s) = z1
#endif
                
                if(idate(3) .lt. 10) &
                     daily_scale=((ch4ratpm(region)%d23(j,l)*fraction_pm) + (ch4rat(region)%d23(j,l)*fraction_cm))
                if(idate(3) .gt. 9 .and. idate(3) .lt. 21) &
                     daily_scale=ch4rat(region)%d23(j,l)
                if(idate(3) .gt. 20) &
                     daily_scale=((ch4ratnm(region)%d23(j,l)*fraction_nm)+ (ch4rat(region)%d23(j,l)*fraction_cm)) 
                
                rmobs2 = rm(i,j,l,ich4)*daily_scale             ! correct za
                rmold2 = rm(i,j,l,ich4)                         ! uncorrected value
                z2 = (rmold2+dtime*gnud*rmobs2)/(1.+dtime*gnud) ! new concentration
                rm(i,j,l,ich4) = z2

#ifdef with_budgets
                ! stratospheric budget....
                nzone=budg_dat(region)%nzong(i,j)    !global budget
                nzone_v=nzon_vg(l) 
#ifndef without_o3_nudging
                budstratg(nzone,nzone_v,io3) = budstratg(nzone,nzone_v,io3) + (z-rmold)/xmo3*1e3    ! mole
                budstratg(nzone,nzone_v,io3s) = budstratg(nzone,nzone_v,io3s) + (z1-rmold1)/xmo3*1e3 ! mole
#endif
                budstratg(nzone,nzone_v,ich4) = budstratg(nzone,nzone_v,ich4) + (z2-rmold2)/xmch4*1e3 ! mole

! #ifndef without_o3_nudging
!                 IF( io3 == 1 ) THEN
!                    sum_stratosphere(region) = sum_stratosphere(region) + z-rmold
!                 END IF
!                 IF ( io3s == 1 ) THEN
!                    sum_stratosphere(region) =  sum_stratosphere(region) + z1-rmold1
!                 END IF
! #endif
!                 IF ( ich4 == 1 ) THEN
!                    sum_stratosphere(region) =  sum_stratosphere(region) + z2-rmold2
!                 END IF
#endif /* BUDGETS */  

             END DO!i
          END DO !l
       END DO !j

    !-------------------------------------------
    ! set constraints for HNO3
    !-------------------------------------------
    IF (use_ODIN) then
       nlevels=3                ! limit nudging to 3 top levels
       pressure=(/550.0, 1122.0, 2820.0, 4365.0/)  ! hPa
       gnud_hno3 = (/gnud_hno3_5hpa, gnud_hno3_10hpa, gnud_hno3_28hpa, 0./)
    ELSE
       nlevels=1
       pressure=(/1000.0, 0., 0., 0./)  ! hPa
       gnud_hno3 = (/gnud_hno3_10hpa, 0., 0., 0./)
    ENDIF

    DO nudge_lev=1,nlevels
    
       lsave = 0   ! locate level closest to pressure[nudge_lev] hPa
       diffmin = 1e10
       DO l=1,lm(region)
          pr = 0.5*(at(l)+at(l+1)) + 0.5*(bt(l)+bt(l+1))*1e5    ! pressure centre box
          IF(ABS(pr-pressure(nudge_lev)) < diffmin) THEN
             diffmin = ABS(pr-pressure(nudge_lev))
             lsave = l
          ENDIF
       ENDDO

       IF (use_ODIN) then
          DO j=j1,j2
             !
             ! Smooth ratio between months to avoid spurious jumps in record
             ! "daily_scale" is the time-interpolated ODIN climatology of the [HNO3]/[O3] ratio of mixing ratios
             !
             if ( idate(3) .lt. 10 ) &
                  daily_scale = odin_hno3_o3_pm(region)%d23(j,nudge_lev) * fraction_pm + &
                                odin_hno3_o3   (region)%d23(j,nudge_lev) * fraction_cm
             
             if ( idate(3) .gt. 9 .and. idate(3) .lt. 21 ) &
                  daily_scale = odin_hno3_o3(region)%d23(j,nudge_lev)

             if ( idate(3) .gt. 20 ) &
                  daily_scale = odin_hno3_o3_nm(region)%d23(j,nudge_lev) * fraction_nm + &
                                odin_hno3_o3   (region)%d23(j,nudge_lev) * fraction_cm 

             DO i=i1,i2
                ! Include ratio of fscale factors, which convert from mass to mixing ratio
                tm5_o3_hno3_mix_ratio = ( rm(i,j,lsave,io3) * fscale(io3) ) / ( rm(i,j,lsave,ihno3) * fscale(ihno3) )
                nudgeFactor = ( 1.0 + dtime*gnud_hno3(NUDGE_LEV)*daily_scale*tm5_o3_hno3_mix_ratio ) / &
                              ( 1.0 + dtime*gnud_hno3(NUDGE_LEV) )
                rmold = rm(i,j,lsave,ihno3)
                z     = rmold * nudgeFactor
                rm(i,j,lsave,ihno3) = z
#ifdef with_budgets
                ! stratospheric budget....
                nzone=budg_dat(region)%nzong(i,j)    !global budget
                nzone_v=nzon_vg(lsave) 
                budstratg(nzone,nzone_v,ihno3)=budstratg(nzone,nzone_v,ihno3)+(z-rmold)/xmhno3*1e3 !mole
                ! IF(ihno3 ==1) sum_stratosphere(region) = sum_stratosphere(region) + z-rmold 
#endif
             ENDDO!i
          ENDDO!j
       ELSE
          DO j=j1,j2
          DO i=i1,i2

             rmold = rm(i,j,lsave,ihno3)
             rmobs = rm(i,j,lsave,io3) * odin_hno3_o3(region)%d23(j,nudge_lev)
             z = (rmold+dtime*gnud_hno3(NUDGE_LEV)*rmobs) / (1.+dtime*gnud_hno3(NUDGE_LEV)) !new concentration
             rm(i,j,lsave,ihno3) = z
#ifdef with_budgets
             ! stratospheric budget....
             nzone=budg_dat(region)%nzong(i,j)    !global budget
             nzone_v=nzon_vg(lsave) 
             budstratg(nzone,nzone_v,ihno3)=budstratg(nzone,nzone_v,ihno3)+(z-rmold)/xmhno3*1e3 !mole
             !          IF(ihno3 ==1) sum_stratosphere(region) = sum_stratosphere(region) + z-rmold 
#endif
          ENDDO!i
          ENDDO!j
       ENDIF
    ENDDO ! nudge levels
       
    ! This is not used (i.e. the global budget is not printed nor saved) so
    ! this is commented. NOTE also that this should be moved to
    ! chemistry/chemistry_done if ever used [PLS, 23-maart-2012]
!#ifdef with_budgets    
!    CALL PAR_REDUCE(sum_stratosphere(region), 'SUM', sum_stratosphere(region), status)
!    IF_NOTOK_RETURN(status=1)
!#endif

    
    !-------------------------------------------
    ! Set constraints for CO nudging
    !-------------------------------------------
    IF (USE_ODIN) THEN    
       
       nlevels=3                ! limit nudging to 3 top levels
       pressure=(/550.0, 1122.0, 2820.0, 4365.0/)  ! hPa
       gnud_co = (/gnud_hno3_5hpa, gnud_hno3_10hpa, gnud_hno3_28hpa, gnud_hno3_28hpa/)

       DO nudge_lev=1,nlevels         
    
          lsave = 0   ! locate level closest to pressure level that is nudged
          diffmin = 1e10
          DO l=1,lm(region)
             pr = 0.5*(at(l)+at(l+1)) + 0.5*(bt(l)+bt(l+1))*1e5    ! pressure centre box
             IF(ABS(pr-1000.0) < diffmin) THEN
                diffmin = ABS(pr-pressure(nudge_lev))
                lsave = l
             ENDIF
          ENDDO

          DO j=j1,j2
            !
            ! Smooth ratio between months to avoid spurious jumps in record
            ! "daily_scale" is the time-interpolated ODIN climatology of the [CO]/[O3] ratio of mixing ratios
            !
            if ( idate(3) .lt. 10 ) &
                 daily_scale = odin_co_o3_pm(region)%d23(j,nudge_lev) * fraction_pm + &
                               odin_co_o3   (region)%d23(j,nudge_lev) * fraction_cm

            if ( idate(3) .gt. 9 .and. idate(3) .lt. 21 ) &
                 daily_scale = odin_co_o3(region)%d23(j,nudge_lev)

            if ( idate(3) .gt. 20 ) &
                 daily_scale = odin_co_o3_nm(region)%d23(j,nudge_lev) * fraction_nm + &
                               odin_co_o3   (region)%d23(j,nudge_lev) * fraction_cm 

            DO i=i1,i2
              ! Include ratio of fscale factors, which convert from mass to mixing ratio
              tm5_o3_co_mix_ratio = ( rm(i,j,lsave,io3) * fscale(io3) ) / ( rm(i,j,lsave,ico) * fscale(ico) )
              nudgeFactor = ( 1.0 + dtime*gnud_co(NUDGE_LEV)*daily_scale*tm5_o3_co_mix_ratio ) / &
                            ( 1.0 + dtime*gnud_co(NUDGE_LEV) )
              rmold = rm(i,j,lsave,ico)
              z     = rmold * nudgeFactor
              rm(i,j,lsave,ico) = z

#ifdef with_budgets
              ! stratospheric budget....
              nzone=budg_dat(region)%nzong(i,j)    !global budget
              nzone_v=nzon_vg(lsave) 
              budstratg(nzone,nzone_v,ico)=budstratg(nzone,nzone_v,ico)+(z-rmold)/xmco*1e3 !mole
              ! IF(ico ==1) sum_stratosphere(region) = sum_stratosphere(region) + z-rmold 
#endif          
            ENDDO!i
          ENDDO!j
       
       ENDDO ! nudge levels

    ENDIF
    
    NULLIFY(rm)
    NULLIFY(m)

    ! end timing:
    call GO_Timer_End( itim_appl, status )
    IF_NOTOK_RETURN(status=1)

    CALL goLabel()
    status = 0

  END SUBROUTINE BOUNDARY_APPLY
  !EOC
  
END MODULE BOUNDARY