! First include the set of model-wide compiler flags #include "tm5.inc" ! macro defs #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_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if !===================================================================== ! ! Central interface to link the assimilation code to the TM5 model ! ! The retrieval-assimilation will be performed on the root node ! The steps performed: ! - Gather all arrays needed by the assimilation ! - call AssimNO2( ntimestep ) ! - Scatter the arrays for the next TM5 time step ! ! The interface is designed such that the same fields are passed as in the ! NRT interface, with TM5 fields read from file ! ! The subroutine AssimNO2 performs a data assimilation time step. ! ! Henk Eskes, KNMI, March 2015 ! Interface to TM5-mp-CB05 !====================================================================== module assim_interf_mod ! Definition of structure to store TM5 data use MTM5Data, only : TTM5Data implicit none private public :: AssimNO2_interface, AssimNO2_interface_init, AssimNO2_interface_done ! public :: TM5Data ! This structure stores all info passed from TM5 to the assimilation type(TTM5Data) :: TM5Data ! scaling factor, the result of the assimilation real, dimension(:,:,:), allocatable :: noxscaling ! This array stores the 3D NOx mass field real, dimension(:,:,:), allocatable :: store_rm_nox ! and these arrays the slopes #ifdef slopes real, dimension(:,:,:), allocatable :: store_rxm_nox real, dimension(:,:,:), allocatable :: store_rym_nox real, dimension(:,:,:), allocatable :: store_rzm_nox #endif ! Arrays used to convert the surface pressure and orography from 3D to 2D real, dimension(:,:,:), allocatable :: sp_temp real, dimension(:,:,:), allocatable :: oro_temp real, dimension(:,:,:), allocatable :: gph_temp ! to store geopotential height character(len=*), parameter :: mname = 'assim_interf_mod' ! ODIN HNO3 / O3 climatology related logical :: doNudgeOdinHNO3 character(len=255) :: OdinHNO3File ! This array stores the 3D HNO3 mass field (to be updated by the HNO3/O3 nudging) real, dimension(:,:,:), allocatable :: store_rm_hno3 ! Used for debug printing logical, parameter :: showDebugPrints = .false. contains subroutine AssimNO2_interface_init ( rcFile ) !====================================================================== ! ! Initialise storage space for TM5 3D global arrays ! ! rcFile : name of the rc (settings) file ! Henk Eskes, KNMI, March 2015 !====================================================================== ! For printing information use GO, only : gol, goPr, goErr ! For reading keys from the .rc file use GO, only : TrcFile, Init, Done, ReadRc ! To identify root proces use partools, only : isRoot ! for dimensions use meteodata, only : global_lli, levi ! to read ODIN HNO3/O3 climatology use MOdinHNO3, only : ReadOdinHNO3 implicit none ! --- in/out character(*), intent(in) :: rcFile ! --- local character(len=*), parameter :: rname = mname//'/AssimNO2_interface_init' type(TrcFile) :: rcF integer :: imr, jmr, lmr, n integer :: status ! --- start code --- if ( isRoot ) then write(gol,'("Assim_interface: initialise interface")') ; call goPr end if ! entire region grid size n = 1 ! global region imr = global_lli(n)%nlon jmr = global_lli(n)%nlat lmr = levi%nlev ! Fill fixed-dim constants TM5Data%im = imr TM5Data%jm = jmr TM5Data%lm = lmr ! Allocate arrays on root, dummy arrays elsewhere if ( isRoot ) then allocate ( TM5Data%lon(imr) ) allocate ( TM5Data%lat(jmr) ) allocate ( TM5Data%hyai(lmr+1) ) allocate ( TM5Data%hybi(lmr+1) ) allocate ( TM5Data%hyam(lmr) ) allocate ( TM5Data%hybm(lmr) ) allocate ( TM5Data%ps(imr,jmr) ) allocate ( TM5Data%oro(imr,jmr) ) allocate ( TM5Data%t(imr,jmr,lmr) ) allocate ( TM5Data%ltropo(imr,jmr) ) allocate ( TM5Data%no2(imr,jmr,lmr) ) allocate ( noxscaling(imr,jmr,lmr) ) allocate ( store_rm_nox(imr,jmr,lmr) ) #ifdef slopes allocate ( store_rxm_nox(imr,jmr,lmr) ) allocate ( store_rym_nox(imr,jmr,lmr) ) allocate ( store_rzm_nox(imr,jmr,lmr) ) #endif allocate ( sp_temp(imr,jmr,1) ) allocate ( oro_temp(imr,jmr,1) ) allocate ( gph_temp(imr,jmr,lmr+1) ) else allocate ( TM5Data%lon(1) ) allocate ( TM5Data%lat(1) ) allocate ( TM5Data%hyai(1) ) allocate ( TM5Data%hybi(1) ) allocate ( TM5Data%hyam(1) ) allocate ( TM5Data%hybm(1) ) allocate ( TM5Data%ps(1,1) ) allocate ( TM5Data%oro(1,1) ) allocate ( TM5Data%t(1,1,1) ) allocate ( TM5Data%ltropo(1,1) ) allocate ( TM5Data%no2(1,1,1) ) allocate ( noxscaling(1,1,1) ) allocate ( store_rm_nox(1,1,1) ) #ifdef slopes allocate ( store_rxm_nox(1,1,1) ) allocate ( store_rym_nox(1,1,1) ) allocate ( store_rzm_nox(1,1,1) ) #endif allocate ( sp_temp(1,1,1) ) allocate ( oro_temp(1,1,1) ) allocate ( gph_temp(1,1,1) ) end if if ( isRoot ) then ! open rcfile: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! Apply the alternative ODIN HNO3/O3 climatology ? call ReadRc( rcF, 'domino.doNudgeOdinHNO3', doNudgeOdinHNO3, status, default=.false. ) IF_ERROR_RETURN(status=1) if ( doNudgeOdinHNO3 ) then ! Read the .nc file name (path/filename) which contains the ODIN HNO3/O3 climatology call ReadRc( rcF, 'domino.odinclimatologyfile', OdinHNO3File, status, default='list' ) IF_ERROR_RETURN(status=1) end if ! close rcfile: call Done( rcF, status ) IF_NOTOK_RETURN(status=1) end if ! Read the ODIN climatology of HNO3 and O3 if ( doNudgeOdinHNO3 ) then if ( isRoot ) then write(gol,'("AssimNO2_interface: initialise ODIN HNO3 climatology")') ; call goPr ! read the ODIN climatology call ReadOdinHNO3 ( OdinHNO3File, jmr, lmr, status ) IF_NOTOK_RETURN(status=1) ! Allocate the memory for the hno3 mass to be adjusted allocate ( store_rm_hno3(imr,jmr,lmr) ) else ! processor is not "root" ! Allocate the memory for the hno3 mass to be adjusted allocate ( store_rm_hno3(1,1,1) ) end if end if end subroutine AssimNO2_interface_init subroutine AssimNO2_interface_done ( ) !====================================================================== ! ! Release storage space for TM5 3D global arrays ! Henk Eskes, KNMI, March 2015 !====================================================================== ! For printing information use GO, only : gol, goPr, goErr ! To identify root proces use partools, only : isRoot ! Odin-HNO3 climatology use MOdinHNO3, only : DoneOdinHNO3 implicit none character(len=*), parameter :: rname = mname//'/AssimNO2_interface_done' ! --- start code --- if ( isRoot ) then write(gol,'("Assim_interface: deallocate arrays")') ; call goPr end if if ( allocated(TM5Data%lon) ) deallocate ( TM5Data%lon ) if ( allocated(TM5Data%lat) ) deallocate ( TM5Data%lat ) if ( allocated(TM5Data%hyai) ) deallocate ( TM5Data%hyai ) if ( allocated(TM5Data%hybi) ) deallocate ( TM5Data%hybi ) if ( allocated(TM5Data%hyam) ) deallocate ( TM5Data%hyam ) if ( allocated(TM5Data%hybm) ) deallocate ( TM5Data%hybm ) if ( allocated(TM5Data%ps) ) deallocate ( TM5Data%ps ) if ( allocated(TM5Data%oro) ) deallocate ( TM5Data%oro ) if ( allocated(TM5Data%t) ) deallocate ( TM5Data%t ) if ( allocated(TM5Data%ltropo) ) deallocate ( TM5Data%ltropo ) if ( allocated(TM5Data%no2) ) deallocate ( TM5Data%no2 ) if ( allocated(noxscaling) ) deallocate ( noxscaling ) if ( allocated(store_rm_nox) ) deallocate ( store_rm_nox ) #ifdef slopes if ( allocated(store_rxm_nox) ) deallocate ( store_rxm_nox ) if ( allocated(store_rym_nox) ) deallocate ( store_rym_nox ) if ( allocated(store_rzm_nox) ) deallocate ( store_rzm_nox ) #endif if ( allocated(sp_temp) ) deallocate ( sp_temp ) if ( allocated(oro_temp) ) deallocate ( oro_temp ) if ( allocated(gph_temp) ) deallocate ( gph_temp ) ! ODIN climatology of HNO3 and O3 if ( doNudgeOdinHNO3 ) then ! Remove ODIN-HNO3 storage call DoneOdinHNO3 ( ) ! De-allocate storage space of the gathered 3D HNO3 field if ( allocated(store_rm_hno3) ) deallocate ( store_rm_hno3 ) end if end subroutine AssimNO2_interface_done subroutine AssimNO2_interface( rcFile, ntimestep, status ) !====================================================================== ! ! Gather - Assimilate - Scatter ! Henk Eskes, KNMI, March 2015 !====================================================================== ! Observation vector dimensions: use tm5_distgrid, only : dgrid, Get_DistGrid, Gather, Scatter ! NOx tracer indices use chem_param, only : fscale, ino, ino2, ino3, in2o5, ihno4, inox use chem_param, only : ihno3, io3 ! Tracer concentrations: transported and chemistry-only use tracer_data, only : mass_dat, chem_dat ! Meteo fields use meteodata, only : temper_dat, oro_dat, sp_dat, m_dat, gph_dat ! grid/level info use meteodata, only : global_lli, levi ! contains "dxyp" use global_data, only : region_dat ! For printing information use GO, only : gol, goPr, goErr ! To identify root proces use partools, only : isRoot ! For computing "dxyp" use binas, only : ae, grav ! earth radius (m), gravitation constant use dims, only : dx, dy, gtor, idate ! Routine to perform an NO2 retrieval + assimilation step use MAssimilation, only : assimNO2 ! to nudge to ODIN HNO3/O3 climatology use MOdinHNO3, only : NudgeOdinHNO3 implicit none real, dimension(:,:,:), allocatable :: tile_no2_field ! For ODIN climatology nudging of HNO3 real, dimension(:,:,:), allocatable :: tile_o3_field, tile_hno3_field real, dimension(:,:,:), allocatable :: mixratio_3d_hno3, mixratio_3d_o3 real, dimension(:,:), allocatable :: mixratio_zonal_hno3, mixratio_zonal_o3 character(len=*), intent(in) :: rcfile ! name of rc file integer, intent(out) :: status integer, intent(in) :: ntimestep ! --- local varables --- character(len=*), parameter :: rname = mname//'/AssimNO2_interface' integer :: imr, jmr, lmr integer :: i1, i2, j1, j2, i, j, l, n integer :: nObs real :: mix_no, mix_no2, mix_no3, mix_n2o5, mix_hno4, mix_nox_transport real :: mix_hno3, mix_o3 real :: noxratio, chemnox real :: dxx, dyy, lat ! --- start code --- ! entire region grid size n = 1 ! global region imr = global_lli(n)%nlon jmr = global_lli(n)%nlat lmr = levi%nlev if ( isRoot ) then ! Fill grid and time properties that are not distributed TM5Data%date(:) = idate(:) TM5Data%im = imr TM5Data%jm = jmr TM5Data%lm = lmr TM5Data%lon(:) = global_lli(n)%lon_deg(:) TM5Data%lat(:) = global_lli(n)%lat_deg(:) TM5Data%hyai(1:lmr+1) = levi%a(0:lmr) TM5Data%hybi(1:lmr+1) = levi%b(0:lmr) TM5Data%hyam(1:lmr) = levi%fa(1:lmr) TM5Data%hybm(1:lmr) = levi%fb(1:lmr) if ( showDebugPrints ) then print*,'date ',TM5Data%date print*,'im, jm, lm ',TM5Data%im,TM5Data%jm,TM5Data%lm print*,'lon range ',TM5Data%lon(1),TM5Data%lon(imr) print*,'lat range ',TM5Data%lat(1),TM5Data%lat(jmr) print*,'hyai range ',TM5Data%hyai(1),TM5Data%hyai(lmr+1) print*,'hybi range ',TM5Data%hybi(1),TM5Data%hybi(lmr+1) print*,'hyam range ',TM5Data%hyam(1),TM5Data%hyam(lmr) print*,'hybm range ',TM5Data%hybm(1),TM5Data%hybm(lmr) end if end if ! *** Gather all arrays needed *** if ( isRoot ) then write(gol,'("Assim_interface: Gather NO2, NOx and meteo fields")') ; call goPr end if ! First, determine updated NO2 after NOx transport (distributed on all arrays) ! All NOx species are scaled with the same factor call Get_DistGrid( dgrid(n), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate( tile_no2_field(i1:i2, j1:j2, lmr) ) do l = 1, lmr do j = j1, j2 do i = i1, i2 ! Get NOx after chemistry from sum of non-transported components ! "fscale(ispec) = xmair/ra(ispec)" mix_no = chem_dat(n)%rm(i,j,l,ino) * fscale(ino) / m_dat(n)%data(i,j,l) mix_no2 = chem_dat(n)%rm(i,j,l,ino2) * fscale(ino2) / m_dat(n)%data(i,j,l) mix_no3 = chem_dat(n)%rm(i,j,l,ino3) * fscale(ino3) / m_dat(n)%data(i,j,l) mix_n2o5 = chem_dat(n)%rm(i,j,l,in2o5) * fscale(in2o5) / m_dat(n)%data(i,j,l) mix_hno4 = chem_dat(n)%rm(i,j,l,ihno4) * fscale(ihno4) / m_dat(n)%data(i,j,l) ! NOx after chemistry, rescaled with molecule mass chemnox = mix_no + mix_no2 + mix_no3 + 2.0*mix_n2o5 + mix_hno4 ! NOx after time stepping, rescaled with molecule mass mix_nox_transport = mass_dat(n)%rm(i,j,l,inox) * fscale(inox) / m_dat(n)%data(i,j,l) ! Get ratio of transported and non-transported NOx if ( chemnox > 1.0e-15 ) then noxratio = mix_nox_transport/chemnox else noxratio = 1.0 end if ! Re-scale NO2 mixing ratio field to concentration after transport tile_no2_field(i,j,l) = noxratio*mix_no2 if ( showDebugPrints ) then if (l==1 .and. j==j1 .and. i==295 ) & print '(a,3i4,f8.4,e13.4)',' i,j,l,nox_ratio,nox = ',i,j,l,noxratio,mix_nox_transport end if end do end do end do if ( showDebugPrints ) then if ( isRoot ) then print '(a,4i4)','i1, i2, j1, j2 ',i1, i2, j1, j2 print '(a,2e10.3)','fscale no, no2 ',fscale(ino),fscale(ino2) print '(a,34e10.2)', 'Tile no2 (i2,j2,:) ', tile_no2_field(i2,j2,:) end if end if ! Store NO2 field tiles into global NO2 array (volume mixing ratio) call Gather( dgrid(n), tile_no2_field, TM5Data%no2, 0, status) IF_NOTOK_RETURN(status=1) deallocate( tile_no2_field ) ! *** Nudging to ODIN climatology, gather the zonal means of HNO3 and O3 *** if ( doNudgeOdinHNO3 ) then write(gol,'("Assim_interface: compute ODIN HNO3 nudging factors ")') ; call goPr n = 1 call Get_DistGrid( dgrid(n), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! allocate the tiles allocate( tile_hno3_field(i1:i2, j1:j2, lmr) ) allocate( tile_o3_field(i1:i2, j1:j2, lmr) ) do l = 1, lmr do j = j1, j2 do i = i1, i2 ! Get O3 and HNO3 mixing ratios ! "fscale(ispec) = xmair/ra(ispec)" tile_hno3_field(i,j,l) = mass_dat(n)%rm(i,j,l,ihno3) * fscale(ihno3) / m_dat(n)%data(i,j,l) tile_o3_field(i,j,l) = mass_dat(n)%rm(i,j,l,io3) * fscale(io3) / m_dat(n)%data(i,j,l) end do end do end do if ( showDebugPrints ) then if ( isRoot ) then print '(a,4i4)','i1, i2, j1, j2 ',i1, i2, j1, j2 print '(a,2e10.3)','fscale hno3, o3 ',fscale(ihno3),fscale(io3) print '(a,34e10.2)', 'Tile hno3 (i2,j2,:) ', tile_hno3_field(i2,j2,:) print '(a,34e10.2)', 'Tile o3 (i2,j2,:) ', tile_o3_field(i2,j2,:) end if end if ! Allocate the memory for the hno3 / o3 3D mixing ratio fields on root if ( isRoot ) then allocate ( mixratio_3d_hno3(imr,jmr,lmr) ) allocate ( mixratio_3d_o3(imr,jmr,lmr) ) else allocate ( mixratio_3d_hno3(1,1,1) ) allocate ( mixratio_3d_o3(1,1,1) ) end if ! Store O3, HNO3 field tiles into global arrays (volume mixing ratio) call Gather( dgrid(n), tile_hno3_field, mixratio_3d_hno3, 0, status) IF_NOTOK_RETURN(status=1) call Gather( dgrid(n), tile_o3_field, mixratio_3d_o3, 0, status) IF_NOTOK_RETURN(status=1) ! remove the tiles deallocate( tile_hno3_field ) deallocate( tile_o3_field ) ! Store HNO3 mass field tiles into global HNO3 array (kg) call Gather( dgrid(n), mass_dat(n)%rm(:,:,:,ihno3), store_rm_hno3, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if ( isRoot ) then ! Allocate zonal mean arrays allocate ( mixratio_zonal_hno3(jmr,lmr) ) allocate ( mixratio_zonal_o3(jmr,lmr) ) ! Determine 2D zonal mean mixing ratio fields of HNO3 and O3 do l = 1, lmr do j = 1, jmr mixratio_zonal_hno3(j,l) = 0.0 mixratio_zonal_o3(j,l) = 0.0 do i = 1, imr mixratio_zonal_hno3(j,l) = mixratio_zonal_hno3(j,l) + mixratio_3d_hno3(i,j,l) mixratio_zonal_o3(j,l) = mixratio_zonal_o3(j,l) + mixratio_3d_o3(i,j,l) end do mixratio_zonal_hno3(j,l) = mixratio_zonal_hno3(j,l) /imr mixratio_zonal_o3(j,l) = mixratio_zonal_o3(j,l) / imr end do end do if ( showDebugPrints ) then print '(a,3i4,i5,5i3,i6)', 'assim_interf: imr, jmr, lmr, idate, timestep = ', imr, jmr, lmr, idate, ntimestep print '(a,34e10.2)', 'assim_interf: mixratio_zonal_hno3 j 120 = ',mixratio_zonal_hno3(120,:) print '(a,34e10.2)', 'assim_interf: mixratio_zonal_o3 j 120 = ',mixratio_zonal_o3(120,:) print '(a,34e10.2)', 'assim_interf: store_rm_hno3 i 1 j 120 = ',store_rm_hno3(1,120,:) print '(a)', 'assim_interf: calling NudgeOdinHNO3' end if ! Determine scaled HNO3 field call NudgeOdinHNO3 ( imr, jmr, lmr, idate, ntimestep, mixratio_zonal_hno3, mixratio_zonal_o3, store_rm_hno3 ) ! Remove zonal mean arrays deallocate ( mixratio_zonal_hno3 ) deallocate ( mixratio_zonal_o3 ) end if ! isRoot ! release storage of the 3D mixing ratio fields deallocate ( mixratio_3d_hno3 ) deallocate ( mixratio_3d_o3 ) end if ! doNudgeOdinHNO3 ! *** Gather other TM5 fields needed *** ! Store NOx mass field tiles into global NOx array (kg) call Gather( dgrid(n), mass_dat(n)%rm(:,:,:,inox), store_rm_nox, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) ! Same for the slopes #ifdef slopes call Gather( dgrid(n), mass_dat(n)%rxm(:,:,:,inox), store_rxm_nox, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call Gather( dgrid(n), mass_dat(n)%rym(:,:,:,inox), store_rym_nox, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call Gather( dgrid(n), mass_dat(n)%rzm(:,:,:,inox), store_rzm_nox, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) #endif ! Surface pressure (Pa): sp(imr,jmr,1) call Gather( dgrid(n), sp_dat(n)%data, sp_temp, sp_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) ! Orography (m2/s2): oro(imr,jmr,1), divide by "g" for meters call Gather( dgrid(n), oro_dat(n)%data, oro_temp, oro_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) ! Temperature field (K): temper(imr,jmr,lmr) call Gather( dgrid(n), temper_dat(n)%data, TM5Data%t, temper_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) ! Geopotential height (m): gph(imr,jmr,lmr+1) call Gather( dgrid(n), gph_dat(n)%data, gph_temp, gph_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if ( isRoot ) then ! Store surface pressure TM5Data%ps(:,:) = sp_temp(:,:,1) ! Store model orography TM5Data%oro(:,:) = oro_temp(:,:,1)/grav ! Determine tropopause level and store do i = 1, imr do j = 1, jmr TM5Data%ltropo(i,j) = ltropo(lmr,TM5Data%hyai,TM5Data%hybi,TM5Data%t(i,j,:),gph_temp(i,j,:)) end do end do ! Initialise assimilation scaling factor noxscaling(:,:,:) = 1.0 ! *** Call the assimilation *** ! Assimilation is done only on the root processor if ( showDebugPrints ) then i = 350 j = 120 print '(a,i4,a,i4,a)', "Gathered data for cell (",i,", ",j,":, 1) " print '(a,10e9.2)', " NO2 ", TM5Data%no2(i, j:j+9, 1) print '(a,10e9.2)', " NOx ", store_rm_nox(i, j:j+9, 1) print '(a,10e9.2)', " NOx slope x ", store_rxm_nox(i, j:j+9, 1) print '(a,10e9.2)', " NOx slope y ", store_rym_nox(i, j:j+9, 1) print '(a,10e9.2)', " NOx slope z ", store_rzm_nox(i, j:j+9, 1) print '(a,10f9.1)', " temp ", TM5Data%t(i, j:j+9, 1) print '(a,10f9.1)', " oro ", TM5Data%oro(i, j:j+9) print '(a,180i3)' , " ltropo_lat ", TM5Data%ltropo(i, 1:180) print '(a,180i3)' , " ltropo_lon ", TM5Data%ltropo(1:360, j) print '(a,180f5.0)', " ptropo ", TM5Data%hyai(TM5Data%ltropo(i, 1:180))*0.01+TM5Data%hybi(TM5Data%ltropo(i, 1:180))*1013.0 print '(a,35f6.0)', " plevs ", TM5Data%hyai(:)*0.01+TM5Data%hybi(:)*1013.0 print '(a,10f9.0)', " sp ", TM5Data%ps(i, j:j+9)*0.01 end if write(gol,'("Assim_interface: Fields have been gathered, now calling assimNO2")') ; call goPr call assimNO2( rcFile, ntimestep, TM5Data, nObs, noxscaling, status ) IF_NOTOK_RETURN(status=1) if ( nObs > 0 ) then ! *** Finally, scatter the analysis NOx field *** write(gol,'("Assim_interface: Update NOx field, slopes, and scatter")') ; call goPr ! Apply assimilation scaling to NOx (only) ... store_rm_nox(:,:,:) = store_rm_nox(:,:,:) * noxscaling(:,:,:) ! Scale also slopes #ifdef slopes store_rxm_nox(:,:,:) = store_rxm_nox(:,:,:) * noxscaling(:,:,:) store_rym_nox(:,:,:) = store_rym_nox(:,:,:) * noxscaling(:,:,:) store_rzm_nox(:,:,:) = store_rzm_nox(:,:,:) * noxscaling(:,:,:) #endif end if ! nObs > 0 end if ! isRoot ! ... and scatter NOx back to the MPI tiles call Scatter( dgrid(n), mass_dat(n)%rm(:,:,:,inox), store_rm_nox, & mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) #ifdef slopes call Scatter( dgrid(n), mass_dat(n)%rxm(:,:,:,inox), store_rxm_nox, & mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call Scatter( dgrid(n), mass_dat(n)%rym(:,:,:,inox), store_rym_nox, & mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call Scatter( dgrid(n), mass_dat(n)%rzm(:,:,:,inox), store_rzm_nox, & mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) #endif ! *** Scatter result of the HNO3 nudging to ODIN climatology *** if ( doNudgeOdinHNO3 ) then write(gol,'("Assim_interface: now scattering updated HNO3 concentration (ODIN nudging) ")') ; call goPr call Scatter( dgrid(n), mass_dat(n)%rm(:,:,:,ihno3), store_rm_hno3, & mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) end if end subroutine AssimNO2_interface !----------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------- !BOP ! ! !FUNCTION: ltropo ! ! !DESCRIPTION: determine tropopause level from temperature gradient ! reference: WMO /Bram Bregman !\\ !\\ ! !INTERFACE: ! integer function ltropo(lmx,at,bt,tx,gphx) ! ! !INPUT PARAMETERS: ! integer,intent(in) :: lmx real,dimension(lmx+1), intent(in) :: at, bt real,dimension(lmx),intent(in) :: tx real,dimension(lmx+1),intent(in) :: gphx ! ! !REVISION HISTORY: ! programmed by fd 01-11-1996 ! changed to function fd 12-07-2002 ! 16 Jul 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC integer :: ltmin,ltmax,l real :: dz,dt ! ltropo is highest tropospheric level ! above is defined as stratosphere ! min tropopause level 450 hPa (ca. 8 km) ltmin=lvlpress(lmx,at,bt,45000.,98400.) ! max tropopause level 70 hPa (ca. 20 km) ltmax=lvlpress(lmx,at,bt,7000.,98400.) ltropo=ltmin do l=ltmin,ltmax dz=(gphx(l)-gphx(l-2))/2. dt=(tx(l)-tx(l-1))/dz if ( dt < -0.002) ltropo=l !wmo 85 criterium for stratosphere ! increase upper tropospheric level end do !l end function ltropo !EOC !--------------------------------------------------------------------------- ! TM5 ! !--------------------------------------------------------------------------- !BOP ! ! !FUNCTION: lvlpress ! ! !DESCRIPTION: lvlpress determines the index of the level with a pressure ! greater than press i.e. in height below press ! determines level independent of vertical resolution lm ! uses hybrid level coefficients to determine lvlpress !\\ !\\ ! !INTERFACE: ! integer function lvlpress(lm,at,bt,press,press0) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: lm real, dimension(lm+1), intent(in) :: at, bt real, intent(in) :: press real, intent(in) :: press0 ! ! !REVISION HISTORY: ! programmed by Peter van Velthoven, 6-november-2000 ! 16 Jul 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC real :: zpress0, zpress integer :: l,lvl if ( press0 == 0.0 ) then zpress0 = 98400. else zpress0 = press0 endif lvl = 1 ! l increases from bottom (l=1) to top (l=lm) do l = 1, lm zpress = (at(l)+at(l+1) + (bt(l)+bt(l+1))*zpress0)*0.5 if ( press < zpress ) then lvl = l endif end do lvlpress = lvl end function lvlpress !EOC end module assim_interf_mod