Changes in src/concoutput_nest_mpi.f90 [92fab65:f3054ea] in flexpart.git

File:

• src/concoutput_nest_mpi.f90 (modified) (9 diffs)

src/concoutput_nest_mpi.f90

@@ -1 @@
-! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
-! SPDX-License-Identifier: GPL-3.0-or-later
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
 
 subroutine concoutput_nest(itime,outnum)
@@ -39 +57 @@
   !                                                                            *
   !*****************************************************************************
-
 
   use unc_mod
@@ -56 +73 @@
   real :: sp_fact
   real :: outnum,densityoutrecept(maxreceptor),xl,yl
+! RLT
+  real :: densitydryrecept(maxreceptor)
+  real :: factor_dryrecept(maxreceptor)
 
   !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
@@ -82 +102 @@
   character :: adate*8,atime*6
   character(len=3) :: anspec
+  logical :: lexist
   integer :: mind 
 ! mind        eso:added to ensure identical results between 2&3-fields versions
@@ -158 +179 @@
         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
              rho(iix,jjy,kzz-1,mind)*dz2)/dz
+! RLT
+        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
+             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
       end do
     end do
@@ -169 +193 @@
     !densityoutrecept(i)=rho(iix,jjy,1,2)
     densityoutrecept(i)=rho(iix,jjy,1,mind)
+! RLT
+    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
   end do
 
+! RLT
+! conversion factor for output relative to dry air
+  factor_drygrid=densityoutgrid/densitydrygrid
+  factor_dryrecept=densityoutrecept/densitydryrecept
 
   ! Output is different for forward and backward simulations
@@ -193 +223 @@
 
   write(anspec,'(i3.3)') ks
+ 
+  if (DRYBKDEP.or.WETBKDEP) then !scavdep output
+      if (DRYBKDEP) &
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_nest_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      if (WETBKDEP) &
+      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_nest_'//adate// &
+           atime//'_'//anspec,form='unformatted')
+      write(unitoutgrid) itime
+  else
   if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     if (ldirect.eq.1) then
@@ -205 +245 @@
      write(unitoutgrid) itime
    endif
+  endif
 
   if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
@@ -535 +576 @@
   end do
 
-
-
-  ! Reinitialization of grid
-  !*************************
-
-  do ks=1,nspec
-    do kp=1,maxpointspec_act
-      do i=1,numreceptor
-        creceptor(i,ks)=0.
-      end do
-      do jy=0,numygridn-1
-        do ix=0,numxgridn-1
-          do l=1,nclassunc
-            do nage=1,nageclass
-              do kz=1,numzgrid
-                griduncn(ix,jy,kz,ks,kp,l,nage)=0.
-              end do
-            end do
-          end do
-        end do
+! RLT Aug 2017
+! Write out conversion factor for dry air
+  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
+  if (lexist) then
+    ! open and append
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='old',action='write',access='append')
+  else
+    ! create new
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='new',action='write')
+  endif
+  sp_count_i=0
+  sp_count_r=0
+  sp_fact=-1.
+  sp_zer=.true.
+  do kz=1,numzgrid
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
+          if (sp_zer.eqv..true.) then ! first value not equal to one
+            sp_count_i=sp_count_i+1
+            sparse_dump_i(sp_count_i)= &
+                  ix+jy*numxgridn+kz*numxgridn*numygridn
+            sp_zer=.false.
+            sp_fact=sp_fact*(-1.)
+          endif
+          sp_count_r=sp_count_r+1
+          sparse_dump_r(sp_count_r)= &
+               sp_fact*factor_drygrid(ix,jy,kz)
+        else ! factor is one
+          sp_zer=.true.
+        endif
       end do
     end do
   end do
-
+  write(unitoutfactor) sp_count_i
+  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
+  write(unitoutfactor) sp_count_r
+  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
+  close(unitoutfactor)
+
+  ! Reinitialization of grid
+  !*************************
+
+  ! do ks=1,nspec
+  !   do kp=1,maxpointspec_act
+  !     do i=1,numreceptor
+  !       creceptor(i,ks)=0.
+  !     end do
+  !     do jy=0,numygridn-1
+  !       do ix=0,numxgridn-1
+  !         do l=1,nclassunc
+  !           do nage=1,nageclass
+  !             do kz=1,numzgrid
+  !               griduncn(ix,jy,kz,ks,kp,l,nage)=0.
+  !             end do
+  !           end do
+  !         end do
+  !       end do
+  !     end do
+  !   end do
+  ! end do
+  creceptor(:,:)=0.
+  griduncn(:,:,:,:,:,:,:)=0.
+  
   if (mp_measure_time) call mpif_mtime('iotime',1)
   ! if (mp_measure_time) then