ecearth3/sources/tm5mp/proj/cbm4/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: MSR2 = 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

#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
  !
  ! !PRIVATE DATA MEMBERS:
  !
  ! Flags for O3, CO and CH4 stratospheric nudging
  LOGICAL            :: use_MSR
  LOGICAL            :: use_HALOE
  LOGICAL            :: use_ODIN   

  character(len=256) :: o3vmr_dirname
  character(len=256) :: o3du_dirname
  character(len=256) :: ch4vmr_dirname
  character(len=256) :: covmr_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)
    
  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
  !
  ! !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 global_data,   ONLY : rcfile, inputdir
    USE partools,      ONLY : isRoot, par_broadcast
    USE MDF,           ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close

    USE binas,         ONLY : p0
    USE GO,            ONLY : TrcFile, Init, Done, ReadRc
    USE Grid,          ONLY : FillGrid, Fill3D
    USE Grid,          ONLY : TLevelInfo
    USE Grid,          ONLY : TllGridInfo, Init, Done
    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


    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)     ::  LeviX, src_lev

    ! some arrays on 1x1 resolution
    REAL, ALLOCATABLE    ::  field3d_pm(:,:,:), field3d(:,:,:), field3d_nm(:,:,:), field3d_h(:,:,:)
    REAL, ALLOCATABLE    ::  sp_z(:,:)
    !
    !  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.climat.MSR', use_MSR, status, default=.FALSE. )
       IF_ERROR_RETURN(status=1)

       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

       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.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 Done( rcF, status )
       IF_NOTOK_RETURN(status=1)

       ! 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
       !--------------------------------------------------

       ! Define vertical grid of original met fields (MSR2 only on 60 levels for now)
       ecstring='ec60'

       IF (isRoot) THEN
          CALL Init(LeviX, ecstring, status)
          IF_NOTOK_RETURN(status=1) 
       ENDIF

       ! to read data on root on the iglbsfc grid
       !nlat180  = dgrid(iglbsfc)%jm_region  ! use get_distgrid to get that
       inp_levs = 60            ! lme

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

          ! read the meridional data       
          IF (use_MSR) THEN

             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)
             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))
             ! 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, 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
                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, 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
                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, 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
                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    )


          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

             CALL Fill3D( lli_z(1),       levi, 'n',       sp_z,       wrk3d_ML, &
                          lli_z(iglbsfc), leviX, 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 = lme

             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)
             
             valid_year = MIN(2035,MAX(idate(1)-1,1999))
             DO k=1,inp_levs
                field3d_pm(1,:,k)=field3d_h(1,:,inp_levs+1-k)*RCP6_CH4_STRAT(valid_year-1998)
             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)      

             valid_year = MIN(2035,MAX(idate(1),1999))
             DO k=1,inp_levs
                field3d(1,:,k)=field3d_h(1,:,inp_levs+1-k)*RCP6_CH4_STRAT(valid_year-1998)
             ENDDO
             !
             ! - Next month
             !             
             nmw=idate(2)+1
             valid_year = MIN(2035,MAX(idate(1),1999))

             if(idate(2) .eq. 12) then
                nmw=1
                valid_year = MIN(2035,MAX(idate(1)+1,1999))
             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)*RCP6_CH4_STRAT(valid_year-1998)
             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, field3d_pm, 'mass-aver', status    )            

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

             CALL Fill3D( lli_z(1),       levi,  'n',     sp_z,        wrk3dnm_ML, &
                          lli_z(iglbsfc), leviX, 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

       ! **
       IF (isRoot) THEN
          CALL Done( LeviX, status )
          IF_NOTOK_RETURN(status=1)
       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 )
    !
    ! !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 (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