#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