Added option to write surface pressure/temperature to netcdf output

src/netcdf_output_mod.f90

@@ -35,9 +35,10 @@ module netcdf_output_mod
                        xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass
   use outg_mod,  only: outheight,oroout,densityoutgrid,factor3d,volume,&
                        wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
-                       area,arean,volumen, orooutn
+                       area,arean,volumen, orooutn, p0out, t0out
   use par_mod,   only: dep_prec, sp, dp, maxspec, maxreceptor, nclassunc,&
-                       unitoutrecept,unitoutreceptppt, nxmax,unittmp
+                       unitoutrecept,unitoutreceptppt, nxmax,unittmp, &
+                       write_p0t0
   use com_mod,   only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
                        loutstep,loutaver,loutsample,outlon0,outlat0,&
                        numxgrid,numygrid,dxout,dyout,numzgrid, height, &
@@ -56,7 +57,7 @@ module netcdf_output_mod
                        itsplit, lsynctime, ctl, ifine, lagespectra, ipin, &
                        ioutputforeachrelease, iflux, mdomainfill, mquasilag, & 
                        nested_output, ipout, surf_only, linit_cond, &
-                       flexversion,mpi_mode,DRYBKDEP,WETBKDEP
+                       flexversion,mpi_mode,DRYBKDEP,WETBKDEP, ps, tt2
 
   use mean_mod
 
@@ -80,6 +81,7 @@ module netcdf_output_mod
   ! netcdf dimension and variable IDs for main and nested output grid
   integer, dimension(maxspec) :: specID,specIDppt, wdspecID,ddspecID
   integer, dimension(maxspec) :: specIDn,specIDnppt, wdspecIDn,ddspecIDn
+  integer :: psID, tt2ID
   integer                     :: timeID, timeIDn
   integer, dimension(6)       :: dimids, dimidsn
   integer, dimension(5)       :: depdimids, depdimidsn
@@ -388,6 +390,14 @@ subroutine writeheader_netcdf(lnest)
   if (write_area) call nf90_err(nf90_def_var(ncid, 'area', nf90_float, &
        &(/ lonDimID, latDimID /), areaID))
 
+    ! surfarce pressure / temperature
+    if (write_p0t0) then
+       call nf90_err(nf90_def_var(ncid, 'pressure', nf90_float, &
+            &(/ lonDimID, latDimID, timeDimID /), psID))
+       call nf90_err(nf90_def_var(ncid, 'temperature', nf90_float, &
+            &(/ lonDimID, latDimID, timeDimID /), tt2ID))
+    end if
+
 
   if (write_releases.eqv..true.) then
     ! release comment
@@ -760,6 +770,8 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
   ! real(sp)            :: gridtotal,gridsigmatotal
   ! real(sp)            :: wetgridtotal,wetgridsigmatotal
   ! real(sp)            :: drygridtotal,drygridsigmatotal
+    real :: ddx,ddy,p0h,t0h,p1,p2,p3,p4,rddx,rddy,xlon,xlat,ylat,xtn,ytn
+    integer :: i1,ixp,j1,jyp,j
 
   real, parameter     :: weightair=28.97
 
@@ -1040,6 +1052,87 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
 
   end do
 
