#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-2099 (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). ! ! For all CMIP6 scenarios, see Gidden et al., 2018; https://doi.org/10.5194/gmd-2018-266 ! !\\ !\\ ! !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 LOGICAL :: L1PCTCO2, LA4xCO2 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 !------------------------------ ! SSP scenario name !------------------------------ character(len=14) :: ssp_name, ssp_name_ch4 logical :: LSSP370_LowCH4 ! 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_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.1pct', L1PCTCO2, status, default=.FALSE. ) IF_ERROR_RETURN(status=1) CALL ReadRc( rcF, 'input.co2.abrupt-4x', LA4xCO2, status, default=.FALSE. ) IF_ERROR_RETURN(status=1) IF (L1PCTCO2 .AND. LA4xCO2) THEN write (gol,'("ERROR - L1PCTCO2 and LA4xCO2 both set")') ; call goErr status=1; TRACEBACK; return ELSE IF (L1PCTCO2 .OR. LA4xCO2) THEN !Use co2_fixyear to set reference year to 1850 co2_fixyear=.TRUE. co2_year=1850 IF (L1PCTCO2) THEN write(gol,*) 'DECK 1pctCO2 simulation:'; call goPr write(gol,*) '...A stepwise change of 1% is applied from one year to the next'; call goPr write(gol,*) '...In the first year a 5% increase is applied compared to the 1850 value'; call goPr write(gol,*) '...The simulation should start in 1850.'; call goPr ENDIF IF (LA4xCO2) THEN write(gol,*) 'DECK abrupt-4xCO2 simulation: CO2 quadrupled compared to 1850'; call goPr ENDIF ELSE 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 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 IF (LCMIP6_O3 .OR. LCMIP6_CH4 .OR. LCMIP6_CO2) THEN call ReadRc( rcF, 'input.CMIP6.SSP', ssp_name, status, default='' ) IF_ERROR_RETURN(status=1) write(gol, '("SSP CMIP6 future scenario for boundary conditions: ",a)') trim(ssp_name); call goPr IF (LCMIP6_CH4) THEN call ReadRc( rcF, 'input.CMIP6.SSP370_LowNTCF.ch4', LSSP370_LowCH4, status, default=.FALSE. ) IF_ERROR_RETURN(status=1) IF (.not.LSSP370_LowCH4) THEN ssp_name_ch4 = ssp_name ELSE ssp_name_ch4 = "SSP3-LowNTCF" write(gol, '("... but for CH4 use instead: ",a)') trim(ssp_name_ch4); call goPr ENDIF ENDIF 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, 137) ! 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),'o3_pi/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_clim_1850.nc' else ! O3 data are only provided up to 2099 if (.not.o3_fixyear) then target_year = MIN(2099,MAX(idate(1),1850)) else target_year = MIN(2099,MAX(o3_year,1850)) endif if ( target_year <= 2014) then ! 2014 is the last year of the historical period. ! The last entry of the 2014 file ! contains the field for January 2015, ! which is taken from the SSP370 scenario. ! This is the only scenario implemented in TM5 anyway. write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'o3_histo/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_',target_year,'.nc' else ! Consistent with EC-Earth: ! - use same names ! - for SSP1-1.9, use ozone from SSP1-2.6 ! - for Extended scenarios, repeat 2100 value of the corresponding non-extended scenario (NOT TESTED) ! - for Tier-2 scenarios, not available ("SSP4-3.4", "SSP4-6.0", "SSP5-3.4-OS" and its extended "SSP5-3.4-OS-Ext") ! On top of it, use SSP370 ozone also for SSP370-LowNTCF select case (trim(ssp_name)) case ("SSP1-1.9", "SSP1-2.6", "SSP1-2.6-Ext") write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'o3_scenarios/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_ssp126_',target_year,'.nc' case ("SSP2-4.5") write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'o3_scenarios/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_ssp245_',target_year,'.nc' case ("SSP3-7.0") write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'o3_scenarios/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_ssp370_',target_year,'.nc' case ("SSP5-8.5","SSP5-8.5-Ext") write(o3vmr_fname,'(a,a,i4,a)') TRIM(o3vmr_dirname),'o3_scenarios/vmro3_input4MIPs_ozone_CMIP6_UReading-CCMI_ssp585_',target_year,'.nc' case default write (gol,'("Ozone fields not implemented for scenario",1x,a)') trim(ssp_name) ; call goErr status=1; TRACEBACK; return end select endif endif write(gol,'(" Boundary - reading O3 file: ",a)') trim(o3vmr_fname); call goPr 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,".nc4")')TRIM(o3vmr_dirname),valid_year_pm,inp_levs WRITE(o3vmr_fname,' (a,"msr2_o3vmr",i4,"_",i2,".nc4")')TRIM(o3vmr_dirname),valid_year, inp_levs WRITE(o3vmr_fnamenm,'(a,"msr2_o3vmr",i4,"_",i2,".nc4")')TRIM(o3vmr_dirname),valid_year_nm,inp_levs ! ! January ! if(nmw .eq. 1) then CALL MDF_Open( TRIM(o3vmr_fnamepm), MDF_NETCDF, 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_NETCDF, 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_NETCDF, 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_NETCDF, 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_NETCDF, 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.nc4) allocate( field3d(1, nlat180, lm(1) )) field3d=0.0 o3vmr_fname=o3vmr_dirname CALL MDF_Open( TRIM(o3vmr_fname), MDF_NETCDF, MDF_READ, hid, status ) IF_NOTOK_RETURN(status=1) !FIXME CALL MDF_Inq_VarID( hid, "idate(2)", varid, status ) !FIXME IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( hid, varid, 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)-1 ! pick out right month in file WRITE(o3du_fname,'(a,"msr2_o3du",i4,"_",i2,".nc4")')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_NETCDF, 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(gol,'(" Boundary - reading CH4 file: ",a)') trim(ch4_hist_fname); call goPr 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) == 2015 .and. idate(2) == 12) .or. (idate(1) > 2015 ) ) ) .or. & ( ch4_fixyear .and. (ch4_year >= 2015) ) ) then select case (trim(ssp_name_ch4)) case ('SSP1-1.9') WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp119-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP1-2.6", "SSP1-2.6-Ext") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp126-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP2-4.5") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-MESSAGE-GLOBIOM-ssp245-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP3-7.0") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-AIM-ssp370-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP3-LowNTCF") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_AerChemMIP_UoM-AIM-ssp370-lowNTCF-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP4-3.4") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp434-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP4-6.0") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp460-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP5-3.4-OS", "SSP5-3.4-OS-Ext") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp534-over-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case ("SSP5-8.5","SSP5-8.5-Ext") WRITE(ch4_sce_fname,'(a,"mole-fraction-of-methane-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp585-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_ch4_dirname_strat) case default write (gol,'("CH4 global annual means not implemented for this scenario")') ; call goErr status=1; TRACEBACK; return end select allocate(ch4_sce(1,486)) ! years 2015-2500 write(gol,'(" Boundary - reading CH4 file: ",a)') trim(ch4_sce_fname); call goPr 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,486/) ) 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(2500,MAX(idate(1)-1,1850)) if(idate(2) .eq. 1) then target_year = MIN(2500,MAX(idate(1)-2,1850)) endif else target_year = MIN(2500,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(2500,MAX(idate(1)-1,1850)) else target_year = MIN(2500,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(2500,MAX(idate(1)-1,1850)) if(idate(2) .eq. 12) then target_year = MIN(2500,MAX(idate(1),1850)) endif else target_year = MIN(2500,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,".nc4")')TRIM(ch4vmr_dirname),inp_levs CALL MDF_Open( TRIM(ch4vmr_fname), MDF_NETCDF, 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.nc4', MDF_NETCDF, 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.nc4', MDF_NETCDF, 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.nc4', MDF_NETCDF, MDF_READ, hid, status ) IF_NOTOK_RETURN(status=1) !FIXME CALL MDF_Inq_VarID( hid, "idate(2)", varid, status ) !FIXME IF_NOTOK_RETURN(status=1) CALL MDF_Get_Var( hid, varid, 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(2500,MAX(idate(1),1850)) else target_year = MIN(2500,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(gol,'(" Boundary - reading CO2 file: ",a)') trim(co2_hist_fname); call goPr 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) IF (L1PCTCO2) THEN ! 1% increase per year starting from the reference value for 1850. ! A stepwise increase is applied from year to year. ! This is consistent with prescribing CO2 mixing ratios as annual means. ! In the first year of the 1pctCO2 simulation ! the mixing ratio is increased by 1.01**0.5, i.e. by almost 0.5%. ! Consistent with the implementation in IFS, ! it is assumed that the 1pctCO2 simulation starts in 1850. co2_glob = co2_glob * EXP((idate(1)-1850+0.5)*LOG(1.01)) ELSE IF (LA4xCO2) THEN ! Constant 4x increase compared to the reference 1850 value. co2_glob = 4. * co2_glob ENDIF deallocate(co2_hist) else select case (trim(ssp_name)) case ('SSP1-1.9') WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp119-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP1-2.6", "SSP1-2.6-Ext") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-IMAGE-ssp126-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP2-4.5") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-MESSAGE-GLOBIOM-ssp245-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP3-7.0") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-AIM-ssp370-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) ! NOT USED IN CMIP6 case ("SSP3-LowNTCF") ! NOT USED IN CMIP6 WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_AerChemMIP_UoM-AIM-ssp370-lowNTCF-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP4-3.4") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp434-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP4-6.0") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-GCAM4-ssp460-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP5-3.4-OS", "SSP5-3.4-OS-Ext") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp534-over-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case ("SSP5-8.5","SSP5-8.5-Ext") WRITE(co2_sce_fname,'(a,"mole-fraction-of-carbon-dioxide-in-air_input4MIPs_GHGConcentrations_ScenarioMIP_UoM-REMIND-MAGPIE-ssp585-1-2-1_gr1-GMNHSH_2015-2500.nc")')TRIM(cmip6_co2_dirname) case default write (gol,'("CO2 global annual means not implemented for this scenario")') ; call goErr status=1; TRACEBACK; return end select allocate(co2_sce(1,486)) ! years 2015-2500 write(gol,'(" Boundary - reading CO2 file: ",a)') trim(co2_sce_fname); call goPr 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,486/) ) 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