+
+    if (write_p0t0) then
+       ! Loop over all output grid cells
+       !********************************
+
+       do jjy=0,numygrid-1
+          do iix=0,numxgrid-1
+             p0h=0.
+             t0h=0.
+
+             ! Take 100 samples of the topography in every grid cell
+             !******************************************************
+
+             do j1=1,10
+                ylat=outlat0+(real(jjy)+real(j1)/10.-0.05)*dyout
+                yl=(ylat-ylat0)/dy
+                do i1=1,10
+                   xlon=outlon0+(real(iix)+real(i1)/10.-0.05)*dxout
+                   xl=(xlon-xlon0)/dx
+
+                   ! Determine the nest we are in
+                   !*****************************
+
+                   ngrid=0
+                   do j=numbnests,1,-1
+                      if ((xl.gt.xln(j)+eps).and.(xl.lt.xrn(j)-eps).and. &
+                           (yl.gt.yln(j)+eps).and.(yl.lt.yrn(j)-eps)) then
+                         ngrid=j
+                         goto 43
+                      endif
+                   end do
+43                 continue
+
+                   ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
+                   !*****************************************************************************
+
+                   if (ngrid.gt.0) then
+                      xtn=(xl-xln(ngrid))*xresoln(ngrid)
+                      ytn=(yl-yln(ngrid))*yresoln(ngrid)
+                      ix=int(xtn)
+                      jy=int(ytn)
+                      ddy=ytn-real(jy)
+                      ddx=xtn-real(ix)
+                   else
+                      ix=int(xl)
+                      jy=int(yl)
+                      ddy=yl-real(jy)
+                      ddx=xl-real(ix)
+                   endif
+                   ixp=ix+1
+                   jyp=jy+1
+                   rddx=1.-ddx
+                   rddy=1.-ddy
+                   p1=rddx*rddy
+                   p2=ddx*rddy
+                   p3=rddx*ddy
+                   p4=ddx*ddy
+
+                   p0h=p0h+p1*ps(ix ,jy,1,memind(1)) &
+                        + p2*ps(ixp,jy,1,memind(1)) &
+                        + p3*ps(ix ,jyp,1,memind(1)) &
+                        + p4*ps(ixp,jyp,1,memind(1))
+                   t0h=t0h+p1*tt2(ix ,jy,1,memind(1)) &
+                        + p2*tt2(ixp,jy,1,memind(1)) &
+                        + p3*tt2(ix ,jyp,1,memind(1)) &
+                        + p4*tt2(ixp,jyp,1,memind(1))
+                end do
+             end do
+
+             ! Divide by the number of samples taken
+             !**************************************
+             p0out(iix,jjy)=p0h/100.
+             t0out(iix,jjy)=t0h/100.
+          end do
+       end do
+
+       call nf90_err(nf90_put_var(ncid, psID, p0out, (/ 1,1,tpointer /), (/ numxgrid,numygrid,1 /)))
+       call nf90_err(nf90_put_var(ncid, tt2ID, t0out,(/ 1,1,tpointer /), (/ numxgrid,numygrid,1 /)))
+    end if
+
+
   ! Close netCDF file
   !**************************
   call nf90_err(nf90_close(ncid))

src/outg_mod.f90

@@ -11,6 +11,8 @@ module outg_mod
   real,allocatable, dimension (:) :: outheighthalf
   real,allocatable, dimension (:,:) :: oroout
   real,allocatable, dimension (:,:) :: orooutn
+  real,allocatable, dimension (:,:) :: t0out
+  real,allocatable, dimension (:,:) :: p0out
   real,allocatable, dimension (:,:) :: area
   real,allocatable, dimension (:,:) :: arean
   real,allocatable, dimension (:,:,:) :: volume

src/par_mod.f90

@@ -32,6 +32,13 @@ module par_mod
 
   integer,parameter :: dep_prec=sp
 
+  !***********************************************************
+  ! Additional output of lowest level pressure and temperature
+  ! (for netcdf output)
+  !***********************************************************
+  logical,parameter :: write_p0t0 = .false.
+
+
   !****************************************************************
   ! Set to F to disable use of kernel for concentrations/deposition
   !****************************************************************

src/readoutgrid.f90

@@ -209,6 +209,13 @@ subroutine readoutgrid
   if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
   allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
   if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
+! ESO: surface pressure/temperature
+  if (write_p0t0) then
+     allocate(t0out(0:numxgrid-1,0:numygrid-1),stat=stat)
+     allocate(p0out(0:numxgrid-1,0:numygrid-1),stat=stat)
+  end if
+
+
   return