Index: /trunk/src/FLEXPART.f90
===================================================================
--- /trunk/src/FLEXPART.f90 (revision 19)
+++ /trunk/src/FLEXPART.f90 (revision 20)
@@ -49,4 +49,6 @@
integer :: i,j,ix,jy,inest
integer :: idummy = -320
+ character(len=256) :: inline_options !pathfile, flexversion, arg2
+
! Generate a large number of random numbers
@@ -58,25 +60,106 @@
call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
+ !
+ flexversion='Version 9.1.8 (2013-12-08)'
+ !verbosity=0
+ ! Read the pathnames where input/output files are stored
+ !*******************************************************
+
+ inline_options='none'
+ select case (iargc())
+ case (2)
+ call getarg(1,arg1)
+ pathfile=arg1
+ call getarg(2,arg2)
+ inline_options=arg2
+ case (1)
+ call getarg(1,arg1)
+ pathfile=arg1
+ verbosity=0
+ if (arg1(1:1).eq.'-') then
+ write(pathfile,'(a11)') './pathnames'
+ inline_options=arg1
+ endif
+ case (0)
+ write(pathfile,'(a11)') './pathnames'
+ verbosity=0
+ end select
+
+ if (inline_options(1:1).eq.'-') then
+ print*, 'inline options=', inline_options
+ if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
+ print*, 'verbose mode 1: additional information will be displayed'
+ verbosity=1
+ endif
+ if (trim(inline_options).eq.'-v2') then
+ print*, 'verbose mode 2: additional information will be displayed'
+ verbosity=2
+ endif
+ if (trim(inline_options).eq.'-i') then
+ print*, 'info mode: will provide run specific information and stop'
+ verbosity=1
+ info_flag=1
+ endif
+ if (trim(inline_options).eq.'-i2') then
+ print*, 'info mode: will provide run specific information and stop'
+ verbosity=2
+ info_flag=1
+ endif
+ endif
+
+
! Print the GPL License statement
!*******************************************************
- print*,'Welcome to FLEXPART Version 9.0'
+ print*,'Welcome to FLEXPART', trim(flexversion)
print*,'FLEXPART is free software released under the GNU Genera'// &
'l Public License.'
-
- ! Read the pathnames where input/output files are stored
- !*******************************************************
-
- call readpaths
+
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call readpaths'
+ endif
+ call readpaths(pathfile)
+
+
+ if (verbosity.gt.1) then !show clock info
+ !print*,'length(4)',length(4)
+ !count=0,count_rate=1000
+ CALL SYSTEM_CLOCK(count_clock0, count_rate, count_max)
+ !WRITE(*,*) 'SYSTEM_CLOCK',count, count_rate, count_max
+ !WRITE(*,*) 'SYSTEM_CLOCK, count_clock0', count_clock0
+ !WRITE(*,*) 'SYSTEM_CLOCK, count_rate', count_rate
+ !WRITE(*,*) 'SYSTEM_CLOCK, count_max', count_max
+ endif
+
! Read the user specifications for the current model run
!*******************************************************
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call readcommand'
+ endif
call readcommand
+ if (verbosity.gt.0) then
+ WRITE(*,*) ' ldirect=', ldirect
+ WRITE(*,*) ' ibdate,ibtime=',ibdate,ibtime
+ WRITE(*,*) ' iedate,ietime=', iedate,ietime
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+ WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+ endif
+ endif
! Read the age classes to be used
!********************************
-
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call readageclasses'
+ endif
call readageclasses
+
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+ WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+ endif
+
@@ -84,12 +167,28 @@
!******************************************************************
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call readavailable'
+ endif
call readavailable
-
! Read the model grid specifications,
! both for the mother domain and eventual nests
!**********************************************
+
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call gridcheck'
+ endif
call gridcheck
+
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+ WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+ endif
+
+
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call gridcheck_nests'
+ endif
call gridcheck_nests
@@ -98,13 +197,24 @@
!************************************
+ if (verbosity.gt.0) then
+ WRITE(*,*) 'call readoutgrid'
+ endif
+
call readoutgrid
- if (nested_output.eq.1) call readoutgrid_nest
-
+
+ if (nested_output.eq.1) then
+ call readoutgrid_nest
+ if (verbosity.gt.0) then
+ WRITE(*,*) '# readoutgrid_nest'
+ endif
+ endif
! Read the receptor points for which extra concentrations are to be calculated
!*****************************************************************************
+ if (verbosity.eq.1) then
+ print*,'call readreceptors'
+ endif
call readreceptors
-
! Read the physico-chemical species property table
@@ -117,4 +227,7 @@
!***************************
+ if (verbosity.gt.0) then
+ print*,'call readlanduse'
+ endif
call readlanduse
@@ -123,4 +236,7 @@
!********************************************************************
+ if (verbosity.gt.0) then
+ print*,'call assignland'
+ endif
call assignland
@@ -130,16 +246,25 @@
!**********************************************
+ if (verbosity.gt.0) then
+ print*,'call readreleases'
+ endif
call readreleases
+
! Read and compute surface resistances to dry deposition of gases
!****************************************************************
+ if (verbosity.gt.0) then
+ print*,'call readdepo'
+ endif
call readdepo
-
! Convert the release point coordinates from geografical to grid coordinates
!***************************************************************************
- call coordtrafo
+ call coordtrafo
+ if (verbosity.gt.0) then
+ print*,'call coordtrafo'
+ endif
@@ -147,4 +272,7 @@
!*****************************************
+ if (verbosity.gt.0) then
+ print*,'Initialize all particles to non-existent'
+ endif
do j=1,maxpart
itra1(j)=-999999999
@@ -155,6 +283,12 @@
if (ipin.eq.1) then
+ if (verbosity.gt.0) then
+ print*,'call readpartpositions'
+ endif
call readpartpositions
else
+ if (verbosity.gt.0) then
+ print*,'numpart=0, numparticlecount=0'
+ endif
numpart=0
numparticlecount=0
@@ -166,4 +300,8 @@
!***************************************************************
+
+ if (verbosity.gt.0) then
+ print*,'call outgrid_init'
+ endif
call outgrid_init
if (nested_output.eq.1) call outgrid_init_nest
@@ -173,6 +311,10 @@
!******************
- if (OHREA.eqv..TRUE.) &
+ if (OHREA.eqv..TRUE.) then
+ if (verbosity.gt.0) then
+ print*,'call readOHfield'
+ endif
call readOHfield
+ endif
! Write basic information on the simulation to a file "header"
@@ -180,7 +322,25 @@
!******************************************************************
+
+ if (verbosity.gt.0) then
+ print*,'call writeheader'
+ endif
+
call writeheader
- if (nested_output.eq.1) call writeheader_nest
- open(unitdates,file=path(2)(1:length(2))//'dates')
+ ! FLEXPART 9.2 ticket ?? write header in ASCII format
+ call writeheader_txt
+ !if (nested_output.eq.1) call writeheader_nest
+ if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
+
+ if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
+ if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
+
+
+
+ !open(unitdates,file=path(2)(1:length(2))//'dates')
+
+ if (verbosity.gt.0) then
+ print*,'call openreceptors'
+ endif
call openreceptors
if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
@@ -190,4 +350,7 @@
!***************************************************************************
+ if (verbosity.gt.0) then
+ print*,'discretize release times'
+ endif
do i=1,numpoint
ireleasestart(i)=nint(real(ireleasestart(i))/ &
@@ -200,4 +363,8 @@
! Initialize cloud-base mass fluxes for the convection scheme
!************************************************************
+
+ if (verbosity.gt.0) then
+ print*,'Initialize cloud-base mass fluxes for the convection scheme'
+ endif
do jy=0,nymin1
@@ -218,4 +385,16 @@
!********************************
+ if (verbosity.gt.0) then
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+ WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+ endif
+ if (info_flag.eq.1) then
+ print*, 'info only mode (stop)'
+ stop
+ endif
+ print*,'call timemanager'
+ endif
+
call timemanager
Index: /trunk/src/com_mod.f90
===================================================================
--- /trunk/src/com_mod.f90 (revision 19)
+++ /trunk/src/com_mod.f90 (revision 20)
@@ -7,5 +7,5 @@
! June 1996 *
! *
-! Last update: 9 August 2000 *
+! Last update:15 August 2013 IP *
! *
!*******************************************************************************
@@ -25,8 +25,11 @@
character :: path(numpath+2*maxnests)*120
integer :: length(numpath+2*maxnests)
-
+ character(len=256) :: pathfile, flexversion, arg1, arg2
+
! path path names needed for trajectory model
! length length of path names needed for trajectory model
-
+ ! pathfile file where pathnames are stored
+ ! flexversion version of flexpart
+ ! arg input arguments from launch at command line
!********************************************************
@@ -65,5 +68,5 @@
integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
integer :: mquasilag,nested_output,ind_source,ind_receptor
- integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond
+ integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
logical :: turbswitch
@@ -91,4 +94,6 @@
! mquasilag 0: normal run
! 1: Particle position output is produced in a condensed format and particles are numbered
+ ! surf_only switch output in grid_time files for surface only or full vertical resolution
+ ! 0=no (full vertical resolution), 1=yes (surface only)
! nested_output: 0 no, 1 yes
! turbswitch determines how the Markov chain is formulated
@@ -139,4 +144,6 @@
real :: decay(maxspec)
real :: weta(maxspec),wetb(maxspec)
+! NIK: 31.01.2013- parameters for in-cloud scavening
+ real :: weta_in(maxspec), wetb_in(maxspec), wetc_in(maxspec), wetd_in(maxspec)
real :: reldiff(maxspec),henry(maxspec),f0(maxspec)
real :: density(maxspec),dquer(maxspec),dsigma(maxspec)
@@ -172,5 +179,7 @@
! WET DEPOSITION
- ! weta, wetb parameters for determining wet scavenging coefficients
+ ! weta, wetb parameters for determining below-cloud wet scavenging coefficients
+ ! weta_in, wetb_in parameters for determining in-cloud wet scavenging coefficients
+ ! wetc_in, wetd_in parameters for determining in-cloud wet scavenging coefficients
! GAS DEPOSITION
@@ -331,6 +340,10 @@
real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,2)
real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,2)
+ !scavenging NIK, PS
integer(kind=1) :: clouds(0:nxmax-1,0:nymax-1,nzmax,2)
integer :: cloudsh(0:nxmax-1,0:nymax-1,2)
+ integer icloudbot(0:nxmax-1,0:nymax-1,2)
+ integer icloudthck(0:nxmax-1,0:nymax-1,2)
+
! uu,vv,ww [m/2] wind components in x,y and z direction
@@ -346,4 +359,8 @@
! rainout conv/lsp dominated 2/3
! washout conv/lsp dominated 4/5
+! PS 2013
+!c icloudbot (m) cloud bottom height
+!c icloudthck (m) cloud thickness
+
! pplev for the GFS version
@@ -662,5 +679,4 @@
-
!********************
! Random number field
@@ -671,4 +687,11 @@
! rannumb field of normally distributed random numbers
+ !********************
+ ! Verbosity, testing flags
+ !********************
+ integer :: verbosity=0
+ integer :: info_flag=0
+ INTEGER :: count_clock, count_clock0, count_rate, count_max
+
end module com_mod
Index: /trunk/src/concoutput.f90
===================================================================
--- /trunk/src/concoutput.f90 (revision 19)
+++ /trunk/src/concoutput.f90 (revision 20)
@@ -108,6 +108,7 @@
write(adate,'(i8.8)') jjjjmmdd
write(atime,'(i6.6)') ihmmss
+ open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
write(unitdates,'(a)') adate//atime
-
+ close(unitdates)
! For forward simulations, output fields have dimension MAXSPEC,
@@ -158,5 +159,5 @@
yl=outlat0+real(jy)*dyout
xl=(xl-xlon0)/dx
- yl=(yl-ylat0)/dx
+ yl=(yl-ylat0)/dy !v9.1.1
iix=max(min(nint(xl),nxmin1),0)
jjy=max(min(nint(yl),nymin1),0)
Index: /trunk/src/concoutput_surf.f90
===================================================================
--- /trunk/src/concoutput_surf.f90 (revision 20)
+++ /trunk/src/concoutput_surf.f90 (revision 20)
@@ -0,0 +1,744 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc, &
+ drygridtotalunc)
+ ! i i o o
+ ! o
+ !*****************************************************************************
+ ! *
+ ! Output of the concentration grid and the receptor concentrations. *
+ ! *
+ ! Author: A. Stohl *
+ ! *
+ ! 24 May 1995 *
+ ! *
+ ! 13 April 1999, Major update: if output size is smaller, dump output *
+ ! in sparse matrix format; additional output of *
+ ! uncertainty *
+ ! *
+ ! 05 April 2000, Major update: output of age classes; output for backward*
+ ! runs is time spent in grid cell times total mass of *
+ ! species. *
+ ! *
+ ! 17 February 2002, Appropriate dimensions for backward and forward runs *
+ ! are now specified in file par_mod *
+ ! *
+ ! June 2006, write grid in sparse matrix with a single write command *
+ ! in order to save disk space *
+ ! *
+ ! 2008 new sparse matrix format *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! outnum number of samples *
+ ! ncells number of cells with non-zero concentrations *
+ ! sparse .true. if in sparse matrix format, else .false. *
+ ! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
+ ! *
+ !*****************************************************************************
+
+ use unc_mod
+ use point_mod
+ use outg_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ real(kind=dp) :: jul
+ integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
+ integer :: sp_count_i,sp_count_r
+ real :: sp_fact
+ real :: outnum,densityoutrecept(maxreceptor),xl,yl
+
+ !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
+ ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
+ ! + maxageclass)
+ !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
+ ! + maxageclass)
+ !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass)
+ !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
+ ! + maxpointspec_act,maxageclass),
+ ! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass),
+ ! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass)
+ !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
+ !real sparse_dump_r(numxgrid*numygrid*numzgrid)
+ !integer sparse_dump_i(numxgrid*numygrid*numzgrid)
+
+ !real sparse_dump_u(numxgrid*numygrid*numzgrid)
+ real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
+ real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
+ real :: drygridtotal,drygridsigmatotal,drygridtotalunc
+ real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
+ real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+ real,parameter :: weightair=28.97
+ logical :: sp_zer
+ character :: adate*8,atime*6
+ character(len=3) :: anspec
+
+
+ if (verbosity.eq.1) then
+ print*,'inside concoutput_surf '
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+
+ ! Determine current calendar date, needed for the file name
+ !**********************************************************
+
+ jul=bdate+real(itime,kind=dp)/86400._dp
+ call caldate(jul,jjjjmmdd,ihmmss)
+ write(adate,'(i8.8)') jjjjmmdd
+ write(atime,'(i6.6)') ihmmss
+ !write(unitdates,'(a)') adate//atime
+
+ open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
+ write(unitdates,'(a)') adate//atime
+ close(unitdates)
+
+ ! For forward simulations, output fields have dimension MAXSPEC,
+ ! for backward simulations, output fields have dimension MAXPOINT.
+ ! Thus, make loops either about nspec, or about numpoint
+ !*****************************************************************
+
+
+ if (ldirect.eq.1) then
+ do ks=1,nspec
+ do kp=1,maxpointspec_act
+ tot_mu(ks,kp)=1
+ end do
+ end do
+ else
+ do ks=1,nspec
+ do kp=1,maxpointspec_act
+ tot_mu(ks,kp)=xmass(kp,ks)
+ end do
+ end do
+ endif
+
+
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf 2'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+
+ !*******************************************************************
+ ! Compute air density: sufficiently accurate to take it
+ ! from coarse grid at some time
+ ! Determine center altitude of output layer, and interpolate density
+ ! data to that altitude
+ !*******************************************************************
+
+ do kz=1,numzgrid
+ if (kz.eq.1) then
+ halfheight=outheight(1)/2.
+ else
+ halfheight=(outheight(kz)+outheight(kz-1))/2.
+ endif
+ do kzz=2,nz
+ if ((height(kzz-1).lt.halfheight).and. &
+ (height(kzz).gt.halfheight)) goto 46
+ end do
+46 kzz=max(min(kzz,nz),2)
+ dz1=halfheight-height(kzz-1)
+ dz2=height(kzz)-halfheight
+ dz=dz1+dz2
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ xl=outlon0+real(ix)*dxout
+ yl=outlat0+real(jy)*dyout
+ xl=(xl-xlon0)/dx
+ yl=(yl-ylat0)/dx
+ iix=max(min(nint(xl),nxmin1),0)
+ jjy=max(min(nint(yl),nymin1),0)
+ densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+ rho(iix,jjy,kzz-1,2)*dz2)/dz
+ end do
+ end do
+ end do
+
+ do i=1,numreceptor
+ xl=xreceptor(i)
+ yl=yreceptor(i)
+ iix=max(min(nint(xl),nxmin1),0)
+ jjy=max(min(nint(yl),nymin1),0)
+ densityoutrecept(i)=rho(iix,jjy,1,2)
+ end do
+
+
+ ! Output is different for forward and backward simulations
+ do kz=1,numzgrid
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (ldirect.eq.1) then
+ factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
+ else
+ factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
+ endif
+ end do
+ end do
+ end do
+
+ !*********************************************************************
+ ! Determine the standard deviation of the mean concentration or mixing
+ ! ratio (uncertainty of the output) and the dry and wet deposition
+ !*********************************************************************
+
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf 3 (sd)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+ gridtotal=0.
+ gridsigmatotal=0.
+ gridtotalunc=0.
+ wetgridtotal=0.
+ wetgridsigmatotal=0.
+ wetgridtotalunc=0.
+ drygridtotal=0.
+ drygridsigmatotal=0.
+ drygridtotalunc=0.
+
+ do ks=1,nspec
+
+ write(anspec,'(i3.3)') ks
+ if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+ if (ldirect.eq.1) then
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
+ atime//'_'//anspec,form='unformatted')
+ else
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
+ atime//'_'//anspec,form='unformatted')
+ endif
+ write(unitoutgrid) itime
+ endif
+
+ if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
+ open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
+ atime//'_'//anspec,form='unformatted')
+
+ write(unitoutgridppt) itime
+ endif
+
+ do kp=1,maxpointspec_act
+ do nage=1,nageclass
+
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+
+ ! WET DEPOSITION
+ if ((WETDEP).and.(ldirect.gt.0)) then
+ do l=1,nclassunc
+ auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,wetgrid(ix,jy), &
+ wetgridsigma(ix,jy),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ wetgrid(ix,jy)=wetgrid(ix,jy) &
+ *nclassunc
+ wetgridtotal=wetgridtotal+wetgrid(ix,jy)
+ ! Calculate standard deviation of the mean
+ wetgridsigma(ix,jy)= &
+ wetgridsigma(ix,jy)* &
+ sqrt(real(nclassunc))
+ wetgridsigmatotal=wetgridsigmatotal+ &
+ wetgridsigma(ix,jy)
+ endif
+
+ ! DRY DEPOSITION
+ if ((DRYDEP).and.(ldirect.gt.0)) then
+ do l=1,nclassunc
+ auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,drygrid(ix,jy), &
+ drygridsigma(ix,jy),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ drygrid(ix,jy)=drygrid(ix,jy)* &
+ nclassunc
+ drygridtotal=drygridtotal+drygrid(ix,jy)
+ ! Calculate standard deviation of the mean
+ drygridsigma(ix,jy)= &
+ drygridsigma(ix,jy)* &
+ sqrt(real(nclassunc))
+125 drygridsigmatotal=drygridsigmatotal+ &
+ drygridsigma(ix,jy)
+ endif
+
+ ! CONCENTRATION OR MIXING RATIO
+ do kz=1,numzgrid
+ do l=1,nclassunc
+ auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,grid(ix,jy,kz), &
+ gridsigma(ix,jy,kz),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ grid(ix,jy,kz)= &
+ grid(ix,jy,kz)*nclassunc
+ gridtotal=gridtotal+grid(ix,jy,kz)
+ ! Calculate standard deviation of the mean
+ gridsigma(ix,jy,kz)= &
+ gridsigma(ix,jy,kz)* &
+ sqrt(real(nclassunc))
+ gridsigmatotal=gridsigmatotal+ &
+ gridsigma(ix,jy,kz)
+ end do
+ end do
+ end do
+
+
+
+
+ !*******************************************************************
+ ! Generate output: may be in concentration (ng/m3) or in mixing
+ ! ratio (ppt) or both
+ ! Output the position and the values alternated multiplied by
+ ! 1 or -1, first line is number of values, number of positions
+ ! For backward simulations, the unit is seconds, stored in grid_time
+ !*******************************************************************
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf 4 (output)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+
+ ! Concentration output
+ !*********************
+ if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf (Wet deposition)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+
+ ! Wet deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(WETDEP)) then
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ !oncentraion greater zero
+ if (wetgrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)=ix+jy*numxgrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf (Dry deposition)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+ ! Dry deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(DRYDEP)) then
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (drygrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)=ix+jy*numxgrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*drygrid(ix,jy)/area(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*drygridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf (Concentrations)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+ ! Concentrations
+
+ ! if surf_only write only 1st layer
+
+ if(surf_only.eq.1) then
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,1
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid+kz*numxgrid*numygrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ grid(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
+ ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
+ sparse_dump_u(sp_count_r)= &
+ gridsigma(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+ else
+
+ ! write full vertical resolution
+ if (verbosity.eq.1) then
+ print*,'concoutput_surf (write full vertical resolution)'
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,numzgrid
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid+kz*numxgrid*numygrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ grid(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
+ ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
+ sparse_dump_u(sp_count_r)= &
+ gridsigma(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+ endif ! surf_only
+
+ endif ! concentration output
+
+ ! Mixing ratio output
+ !********************
+
+ if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
+
+ ! Wet deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(WETDEP)) then
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (wetgrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*wetgrid(ix,jy)/area(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+ ! Dry deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(DRYDEP)) then
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (drygrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*drygrid(ix,jy)/area(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*drygridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+ ! Mixing ratios
+
+ ! if surf_only write only 1st layer
+
+ if(surf_only.eq.1) then
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,1
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid+kz*numxgrid*numygrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*grid(ix,jy,kz) &
+ /volume(ix,jy,kz)/outnum* &
+ weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
+ outnum*weightair/weightmolar(ks)/ &
+ densityoutgrid(ix,jy,kz)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+ else
+
+ ! write full vertical resolution
+
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,numzgrid
+ do jy=0,numygrid-1
+ do ix=0,numxgrid-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgrid+kz*numxgrid*numygrid
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*grid(ix,jy,kz) &
+ /volume(ix,jy,kz)/outnum* &
+ weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
+ outnum*weightair/weightmolar(ks)/ &
+ densityoutgrid(ix,jy,kz)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+! write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+ endif ! surf_only
+
+ endif ! output for ppt
+
+ end do
+ end do
+
+ close(unitoutgridppt)
+ close(unitoutgrid)
+
+ end do
+
+ if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
+ if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
+ wetgridtotal
+ if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
+ drygridtotal
+
+ ! Dump of receptor concentrations
+
+ if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3) ) then
+ write(unitoutreceptppt) itime
+ do ks=1,nspec
+ write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
+ weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
+ end do
+ endif
+
+ ! Dump of receptor concentrations
+
+ if (numreceptor.gt.0) then
+ write(unitoutrecept) itime
+ do ks=1,nspec
+ write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
+ i=1,numreceptor)
+ end do
+ endif
+
+
+
+ ! 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,numygrid-1
+ do ix=0,numxgrid-1
+ do l=1,nclassunc
+ do nage=1,nageclass
+ do kz=1,numzgrid
+ gridunc(ix,jy,kz,ks,kp,l,nage)=0.
+ end do
+ end do
+ end do
+ end do
+ end do
+ end do
+ end do
+
+
+end subroutine concoutput_surf
Index: /trunk/src/concoutput_surf_nest.f90
===================================================================
--- /trunk/src/concoutput_surf_nest.f90 (revision 20)
+++ /trunk/src/concoutput_surf_nest.f90 (revision 20)
@@ -0,0 +1,652 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine concoutput_surf_nest(itime,outnum)
+ ! i i
+ !*****************************************************************************
+ ! *
+ ! Output of the concentration grid and the receptor concentrations. *
+ ! *
+ ! Author: A. Stohl *
+ ! *
+ ! 24 May 1995 *
+ ! *
+ ! 13 April 1999, Major update: if output size is smaller, dump output *
+ ! in sparse matrix format; additional output of *
+ ! uncertainty *
+ ! *
+ ! 05 April 2000, Major update: output of age classes; output for backward*
+ ! runs is time spent in grid cell times total mass of *
+ ! species. *
+ ! *
+ ! 17 February 2002, Appropriate dimensions for backward and forward runs *
+ ! are now specified in file par_mod *
+ ! *
+ ! June 2006, write grid in sparse matrix with a single write command *
+ ! in order to save disk space *
+ ! *
+ ! 2008 new sparse matrix format *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! outnum number of samples *
+ ! ncells number of cells with non-zero concentrations *
+ ! sparse .true. if in sparse matrix format, else .false. *
+ ! tot_mu 1 for forward, initial mass mixing ration for backw. runs *
+ ! *
+ !*****************************************************************************
+
+ use unc_mod
+ use point_mod
+ use outg_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ real(kind=dp) :: jul
+ integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
+ integer :: sp_count_i,sp_count_r
+ real :: sp_fact
+ real :: outnum,densityoutrecept(maxreceptor),xl,yl
+
+ !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
+ ! +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
+ ! + maxageclass)
+ !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
+ ! + maxageclass)
+ !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass)
+ !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
+ ! + maxpointspec_act,maxageclass),
+ ! + drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass),
+ ! + wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
+ ! + maxpointspec_act,maxageclass)
+ !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
+ !real sparse_dump_r(numxgrid*numygrid*numzgrid)
+ !integer sparse_dump_i(numxgrid*numygrid*numzgrid)
+
+ !real sparse_dump_u(numxgrid*numygrid*numzgrid)
+ real :: auxgrid(nclassunc)
+ real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
+ real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+ real,parameter :: weightair=28.97
+ logical :: sp_zer
+ character :: adate*8,atime*6
+ character(len=3) :: anspec
+
+
+ ! Determine current calendar date, needed for the file name
+ !**********************************************************
+
+ jul=bdate+real(itime,kind=dp)/86400._dp
+ call caldate(jul,jjjjmmdd,ihmmss)
+ write(adate,'(i8.8)') jjjjmmdd
+ write(atime,'(i6.6)') ihmmss
+
+
+ ! For forward simulations, output fields have dimension MAXSPEC,
+ ! for backward simulations, output fields have dimension MAXPOINT.
+ ! Thus, make loops either about nspec, or about numpoint
+ !*****************************************************************
+
+
+ if (ldirect.eq.1) then
+ do ks=1,nspec
+ do kp=1,maxpointspec_act
+ tot_mu(ks,kp)=1
+ end do
+ end do
+ else
+ do ks=1,nspec
+ do kp=1,maxpointspec_act
+ tot_mu(ks,kp)=xmass(kp,ks)
+ end do
+ end do
+ endif
+
+
+ !*******************************************************************
+ ! Compute air density: sufficiently accurate to take it
+ ! from coarse grid at some time
+ ! Determine center altitude of output layer, and interpolate density
+ ! data to that altitude
+ !*******************************************************************
+
+ do kz=1,numzgrid
+ if (kz.eq.1) then
+ halfheight=outheight(1)/2.
+ else
+ halfheight=(outheight(kz)+outheight(kz-1))/2.
+ endif
+ do kzz=2,nz
+ if ((height(kzz-1).lt.halfheight).and. &
+ (height(kzz).gt.halfheight)) goto 46
+ end do
+46 kzz=max(min(kzz,nz),2)
+ dz1=halfheight-height(kzz-1)
+ dz2=height(kzz)-halfheight
+ dz=dz1+dz2
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ xl=outlon0n+real(ix)*dxoutn
+ yl=outlat0n+real(jy)*dyoutn
+ xl=(xl-xlon0)/dx
+ yl=(yl-ylat0)/dy
+ iix=max(min(nint(xl),nxmin1),0)
+ jjy=max(min(nint(yl),nymin1),0)
+ densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+ rho(iix,jjy,kzz-1,2)*dz2)/dz
+ end do
+ end do
+ end do
+
+ do i=1,numreceptor
+ xl=xreceptor(i)
+ yl=yreceptor(i)
+ iix=max(min(nint(xl),nxmin1),0)
+ jjy=max(min(nint(yl),nymin1),0)
+ densityoutrecept(i)=rho(iix,jjy,1,2)
+ end do
+
+
+ ! Output is different for forward and backward simulations
+ do kz=1,numzgrid
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (ldirect.eq.1) then
+ factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
+ else
+ factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
+ endif
+ end do
+ end do
+ end do
+
+ !*********************************************************************
+ ! Determine the standard deviation of the mean concentration or mixing
+ ! ratio (uncertainty of the output) and the dry and wet deposition
+ !*********************************************************************
+
+ do ks=1,nspec
+
+ write(anspec,'(i3.3)') ks
+ if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+ if (ldirect.eq.1) then
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
+ //adate// &
+ atime//'_'//anspec,form='unformatted')
+ else
+ open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_' &
+ //adate// &
+ atime//'_'//anspec,form='unformatted')
+ endif
+ write(unitoutgrid) itime
+ endif
+
+ if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
+ open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' &
+ //adate// &
+ atime//'_'//anspec,form='unformatted')
+
+ write(unitoutgridppt) itime
+ endif
+
+ do kp=1,maxpointspec_act
+ do nage=1,nageclass
+
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+
+ ! WET DEPOSITION
+ if ((WETDEP).and.(ldirect.gt.0)) then
+ do l=1,nclassunc
+ auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,wetgrid(ix,jy), &
+ wetgridsigma(ix,jy),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ wetgrid(ix,jy)=wetgrid(ix,jy) &
+ *nclassunc
+ ! Calculate standard deviation of the mean
+ wetgridsigma(ix,jy)= &
+ wetgridsigma(ix,jy)* &
+ sqrt(real(nclassunc))
+ endif
+
+ ! DRY DEPOSITION
+ if ((DRYDEP).and.(ldirect.gt.0)) then
+ do l=1,nclassunc
+ auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,drygrid(ix,jy), &
+ drygridsigma(ix,jy),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ drygrid(ix,jy)=drygrid(ix,jy)* &
+ nclassunc
+ ! Calculate standard deviation of the mean
+ drygridsigma(ix,jy)= &
+ drygridsigma(ix,jy)* &
+ sqrt(real(nclassunc))
+ endif
+
+ ! CONCENTRATION OR MIXING RATIO
+ do kz=1,numzgrid
+ do l=1,nclassunc
+ auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
+ end do
+ call mean(auxgrid,grid(ix,jy,kz), &
+ gridsigma(ix,jy,kz),nclassunc)
+ ! Multiply by number of classes to get total concentration
+ grid(ix,jy,kz)= &
+ grid(ix,jy,kz)*nclassunc
+ ! Calculate standard deviation of the mean
+ gridsigma(ix,jy,kz)= &
+ gridsigma(ix,jy,kz)* &
+ sqrt(real(nclassunc))
+ end do
+ end do
+ end do
+
+
+ !*******************************************************************
+ ! Generate output: may be in concentration (ng/m3) or in mixing
+ ! ratio (ppt) or both
+ ! Output the position and the values alternated multiplied by
+ ! 1 or -1, first line is number of values, number of positions
+ ! For backward simulations, the unit is seconds, stored in grid_time
+ !*******************************************************************
+
+ ! Concentration output
+ !*********************
+ if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+
+ ! Wet deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(WETDEP)) then
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ !oncentraion greater zero
+ if (wetgrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)=ix+jy*numxgridn
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact*1.e12*wetgrid(ix,jy)/arean(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+
+ ! Dry deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(DRYDEP)) then
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (drygrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)=ix+jy*numxgridn
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*drygrid(ix,jy)/arean(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*drygridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+
+ ! Concentrations
+
+ ! if surf_only write only 1st layer
+
+ if(surf_only.eq.1) then
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,1
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ 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* &
+ grid(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
+ ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
+ sparse_dump_u(sp_count_r)= &
+ gridsigma(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+ else
+
+ ! write full vertical resolution
+
+ 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 (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ 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* &
+ grid(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ ! if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
+ ! + write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
+ sparse_dump_u(sp_count_r)= &
+ gridsigma(ix,jy,kz)* &
+ factor3d(ix,jy,kz)/tot_mu(ks,kp)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgrid) sp_count_i
+ write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgrid) sp_count_r
+ write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
+ endif ! surf_only
+
+
+ endif ! concentration output
+
+ ! Mixing ratio output
+ !********************
+
+ if ((iout.eq.2).or.(iout.eq.3)) then ! mixing ratio
+
+ ! Wet deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(WETDEP)) then
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (wetgrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgridn
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1.)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*wetgrid(ix,jy)/arean(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*wetgridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+ ! Dry deposition
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ if ((ldirect.eq.1).and.(DRYDEP)) then
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (drygrid(ix,jy).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ sp_count_i=sp_count_i+1
+ sparse_dump_i(sp_count_i)= &
+ ix+jy*numxgridn
+ sp_zer=.false.
+ sp_fact=sp_fact*(-1)
+ endif
+ sp_count_r=sp_count_r+1
+ sparse_dump_r(sp_count_r)= &
+ sp_fact* &
+ 1.e12*drygrid(ix,jy)/arean(ix,jy)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*drygridsigma(ix,jy)/area(ix,jy)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ else
+ sp_count_i=0
+ sp_count_r=0
+ endif
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+
+
+ ! Mixing ratios
+
+ ! if surf_only write only 1st layer
+
+ if(surf_only.eq.1) then
+ sp_count_i=0
+ sp_count_r=0
+ sp_fact=-1.
+ sp_zer=.true.
+ do kz=1,1
+ do jy=0,numygridn-1
+ do ix=0,numxgridn-1
+ if (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ 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* &
+ 1.e12*grid(ix,jy,kz) &
+ /volumen(ix,jy,kz)/outnum* &
+ weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*gridsigma(ix,jy,kz)/volumen(ix,jy,kz)/ &
+ outnum*weightair/weightmolar(ks)/ &
+ densityoutgrid(ix,jy,kz)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+ else
+
+ ! write full vertical resolution
+
+ 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 (grid(ix,jy,kz).gt.smallnum) then
+ if (sp_zer.eqv..true.) then ! first non zero value
+ 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* &
+ 1.e12*grid(ix,jy,kz) &
+ /volumen(ix,jy,kz)/outnum* &
+ weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
+ sparse_dump_u(sp_count_r)= &
+ 1.e12*gridsigma(ix,jy,kz)/volumen(ix,jy,kz)/ &
+ outnum*weightair/weightmolar(ks)/ &
+ densityoutgrid(ix,jy,kz)
+ else ! concentration is zero
+ sp_zer=.true.
+ endif
+ end do
+ end do
+ end do
+ write(unitoutgridppt) sp_count_i
+ write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+ write(unitoutgridppt) sp_count_r
+ write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
+ endif ! surf_only
+
+ endif ! output for ppt
+
+ end do
+ end do
+
+ close(unitoutgridppt)
+ close(unitoutgrid)
+
+ 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
+ end do
+ end do
+ end do
+
+
+end subroutine concoutput_surf_nest
+
Index: /trunk/src/gridcheck.f90
===================================================================
--- /trunk/src/gridcheck.f90 (revision 19)
+++ /trunk/src/gridcheck.f90 (revision 20)
@@ -98,5 +98,6 @@
integer :: isec1(56),isec2(22+nxmax+nymax)
- real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+ !real(kind=4) :: zsec2(60+2*nuvzmax),zsec4(jpunp)
+ real(kind=4) :: zsec2(184),zsec4(jpunp)
character(len=1) :: opt
@@ -466,5 +467,4 @@
k=nlev_ec+1+numskip+i
akm(nwz-i+1)=zsec2(j)
- ! write (*,*) 'ifield:',ifield,k,j,zsec2(10+j)
bkm(nwz-i+1)=zsec2(k)
end do
Index: /trunk/src/gridcheck_fnl.f90
===================================================================
--- /trunk/src/gridcheck_fnl.f90 (revision 20)
+++ /trunk/src/gridcheck_fnl.f90 (revision 20)
@@ -0,0 +1,545 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine gridcheck
+
+ !**********************************************************************
+ ! *
+ ! FLEXPART MODEL SUBROUTINE GRIDCHECK *
+ ! *
+ !**********************************************************************
+ ! *
+ ! AUTHOR: G. WOTAWA *
+ ! DATE: 1997-08-06 *
+ ! LAST UPDATE: 1997-10-10 *
+ ! *
+ ! Update: 1999-02-08, global fields allowed, A. Stohl*
+ ! CHANGE: 17/11/2005, Caroline Forster, GFS data *
+ ! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
+ ! ECMWF grib_api *
+ ! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
+ ! ECMWF grib_api *
+ ! *
+ !**********************************************************************
+ ! *
+ ! DESCRIPTION: *
+ ! *
+ ! THIS SUBROUTINE DETERMINES THE GRID SPECIFICATIONS (LOWER LEFT *
+ ! LONGITUDE, LOWER LEFT LATITUDE, NUMBER OF GRID POINTS, GRID DIST- *
+ ! ANCE AND VERTICAL DISCRETIZATION OF THE ECMWF MODEL) FROM THE *
+ ! GRIB HEADER OF THE FIRST INPUT FILE. THE CONSISTANCY (NO CHANGES *
+ ! WITHIN ONE FLEXPART RUN) IS CHECKED IN THE ROUTINE "READWIND" AT *
+ ! ANY CALL. *
+ ! *
+ ! XLON0 geographical longitude of lower left gridpoint *
+ ! YLAT0 geographical latitude of lower left gridpoint *
+ ! NX number of grid points x-direction *
+ ! NY number of grid points y-direction *
+ ! DX grid distance x-direction *
+ ! DY grid distance y-direction *
+ ! NUVZ number of grid points for horizontal wind *
+ ! components in z direction *
+ ! NWZ number of grid points for vertical wind *
+ ! component in z direction *
+ ! sizesouth, sizenorth give the map scale (i.e. number of virtual grid*
+ ! points of the polar stereographic grid): *
+ ! used to check the CFL criterion *
+ ! VHEIGHT(1)- heights of gridpoints where u and v are *
+ ! UVHEIGHT(NUVZ) given *
+ ! WHEIGHT(1)- heights of gridpoints where w is given *
+ ! WHEIGHT(NWZ) *
+ ! *
+ !**********************************************************************
+
+ use grib_api
+ use par_mod
+ use com_mod
+ use conv_mod
+ use cmapf_mod, only: stlmbr,stcm2p
+
+ implicit none
+
+ !HSO parameters for grib_api
+ integer :: ifile
+ integer :: iret
+ integer :: igrib
+ real(kind=4) :: xaux1,xaux2,yaux1,yaux2
+ real(kind=8) :: xaux1in,xaux2in,yaux1in,yaux2in
+ integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
+ !HSO end
+ integer :: ix,jy,i,ifn,ifield,j,k,iumax,iwmax,numskip
+ real :: sizesouth,sizenorth,xauxa,pint
+ real :: akm_usort(nwzmax)
+ real,parameter :: eps=0.0001
+
+ ! NCEP GFS
+ real :: pres(nwzmax), help
+
+ integer :: i179,i180,i181
+
+ ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
+
+ integer :: isec1(8),isec2(3)
+ real(kind=4) :: zsec4(jpunp)
+ character(len=1) :: opt
+
+ !HSO grib api error messages
+ character(len=24) :: gribErrorMsg = 'Error reading grib file'
+ character(len=20) :: gribFunction = 'gridcheckwind_gfs'
+ !
+ if (numbnests.ge.1) then
+ write(*,*) ' ###########################################'
+ write(*,*) ' FLEXPART ERROR SUBROUTINE GRIDCHECK:'
+ write(*,*) ' NO NESTED WINDFIELDAS ALLOWED FOR GFS! '
+ write(*,*) ' ###########################################'
+ stop
+ endif
+
+ iumax=0
+ iwmax=0
+
+ if(ideltas.gt.0) then
+ ifn=1
+ else
+ ifn=numbwf
+ endif
+ !
+ ! OPENING OF DATA FILE (GRIB CODE)
+ !
+5 call grib_open_file(ifile,path(3)(1:length(3)) &
+ //trim(wfname(ifn)),'r',iret)
+ if (iret.ne.GRIB_SUCCESS) then
+ goto 999 ! ERROR DETECTED
+ endif
+ !turn on support for multi fields messages
+ call grib_multi_support_on
+
+ ifield=0
+10 ifield=ifield+1
+ !
+ ! GET NEXT FIELDS
+ !
+ call grib_new_from_file(ifile,igrib,iret)
+ if (iret.eq.GRIB_END_OF_FILE ) then
+ goto 30 ! EOF DETECTED
+ elseif (iret.ne.GRIB_SUCCESS) then
+ goto 999 ! ERROR DETECTED
+ endif
+
+ !first see if we read GRIB1 or GRIB2
+ call grib_get_int(igrib,'editionNumber',gribVer,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ if (gribVer.eq.1) then ! GRIB Edition 1
+
+ !read the grib1 identifiers
+ call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'indicatorOfTypeOfLevel',isec1(7),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'level',isec1(8),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ !get the size and data of the values array
+ call grib_get_real4_array(igrib,'values',zsec4,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ else ! GRIB Edition 2
+
+ !read the grib2 identifiers
+ call grib_get_int(igrib,'discipline',discipl,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'parameterCategory',parCat,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'parameterNumber',parNum,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'scaledValueOfFirstFixedSurface', &
+ valSurf,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ !convert to grib1 identifiers
+ isec1(6)=-1
+ isec1(7)=-1
+ isec1(8)=-1
+ if ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.100)) then ! U
+ isec1(6)=33 ! indicatorOfParameter
+ isec1(7)=100 ! indicatorOfTypeOfLevel
+ isec1(8)=valSurf/100 ! level, convert to hPa
+ elseif ((parCat.eq.3).and.(parNum.eq.5).and.(typSurf.eq.1)) then ! TOPO
+ isec1(6)=7 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.1) &
+ .and.(discipl.eq.2)) then ! LSM
+ isec1(6)=81 ! indicatorOfParameter
+ isec1(7)=1 ! indicatorOfTypeOfLevel
+ isec1(8)=0
+ endif
+
+ if (isec1(6).ne.-1) then
+ ! get the size and data of the values array
+ call grib_get_real4_array(igrib,'values',zsec4,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ endif
+
+ endif ! gribVer
+
+ if(ifield.eq.1) then
+
+ !get the required fields from section 2
+ !store compatible to gribex input
+ call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
+ isec2(2),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
+ isec2(3),iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
+ xaux1in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'longitudeOfLastGridPointInDegrees', &
+ xaux2in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
+ yaux1in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+ call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
+ yaux2in,iret)
+ call grib_check(iret,gribFunction,gribErrorMsg)
+
+ xaux1=xaux1in
+! xaux2=xaux2in
+ xaux2=xaux2in+360. !!! Transform FNL xauu2in from -1 to 359 (Xuekue Fang, 31 Jan 2013)
+ yaux1=yaux1in
+ yaux2=yaux2in
+ write(*,*) 'xaux1,xaux2,yaux1,yaux2',xaux1,xaux2,yaux1,yaux2
+ nxfield=isec2(2)
+ ny=isec2(3)
+
+
+
+ if((abs(xaux1).lt.eps).and.(xaux2.ge.359)) then ! NCEP DATA FROM 0 TO
+ xaux1=-179.0 ! 359 DEG EAST ->
+ xaux2=-179.0+360.-360./real(nxfield) ! TRANSFORMED TO -179
+ endif ! TO 180 DEG EAST
+ if (xaux1.gt.180) xaux1=xaux1-360.0
+ if (xaux2.gt.180) xaux2=xaux2-360.0
+ if (xaux1.lt.-180) xaux1=xaux1+360.0
+ if (xaux2.lt.-180) xaux2=xaux2+360.0
+ if (xaux2.lt.xaux1) xaux2=xaux2+360.
+ xlon0=xaux1
+ ylat0=yaux1
+ write(*,*) 'xlon0,ylat0',xlon0,ylat0
+
+ dx=(xaux2-xaux1)/real(nxfield-1)
+ dy=(yaux2-yaux1)/real(ny-1)
+ dxconst=180./(dx*r_earth*pi)
+ dyconst=180./(dy*r_earth*pi)
+ write(*,*) 'dx,dy,dxconst,dyconst',dx,dy,dxconst,dyconst
+ !HSO end edits
+
+
+ ! Check whether fields are global
+ ! If they contain the poles, specify polar stereographic map
+ ! projections using the stlmbr- and stcm2p-calls
+ !***********************************************************
+
+ xauxa=abs(xaux2+dx-360.-xaux1)
+ if (xauxa.lt.0.001) then
+ nx=nxfield+1 ! field is cyclic
+ xglobal=.true.
+ if (abs(nxshift).ge.nx) &
+ stop 'nxshift in file par_mod is too large'
+ xlon0=xlon0+real(nxshift)*dx
+ else
+ nx=nxfield
+ xglobal=.false.
+ if (nxshift.ne.0) &
+ stop 'nxshift (par_mod) must be zero for non-global domain'
+ endif
+ nxmin1=nx-1
+ nymin1=ny-1
+ if (xlon0.gt.180.) xlon0=xlon0-360.
+ xauxa=abs(yaux1+90.)
+ if (xglobal.and.xauxa.lt.0.001) then
+ sglobal=.true. ! field contains south pole
+ ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+ ! map scale
+ sizesouth=6.*(switchsouth+90.)/dy
+ call stlmbr(southpolemap,-90.,0.)
+ call stcm2p(southpolemap,0.,0.,switchsouth,0.,sizesouth, &
+ sizesouth,switchsouth,180.)
+ switchsouthg=(switchsouth-ylat0)/dy
+ else
+ sglobal=.false.
+ switchsouthg=999999.
+ endif
+ xauxa=abs(yaux2-90.)
+ if (xglobal.and.xauxa.lt.0.001) then
+ nglobal=.true. ! field contains north pole
+ ! Enhance the map scale by factor 3 (*2=6) compared to north-south
+ ! map scale
+ sizenorth=6.*(90.-switchnorth)/dy
+ call stlmbr(northpolemap,90.,0.)
+ call stcm2p(northpolemap,0.,0.,switchnorth,0.,sizenorth, &
+ sizenorth,switchnorth,180.)
+ switchnorthg=(switchnorth-ylat0)/dy
+ else
+ nglobal=.false.
+ switchnorthg=999999.
+ endif
+ endif ! ifield.eq.1
+
+ if (nxshift.lt.0) stop 'nxshift (par_mod) must not be negative'
+ if (nxshift.ge.nxfield) stop 'nxshift (par_mod) too large'
+
+ ! NCEP ISOBARIC LEVELS
+ !*********************
+
+ if((isec1(6).eq.33).and.(isec1(7).eq.100)) then ! check for U wind
+ iumax=iumax+1
+ pres(iumax)=real(isec1(8))*100.0
+ endif
+
+
+ i179=nint(179./dx)
+ if (dx.lt.0.7) then
+ i180=nint(180./dx)+1 ! 0.5 deg data
+ else
+ i180=nint(179./dx)+1 ! 1 deg data
+ endif
+ i181=i180+1
+
+
+ ! NCEP TERRAIN
+ !*************
+
+ if((isec1(6).eq.007).and.(isec1(7).eq.001)) then
+ do jy=0,ny-1
+ do ix=0,nxfield-1
+ help=zsec4(nxfield*(ny-jy-1)+ix+1)
+ if(ix.le.i180) then
+ oro(i179+ix,jy)=help
+ excessoro(i179+ix,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ else
+ oro(ix-i181,jy)=help
+ excessoro(ix-i181,jy)=0.0 ! ISOBARIC SURFACES: SUBGRID TERRAIN DISREGARDED
+ endif
+ end do
+ end do
+ endif
+
+ ! NCEP LAND SEA MASK
+ !*******************
+
+ if((isec1(6).eq.081).and.(isec1(7).eq.001)) then
+ do jy=0,ny-1
+ do ix=0,nxfield-1
+ help=zsec4(nxfield*(ny-jy-1)+ix+1)
+ if(ix.le.i180) then
+ lsm(i179+ix,jy)=help
+ else
+ lsm(ix-i181,jy)=help
+ endif
+ end do
+ end do
+ endif
+
+ call grib_release(igrib)
+
+ goto 10 !! READ NEXT LEVEL OR PARAMETER
+ !
+ ! CLOSING OF INPUT DATA FILE
+ !
+
+ ! HSO
+30 continue
+ call grib_close_file(ifile)
+ ! HSO end edits
+
+ nuvz=iumax
+ nwz =iumax
+ nlev_ec=iumax
+
+ if (nx.gt.nxmax) then
+ write(*,*) 'FLEXPART error: Too many grid points in x direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) nx,nxmax
+ stop
+ endif
+
+ if (ny.gt.nymax) then
+ write(*,*) 'FLEXPART error: Too many grid points in y direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) ny,nymax
+ stop
+ endif
+
+ if (nuvz.gt.nuvzmax) then
+ write(*,*) 'FLEXPART error: Too many u,v grid points in z '// &
+ 'direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) nuvz,nuvzmax
+ stop
+ endif
+
+ if (nwz.gt.nwzmax) then
+ write(*,*) 'FLEXPART error: Too many w grid points in z '// &
+ 'direction.'
+ write(*,*) 'Reduce resolution of wind fields.'
+ write(*,*) 'Or change parameter settings in file par_mod.'
+ write(*,*) nwz,nwzmax
+ stop
+ endif
+
+ ! If desired, shift all grids by nxshift grid cells
+ !**************************************************
+
+ if (xglobal) then
+ call shift_field_0(oro,nxfield,ny)
+ call shift_field_0(lsm,nxfield,ny)
+ call shift_field_0(excessoro,nxfield,ny)
+ endif
+
+ ! Output of grid info
+ !********************
+
+ write(*,*)
+ write(*,*)
+ write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
+ nuvz,nwz
+ write(*,*)
+ write(*,'(a)') 'other domain:'
+ write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Longitude range: ', &
+ xlon0,' to ',xlon0+(nx-1)*dx,' Grid distance: ',dx
+ write(*,'(a,f10.2,a1,f10.2,a,f10.2)') ' Latitude range: ', &
+ ylat0,' to ',ylat0+(ny-1)*dy,' Grid distance: ',dy
+ write(*,*)
+
+
+ ! CALCULATE VERTICAL DISCRETIZATION OF ECMWF MODEL
+ ! PARAMETER akm,bkm DESCRIBE THE HYBRID "ETA" COORDINATE SYSTEM
+
+ numskip=nlev_ec-nuvz ! number of ecmwf model layers not used
+ ! by trajectory model
+ do i=1,nwz
+ j=numskip+i
+ k=nlev_ec+1+numskip+i
+ akm_usort(nwz-i+1)=pres(nwz-i+1)
+ bkm(nwz-i+1)=0.0
+ end do
+
+ !******************************
+ ! change Sabine Eckhardt: akm should always be in descending order ... readwind adapted!
+ !******************************
+ do i=1,nwz
+ if (akm_usort(1).gt.akm_usort(2)) then
+ akm(i)=akm_usort(i)
+ else
+ akm(i)=akm_usort(nwz-i+1)
+ endif
+ end do
+
+ !
+ ! CALCULATION OF AKZ, BKZ
+ ! AKZ,BKZ: model discretization parameters at the center of each model
+ ! layer
+ !
+ ! Assign the 10 m winds to an artificial model level with akz=0 and bkz=1.0,
+ ! i.e. ground level
+ !*****************************************************************************
+
+ do i=1,nuvz
+ akz(i)=akm(i)
+ bkz(i)=bkm(i)
+ end do
+
+ ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+ ! upon the transformation to z levels. In order to save computer memory, this is
+ ! not done anymore in the standard version. However, this option can still be
+ ! switched on by replacing the following lines with those below, that are
+ ! currently commented out. For this, similar changes are necessary in
+ ! verttransform.f and verttranform_nests.f
+ !*****************************************************************************
+
+ nz=nuvz
+ if (nz.gt.nzmax) stop 'nzmax too small'
+ do i=1,nuvz
+ aknew(i)=akz(i)
+ bknew(i)=bkz(i)
+ end do
+
+ ! Switch on following lines to use doubled vertical resolution
+ !*************************************************************
+ !nz=nuvz+nwz-1
+ !if (nz.gt.nzmax) stop 'nzmax too small'
+ !do 100 i=1,nwz
+ ! aknew(2*(i-1)+1)=akm(i)
+ !00 bknew(2*(i-1)+1)=bkm(i)
+ !do 110 i=2,nuvz
+ ! aknew(2*(i-1))=akz(i)
+ !10 bknew(2*(i-1))=bkz(i)
+ ! End doubled vertical resolution
+
+
+ ! Determine the uppermost level for which the convection scheme shall be applied
+ ! by assuming that there is no convection above 50 hPa (for standard SLP)
+ !*****************************************************************************
+
+ do i=1,nuvz-2
+ pint=akz(i)+bkz(i)*101325.
+ if (pint.lt.5000.) goto 96
+ end do
+96 nconvlev=i
+ if (nconvlev.gt.nconvlevmax-1) then
+ nconvlev=nconvlevmax-1
+ write(*,*) 'Attention, convection only calculated up to ', &
+ akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
+ endif
+
+ return
+
+999 write(*,*)
+ write(*,*) ' ###########################################'// &
+ '###### '
+ write(*,*) ' TRAJECTORY MODEL SUBROUTINE GRIDCHECK:'
+ write(*,*) ' CAN NOT OPEN INPUT DATA FILE '//wfname(ifn)
+ write(*,*) ' ###########################################'// &
+ '###### '
+ write(*,*)
+ write(*,'(a)') '!!! PLEASE INSERT A NEW CD-ROM AND !!!'
+ write(*,'(a)') '!!! PRESS ANY KEY TO CONTINUE... !!!'
+ write(*,'(a)') '!!! ...OR TERMINATE FLEXPART PRESSING!!!'
+ write(*,'(a)') '!!! THE "X" KEY... !!!'
+ write(*,*)
+ read(*,'(a)') opt
+ if(opt.eq.'X') then
+ stop
+ else
+ goto 5
+ endif
+
+end subroutine gridcheck
Index: /trunk/src/interpol_rain.f90
===================================================================
--- /trunk/src/interpol_rain.f90 (revision 19)
+++ /trunk/src/interpol_rain.f90 (revision 20)
@@ -20,8 +20,15 @@
!**********************************************************************
-subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
- ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
- ! i i i i i i i
- !i i i i i i i i o o o
+!subroutine interpol_rain(yy1,yy2,yy3,nxmax,nymax,nzmax,nx, &
+! ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3)
+! ! i i i i i i i
+! !i i i i i i i i o o o
+
+ subroutine interpol_rain(yy1,yy2,yy3,iy1,iy2,nxmax,nymax,nzmax,nx, &
+ ny,memind,xt,yt,level,itime1,itime2,itime,yint1,yint2,yint3, &
+ intiy1,intiy2,icmv)
+! i i i i i i i
+! i i i i i i i i o o o
+
!****************************************************************************
! *
@@ -40,5 +47,8 @@
! *
! 30 August 1996 *
- ! *
+ !
+ !* Petra Seibert, 2011/2012:
+ !* Add interpolation of cloud bottom and cloud thickness
+ !* for fix to SE's new wet scavenging scheme *
!****************************************************************************
! *
@@ -70,10 +80,15 @@
integer :: nx,ny,nxmax,nymax,nzmax,memind(2),m,ix,jy,ixp,jyp
- integer :: itime,itime1,itime2,level,indexh
+ !integer :: itime,itime1,itime2,level,indexh
+ integer :: itime,itime1,itime2,level,indexh,ip1,ip2,ip3,ip4
+ integer :: intiy1,intiy2,ipsum,icmv
real :: yy1(0:nxmax-1,0:nymax-1,nzmax,2)
real :: yy2(0:nxmax-1,0:nymax-1,nzmax,2)
real :: yy3(0:nxmax-1,0:nymax-1,nzmax,2)
- real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
- real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
+ integer iy1(0:nxmax-1,0:nymax-1,2),iy2(0:nxmax-1,0:nymax-1,2)
+ real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2),yi1(2),yi2(2)
+ !real :: ddx,ddy,rddx,rddy,dt1,dt2,dt,y1(2),y2(2),y3(2)
+ real :: xt,yt,yint1,yint2,yint3,yint4,p1,p2,p3,p4
+ !real :: xt,yt,yint1,yint2,yint3,p1,p2,p3,p4
@@ -127,4 +142,43 @@
+ p3*yy3(ix ,jyp,level,indexh) &
+ p4*yy3(ixp,jyp,level,indexh)
+
+!CPS clouds:
+ ip1=1
+ ip2=1
+ ip3=1
+ ip4=1
+ if (iy1(ix ,jy ,indexh) .eq. icmv) ip1=0
+ if (iy1(ixp,jy ,indexh) .eq. icmv) ip2=0
+ if (iy1(ix ,jyp,indexh) .eq. icmv) ip3=0
+ if (iy1(ixp,jyp,indexh) .eq. icmv) ip4=0
+ ipsum= ip1+ip2+ip3+ip4
+ if (ipsum .eq. 0) then
+ yi1(m)=icmv
+ else
+ yi1(m)=(ip1*p1*iy1(ix ,jy ,indexh) &
+ + ip2*p2*iy1(ixp,jy ,indexh) &
+ + ip3*p3*iy1(ix ,jyp,indexh) &
+ + ip4*p4*iy1(ixp,jyp,indexh))/ipsum
+ endif
+
+ ip1=1
+ ip2=1
+ ip3=1
+ ip4=1
+ if (iy2(ix ,jy ,indexh) .eq. icmv) ip1=0
+ if (iy2(ixp,jy ,indexh) .eq. icmv) ip2=0
+ if (iy2(ix ,jyp,indexh) .eq. icmv) ip3=0
+ if (iy2(ixp,jyp,indexh) .eq. icmv) ip4=0
+ ipsum= ip1+ip2+ip3+ip4
+ if (ipsum .eq. 0) then
+ yi2(m)=icmv
+ else
+ yi2(m)=(ip1*p1*iy2(ix ,jy ,indexh) &
+ + ip2*p2*iy2(ixp,jy ,indexh) &
+ + ip3*p3*iy2(ix ,jyp,indexh) &
+ + ip4*p4*iy2(ixp,jyp,indexh))/ipsum
+ endif
+!CPS end clouds
+
end do
@@ -134,4 +188,11 @@
!************************************
+ if (abs(itime) .lt. abs(itime1)) then
+ print*,'interpol_rain.f90'
+ print*,itime,itime1,itime2
+ stop 'ITIME PROBLEM'
+ endif
+
+
dt1=real(itime-itime1)
dt2=real(itime2-itime)
@@ -143,3 +204,20 @@
+!PS clouds:
+ intiy1=(yi1(1)*dt2 + yi1(2)*dt1)/dt
+ if (yi1(1) .eq. float(icmv)) intiy1=yi1(2)
+ if (yi1(2) .eq. float(icmv)) intiy1=yi1(1)
+
+ intiy2=(yi2(1)*dt2 + yi2(2)*dt1)/dt
+ if (yi2(1) .eq. float(icmv)) intiy2=yi2(2)
+ if (yi2(2) .eq. float(icmv)) intiy2=yi2(1)
+
+ if (intiy1 .ne. icmv .and. intiy2 .ne. icmv) then
+ intiy2 = intiy1 + intiy2 ! convert cloud thickness to cloud top
+ else
+ intiy1=icmv
+ intiy2=icmv
+ endif
+!PS end clouds
+
end subroutine interpol_rain
Index: /trunk/src/makefile
===================================================================
--- /trunk/src/makefile (revision 20)
+++ /trunk/src/makefile (revision 20)
@@ -0,0 +1,119 @@
+SHELL = /bin/bash
+TARGET = local
+WINDS=ecmwf
+#WINDS=gfs
+#WINDS=fnl
+
+FC = gfortran
+FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+
+ifeq ($(TARGET),dmz)
+ # options for ganglia
+ INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+ LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+ LIBPATH2 = /usr/lib/x86_64-linux-gnu/
+ MAIN = FLEXPART_dmz_
+endif
+ifeq ($(TARGET),local)
+ # local options
+ #libs_dir=/.../flexpart/libs/
+ libs_dir=/Users/ignacio/flexpart/libs/
+ INCPATH = $(libs_dir)/grib_api-1.9.9_dir/include
+ LIBPATH1 = $(libs_dir)/grib_api-1.9.9_dir/lib
+ LIBPATH2 = $(libs_dir)/jasper_dir/lib
+ MAIN = FLEXPART_local_
+endif
+
+LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
+
+MODOBJS = \
+par_mod.o com_mod.o \
+conv_mod.o hanna_mod.o \
+interpol_mod.o cmapf_mod.o \
+unc_mod.o oh_mod.o \
+xmass_mod.o flux_mod.o \
+point_mod.o outg_mod.o
+
+OBJECTS = \
+writeheader.o writeheader_txt.o writeheader_surf.o assignland.o\
+ part0.o \
+caldate.o partdep.o \
+coordtrafo.o psih.o \
+raerod.o \
+drydepokernel.o random.o \
+erf.o readavailable.o \
+ew.o readcommand.o \
+advance.o readdepo.o \
+releaseparticles.o psim.o \
+FLEXPART.o readlanduse.o \
+getfields.o init_domainfill.o\
+interpol_wind.o readoutgrid.o \
+interpol_all.o readpaths.o \
+getrb.o readreceptors.o \
+getrc.o readreleases.o \
+getvdep.o readspecies.o \
+interpol_misslev.o \
+conccalc.o \
+concoutput.o concoutput_surf.o scalev.o \
+pbl_profile.o readOHfield.o\
+juldate.o timemanager.o \
+interpol_vdep.o interpol_rain.o \
+partoutput.o \
+hanna.o wetdepokernel.o \
+mean.o wetdepo.o \
+hanna_short.o windalign.o \
+hanna1.o initialize.o \
+ calcpar_nests.o \
+verttransform_nests.o interpol_all_nests.o \
+interpol_wind_nests.o interpol_misslev_nests.o \
+interpol_vdep_nests.o interpol_rain_nests.o \
+getvdep_nests.o gridcheck_nests.o \
+readwind_nests.o \
+readageclasses.o readpartpositions.o \
+calcfluxes.o fluxoutput.o \
+qvsat.o skplin.o \
+convect43c.o redist.o \
+sort2.o distance.o \
+centerofmass.o plumetraj.o \
+openouttraj.o calcpv.o \
+calcpv_nests.o distance2.o \
+clustering.o interpol_wind_short.o \
+interpol_wind_short_nests.o shift_field_0.o \
+shift_field.o outgrid_init.o \
+openreceptors.o boundcond_domainfill.o\
+partoutput_short.o readoutgrid_nest.o \
+outgrid_init_nest.o writeheader_nest.o writeheader_nest_surf.o \
+concoutput_nest.o concoutput_surf_nest.o wetdepokernel_nest.o \
+drydepokernel_nest.o zenithangle.o \
+ohreaction.o getvdep_nests.o \
+initial_cond_calc.o initial_cond_output.o \
+dynamic_viscosity.o get_settling.o
+
+
+ifeq ($(WINDS),ecmwf)
+ OBJECTS_WINDS = \
+ calcpar.o readwind.o \
+ richardson.o verttransform.o \
+ obukhov.o gridcheck.o \
+ convmix.o calcmatrix.o
+endif
+
+ifeq ($(WINDS),gfs)
+ OBJECTS_WINDS = \
+ calcpar_gfs.o readwind_gfs.o \
+ richardson_gfs.o verttransform_gfs.o \
+ obukhov_gfs.o gridcheck_gfs.o \
+ convmix_gfs.o calcmatrix_gfs.o
+endif
+
+$(MAIN): $(MODOBJS) $(OBJECTS) $(OBJECTS_WINDS)
+ $(FC) *.o -o $(MAIN)_$(WINDS) $(LDFLAGS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+ $(FC) -c $(FFLAGS) $<
+
+clean:
+ rm *.o *.mod
+
Index: unk/src/makefile.ecmwf_absoft
===================================================================
--- /trunk/src/makefile.ecmwf_absoft (revision 19)
+++ (revision )
@@ -1,92 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART
-INCF = incl*
-#
-
-FC = /opt/absoft10.1/bin/f90
-#FC = /opt/absoft11.5/bin/f95
-INCPATH = /nilu2/home/flexpart/lib64/absoft/include
-LIBPATH1 = /nilu2/home/flexpart/lib64/absoft/lib
-LIBPATH2 = /nilu2/home/flexpart/lib64/absoft/lib
-FFLAGS = -O2 -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-#FFLAGS = -g -Rb -Rc -Rs -s -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
- raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
-
Index: unk/src/makefile.ecmwf_gfortran
===================================================================
--- /trunk/src/makefile.ecmwf_gfortran (revision 19)
+++ (revision )
@@ -1,94 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC = gfortran
-INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-LIBPATH2 = /usr/lib/x86_64-linux-gnu/
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-#FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: unk/src/makefile.ecmwf_gfortran_jfb
===================================================================
--- /trunk/src/makefile.ecmwf_gfortran_jfb (revision 19)
+++ (revision )
@@ -1,91 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC = gfortran
-INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /xnilu_wrk/flex_wrk/lib64/jasper/lib
-#FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: unk/src/makefile.gfs_absoft
===================================================================
--- /trunk/src/makefile.gfs_absoft (revision 19)
+++ (revision )
@@ -1,93 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_GFS
-INCF = incl*
-#
-
-
-FC = /opt/absoft10.1/bin/f95
-INCPATH = /nilu2/home/flexpart/lib64/absoft/include
-LIBPATH1 = /nilu2/home/flexpart/lib64/absoft/lib
-LIBPATH2 = /nilu2/home/flexpart/lib64/absoft/lib
-FFLAGS = -s -O2 -I$(INCPATH) -p$(INCPATH) -m64 -mcmodel=medium
-#FFLAGS = -g -Rb -Rc -Rs -s -B108 -YEXT_NAMES=LCS -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: unk/src/makefile.gfs_gfortran
===================================================================
--- /trunk/src/makefile.gfs_gfortran (revision 19)
+++ (revision )
@@ -1,90 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_GFS_GFORTRAN
-INCF = incl*
-#
-
-FC = gfortran
-INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar_gfs.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind_gfs.o \
-conccalc.o richardson_gfs.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform_gfs.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov_gfs.o gridcheck_gfs.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix_gfs.o calcmatrix_gfs.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: unk/src/makefile.henrik
===================================================================
--- /trunk/src/makefile.henrik (revision 19)
+++ (revision )
@@ -1,100 +1,0 @@
-SHELL = /bin/bash
-MAIN = FLEXPART_ICEALOT
-#
-
-FC = gfortran
-## OUTSIDE TEST PATH
-#INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api_test/include
-#LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api_test/lib
-#LIBPATH2 = /xnilu_wrk/flex_wrk/bin64/jasper/lib
-# OUTSIDE PATH
-INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-LIBPATH2 = /xnilu_wrk/flex_wrk/bin64/jasper/lib
-# INSIDE PATH
-#INCPATH = /nilu2/home/flexpart/lib64/gfortran/include
-#LIBPATH1 = /nilu2/home/flexpart/lib64/gfortran/lib
-#LIBPATH2 = /nilu2/home/flexpart/lib64/gfortran/lib
-#FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: unk/src/makefile.test
===================================================================
--- /trunk/src/makefile.test (revision 19)
+++ (revision )
@@ -1,96 +1,0 @@
-SHELL = /bin/bash
-MAIN = FP_ecmwf_gfortran
-#
-
-FC = gfortran
-#INCPATH = /xnilu_wrk/flex_wrk/bin64/grib_api/include
-#LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-INCPATH = /usr/lib/
-LIBPATH1 = /usr/lib/
-LIBPATH2 = /usr/lib/x86_64-linux-gnu/
-#LIBPATH2 = /flex_wrk/flexpart/lib64/gfortran/lib/
-#FFLAGS = -O3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-#FFLAGS = -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
-LDFLAGS = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-#
-
-
-MODOBJS = \
-par_mod.o com_mod.o \
-conv_mod.o hanna_mod.o \
-interpol_mod.o cmapf_mod.o \
-unc_mod.o oh_mod.o \
-xmass_mod.o flux_mod.o \
-point_mod.o outg_mod.o
-
-OBJECTS = \
-writeheader.o assignland.o\
-calcpar.o part0.o \
-caldate.o partdep.o \
-coordtrafo.o psih.o \
-raerod.o \
-drydepokernel.o random.o \
-erf.o readavailable.o \
-ew.o readcommand.o \
-advance.o readdepo.o \
-releaseparticles.o psim.o \
-FLEXPART.o readlanduse.o \
-getfields.o init_domainfill.o\
-interpol_wind.o readoutgrid.o \
-interpol_all.o readpaths.o \
-getrb.o readreceptors.o \
-getrc.o readreleases.o \
-getvdep.o readspecies.o \
-interpol_misslev.o readwind.o \
-conccalc.o richardson.o \
-concoutput.o scalev.o \
-pbl_profile.o readOHfield.o\
-juldate.o timemanager.o \
-interpol_vdep.o interpol_rain.o \
-verttransform.o partoutput.o \
-hanna.o wetdepokernel.o \
-mean.o wetdepo.o \
-hanna_short.o windalign.o \
-obukhov.o gridcheck.o \
-hanna1.o initialize.o \
- gridcheck_nests.o \
-readwind_nests.o calcpar_nests.o \
-verttransform_nests.o interpol_all_nests.o \
-interpol_wind_nests.o interpol_misslev_nests.o \
-interpol_vdep_nests.o interpol_rain_nests.o \
-getvdep_nests.o \
-readageclasses.o readpartpositions.o \
-calcfluxes.o fluxoutput.o \
-qvsat.o skplin.o \
-convmix.o calcmatrix.o \
-convect43c.o redist.o \
-sort2.o distance.o \
-centerofmass.o plumetraj.o \
-openouttraj.o calcpv.o \
-calcpv_nests.o distance2.o \
-clustering.o interpol_wind_short.o \
-interpol_wind_short_nests.o shift_field_0.o \
-shift_field.o outgrid_init.o \
-openreceptors.o boundcond_domainfill.o\
-partoutput_short.o readoutgrid_nest.o \
-outgrid_init_nest.o writeheader_nest.o \
-concoutput_nest.o wetdepokernel_nest.o \
-drydepokernel_nest.o zenithangle.o \
-ohreaction.o getvdep_nests.o \
-initial_cond_calc.o initial_cond_output.o \
-dynamic_viscosity.o get_settling.o
-
-
-$(MAIN): $(MODOBJS) $(OBJECTS)
- $(FC) *.o -o $(MAIN) $(LDFLAGS)
-
-$(OBJECTS): $(MODOBJS)
-
-%.o: %.f90
- $(FC) -c $(FFLAGS) $<
-
-clean:
- rm *.o *.mod
-
Index: /trunk/src/openouttraj.f90
===================================================================
--- /trunk/src/openouttraj.f90 (revision 19)
+++ /trunk/src/openouttraj.f90 (revision 20)
@@ -54,7 +54,7 @@
if (ldirect.eq.1) then
- write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime,'FLEXPART V8.2'
+ write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime, trim(flexversion)
else
- write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime,'FLEXPART V8.2'
+ write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime, trim(flexversion)
endif
write(unitouttraj,*) method,lsubgrid,lconvection
Index: /trunk/src/outg_mod.f90
===================================================================
--- /trunk/src/outg_mod.f90 (revision 19)
+++ /trunk/src/outg_mod.f90 (revision 20)
@@ -43,4 +43,5 @@
real,allocatable, dimension (:,:) :: wetgridsigma
real,allocatable, dimension (:) :: sparse_dump_r
+ real,allocatable, dimension (:) :: sparse_dump_u
integer,allocatable, dimension (:) :: sparse_dump_i
Index: /trunk/src/outgrid_init.f90
===================================================================
--- /trunk/src/outgrid_init.f90 (revision 19)
+++ /trunk/src/outgrid_init.f90 (revision 20)
@@ -247,4 +247,9 @@
max(numygrid,numygridn)*numzgrid),stat=stat)
if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
+ allocate(sparse_dump_u(max(numxgrid,numxgridn)* &
+ max(numygrid,numygridn)*numzgrid),stat=stat)
+ if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
+
allocate(sparse_dump_i(max(numxgrid,numxgridn)* &
max(numygrid,numygridn)*numzgrid),stat=stat)
Index: /trunk/src/par_mod.f90
===================================================================
--- /trunk/src/par_mod.f90 (revision 19)
+++ /trunk/src/par_mod.f90 (revision 20)
@@ -28,5 +28,5 @@
! 1997 *
! *
-! Last update 10 August 2000 *
+! Last update 15 August 2013 IP *
! *
!*******************************************************************************
@@ -120,10 +120,13 @@
! Maximum dimensions of the input mother grids
!*********************************************
-
- integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92
- !integer,parameter :: nxmax=361,nymax=181,nuvzmax=61,nwzmax=61,nzmax=61
+
+ integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL
+ !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
+ !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
!integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
- integer,parameter :: nxshift=359 ! for ECMWF
- !integer,parameter :: nxshift=0 ! for GFS
+ !integer,parameter :: nxmax=1201,nymax=235,nuvzmax=58,nwzmax=58,nzmax=58
+
+ integer,parameter :: nxshift=359 ! for ECMWF
+ !integer,parameter :: nxshift=0 ! for GFS or FNL (XF)
integer,parameter :: nconvlevmax = nuvzmax-1
@@ -150,7 +153,7 @@
!*********************************************
- integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
- !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151
-
+ !integer,parameter :: maxnests=0, nxmaxn=0, nymaxn=0
+ !integer,parameter :: maxnests=1,nxmaxn=251,nymaxn=151 !ECMWF
+ integer,parameter :: maxnests=1, nxmaxn=201, nymaxn=161 ! FNL XF
! maxnests maximum number of nested grids
! nxmaxn,nymaxn maximum dimension of nested wind fields in
@@ -194,6 +197,7 @@
!**************************************************
- integer,parameter :: maxpart=2000000
- integer,parameter :: maxspec=4
+ !integer,parameter :: maxpart=4000000
+ integer,parameter :: maxpart=2000
+ integer,parameter :: maxspec=6
@@ -255,6 +259,18 @@
integer,parameter :: unitlsm=1, unitsurfdata=1, unitland=1, unitwesely=1
integer,parameter :: unitOH=1
- integer,parameter :: unitdates=94, unitheader=90, unitshortpart=95
+ integer,parameter :: unitdates=94, unitheader=90,unitheader_txt=100, unitshortpart=95
integer,parameter :: unitboundcond=89
+!******************************************************
+! integer code for missing values, used in wet scavenging (PS, 2012)
+!******************************************************
+
+ ! integer icmv
+ ! parameter(icmv=-9999)
+ integer,parameter :: icmv=-9999
+
+! Parameters for testing
+!*******************************************
+! integer :: verbosity=0
+
end module par_mod
Index: /trunk/src/readavailable.f90
===================================================================
--- /trunk/src/readavailable.f90 (revision 19)
+++ /trunk/src/readavailable.f90 (revision 20)
@@ -123,4 +123,6 @@
do k=1,numbnests
+ print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
+ print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
open(unitavailab,file=path(numpath+2*(k-1)+2) &
(1:length(numpath+2*(k-1)+2)),status='old',err=998)
@@ -275,5 +277,5 @@
return
-998 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
+998 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE FILE #### '
write(*,'(a)') ' '//path(numpath+2*(k-1)+2) &
(1:length(numpath+2*(k-1)+2))
@@ -281,5 +283,5 @@
stop
-999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! AVAILABLE IILE #### '
write(*,'(a)') ' '//path(4)(1:length(4))
write(*,*) ' #### CANNOT BE OPENED #### '
Index: /trunk/src/readcommand.f90
===================================================================
--- /trunk/src/readcommand.f90 (revision 19)
+++ /trunk/src/readcommand.f90 (revision 20)
@@ -80,85 +80,171 @@
character(len=50) :: line
logical :: old
-
+ logical :: nmlout=.true. !.false.
+ integer :: readerror
+
+ namelist /command/ &
+ ldirect, &
+ ibdate,ibtime, &
+ iedate,ietime, &
+ loutstep, &
+ loutaver, &
+ loutsample, &
+ itsplit, &
+ lsynctime, &
+ ctl, &
+ ifine, &
+ iout, &
+ ipout, &
+ lsubgrid, &
+ lconvection, &
+ lagespectra, &
+ ipin, &
+ ioutputforeachrelease, &
+ iflux, &
+ mdomainfill, &
+ ind_source, &
+ ind_receptor, &
+ mquasilag, &
+ nested_output, &
+ linit_cond, &
+ surf_only
+
+ ! Presetting namelist command
+ ldirect=1
+ ibdate=20000101
+ ibtime=0
+ iedate=20000102
+ ietime=0
+ loutstep=10800
+ loutaver=10800
+ loutsample=900
+ itsplit=999999999
+ lsynctime=900
+ ctl=-5.0
+ ifine=4
+ iout=3
+ ipout=0
+ lsubgrid=1
+ lconvection=1
+ lagespectra=0
+ ipin=1
+ ioutputforeachrelease=1
+ iflux=1
+ mdomainfill=0
+ ind_source=1
+ ind_receptor=1
+ mquasilag=0
+ nested_output=0
+ linit_cond=0
+ surf_only=0
! Open the command file and read user options
- !********************************************
-
-
+ ! Namelist input first: try to read as namelist file
+ !**************************************************************************
open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
- err=999)
-
- ! Check the format of the COMMAND file (either in free format,
- ! or using formatted mask)
- ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
- !**************************************************************************
-
- call skplin(9,unitcommand)
- read (unitcommand,901) line
-901 format (a)
- if (index(line,'LDIRECT') .eq. 0) then
- old = .false.
- else
- old = .true.
- endif
- rewind(unitcommand)
-
- ! Read parameters
- !****************
-
- call skplin(7,unitcommand)
- if (old) call skplin(1,unitcommand)
-
- read(unitcommand,*) ldirect
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ibdate,ibtime
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) iedate,ietime
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) loutstep
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) loutaver
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) loutsample
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) itsplit
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) lsynctime
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ctl
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ifine
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) iout
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ipout
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) lsubgrid
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) lconvection
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) lagespectra
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ipin
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ioutputforeachrelease
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) iflux
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) mdomainfill
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ind_source
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) ind_receptor
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) mquasilag
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) nested_output
- if (old) call skplin(3,unitcommand)
- read(unitcommand,*) linit_cond
+ form='formatted',iostat=readerror)
+ ! If fail, check if file does not exist
+ if (readerror.ne.0) then
+
+ print*,'***ERROR: file COMMAND not found in '
+ print*, path(1)(1:length(1))//'COMMAND'
+ print*, 'Check your pathnames file.'
+ stop
+
+ endif
+
+ read(unitcommand,command,iostat=readerror)
close(unitcommand)
+ ! If error in namelist format, try to open with old input code
+ if (readerror.ne.0) then
+
+ open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
+ err=999)
+
+ ! Check the format of the COMMAND file (either in free format,
+ ! or using formatted mask)
+ ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+ !**************************************************************************
+
+ call skplin(9,unitcommand)
+ read (unitcommand,901) line
+ 901 format (a)
+ if (index(line,'LDIRECT') .eq. 0) then
+ old = .false.
+ else
+ old = .true.
+ endif
+ rewind(unitcommand)
+
+ ! Read parameters
+ !****************
+
+ call skplin(7,unitcommand)
+ if (old) call skplin(1,unitcommand)
+
+ read(unitcommand,*) ldirect
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ibdate,ibtime
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) iedate,ietime
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) loutstep
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) loutaver
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) loutsample
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) itsplit
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) lsynctime
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ctl
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ifine
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) iout
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ipout
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) lsubgrid
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) lconvection
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) lagespectra
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ipin
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ioutputforeachrelease
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) iflux
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) mdomainfill
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ind_source
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) ind_receptor
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) mquasilag
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) nested_output
+ if (old) call skplin(3,unitcommand)
+ read(unitcommand,*) linit_cond
+ close(unitcommand)
+
+ endif ! input format
+
+ ! write command file in namelist format to output directory if requested
+ if (nmlout.eqv..true.) then
+ !open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist.out',status='new',err=1000)
+ open(unitcommand,file=path(2)(1:length(2))//'COMMAND.namelist',err=1000)
+ write(unitcommand,nml=command)
+ close(unitcommand)
+ ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
+ ! write(unitheader,NML=COMMAND)
+ !close(unitheader)
+ endif
+
ifine=max(ifine,1)
-
! Determine how Markov chain is formulated (for w or for w/sigw)
@@ -371,5 +457,5 @@
endif
- if(lsubgrid.ne.1) then
+ if(lsubgrid.ne.1.and.verbosity.eq.0) then
write(*,*) ' ---------------- '
write(*,*) ' INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS'
@@ -505,3 +591,7 @@
stop
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND" #### '
+ write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### '
+ write(*,'(a)') path(2)(1:length(1))
+ stop
end subroutine readcommand
Index: /trunk/src/readoutgrid.f90
===================================================================
--- /trunk/src/readoutgrid.f90 (revision 19)
+++ /trunk/src/readoutgrid.f90 (revision 20)
@@ -91,5 +91,7 @@
if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
.or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
+ write(*,*) 'outlon0,outlat0:'
write(*,*) outlon0,outlat0
+ write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT ####'
Index: /trunk/src/readpartpositions.f90
===================================================================
--- /trunk/src/readpartpositions.f90 (revision 19)
+++ /trunk/src/readpartpositions.f90 (revision 20)
@@ -77,11 +77,6 @@
read(unitpartin) numpointin
- if (numpointin.ne.numpoint) then
- write(*,*) ' #### WARNING #### '
- write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT #### '
- write(*,*) ' #### AGREE WITH CURRENT SETTINGS! #### '
- write(*,*) ' #### THIS IS JUST A CHECK TO SEE IF YOU READ #### '
- write(*,*) ' #### THE CORRECT PARTICLE DUMP FILE. #### '
- endif
+ if (numpointin.ne.numpoint) goto 995
+999 continue
do i=1,numpointin
read(unitpartin)
@@ -152,4 +147,10 @@
stop
+!995 write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
+995 write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
+ write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT #### '
+ write(*,*) ' #### AGREE WITH CURRENT SETTINGS! #### '
+! stop
+ goto 999
996 write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
Index: /trunk/src/readpaths.f90
===================================================================
--- /trunk/src/readpaths.f90 (revision 19)
+++ /trunk/src/readpaths.f90 (revision 20)
@@ -20,5 +20,5 @@
!**********************************************************************
-subroutine readpaths
+subroutine readpaths !(pathfile)
!*****************************************************************************
@@ -30,4 +30,7 @@
! *
! 1 February 1994 *
+ ! last modified *
+ ! HS, 7.9.2012 *
+ ! option to give pathnames file as command line option *
! *
!*****************************************************************************
@@ -47,15 +50,30 @@
implicit none
- integer :: i
+ integer :: i
+ character(256) :: string_test
+ character(1) :: character_test
! Read the pathname information stored in unitpath
!*************************************************
-
- open(unitpath,file='pathnames',status='old',err=999)
+ open(unitpath,file=trim(pathfile),status='old',err=999)
do i=1,numpath
read(unitpath,'(a)',err=998) path(i)
length(i)=index(path(i),' ')-1
+
+
+ string_test = path(i)
+ character_test = string_test(length(i):length(i))
+ !print*, 'character_test, string_test ', character_test, string_test
+ if ((character_test .NE. '/') .AND. (i .LT. 4)) then
+ print*, 'WARNING: path not ending in /'
+ print*, path(i)
+ path(i) = string_test(1:length(i)) // '/'
+ length(i)=length(i)+1
+ print*, 'fix: padded with /'
+ print*, path(i)
+ print*, 'length(i) increased 1'
+ endif
end do
@@ -70,4 +88,5 @@
length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
end do
+ print*,length(5),length(6)
Index: /trunk/src/readreleases.f90
===================================================================
--- /trunk/src/readreleases.f90 (revision 19)
+++ /trunk/src/readreleases.f90 (revision 20)
@@ -57,4 +57,6 @@
! weta, wetb parameters to determine the wet scavenging coefficient *
! zpoint1,zpoint2 height range, over which release takes place *
+ ! num_min_discrete if less, release cannot be randomized and happens at *
+ ! time mid-point of release interval *
! *
!*****************************************************************************
@@ -68,6 +70,7 @@
integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+ integer,parameter :: num_min_discrete=100
real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
- real(kind=dp) :: jul1,jul2,juldate
+ real(kind=dp) :: jul1,jul2,julm,juldate
character(len=50) :: line
logical :: old
@@ -376,4 +379,5 @@
jul1=juldate(id1,it1)
jul2=juldate(id2,it2)
+ julm=(jul1+jul2)/2.
if (jul1.gt.jul2) then
write(*,*) 'FLEXPART MODEL ERROR'
@@ -391,6 +395,11 @@
stop
endif
- ireleasestart(numpoint)=int((jul1-bdate)*86400.)
- ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+ if (npart(numpoint).gt.num_min_discrete) then
+ ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+ ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+ else
+ ireleasestart(numpoint)=int((julm-bdate)*86400.)
+ ireleaseend(numpoint)=int((julm-bdate)*86400.)
+ endif
else if (ldirect.eq.-1) then
if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
@@ -401,6 +410,11 @@
stop
endif
- ireleasestart(numpoint)=int((jul1-bdate)*86400.)
- ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+ if (npart(numpoint).gt.num_min_discrete) then
+ ireleasestart(numpoint)=int((jul1-bdate)*86400.)
+ ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+ else
+ ireleasestart(numpoint)=int((julm-bdate)*86400.)
+ ireleaseend(numpoint)=int((julm-bdate)*86400.)
+ endif
endif
endif
Index: /trunk/src/readwind.f90
===================================================================
--- /trunk/src/readwind.f90 (revision 19)
+++ /trunk/src/readwind.f90 (revision 20)
@@ -101,4 +101,9 @@
character(len=24) :: gribErrorMsg = 'Error reading grib file'
character(len=20) :: gribFunction = 'readwind'
+
+ !HSO conversion of ECMWF etadot to etadot*dp/deta
+ logical :: etacon=.false.
+ real,parameter :: p00=101325.
+ real :: dak,dbk
hflswitch=.false.
@@ -365,4 +370,20 @@
endif
+ ! convert from ECMWF etadot to etadot*dp/deta as needed by FLEXPART
+ if(etacon.eqv..true.) then
+ do k=1,nwzmax
+ dak=akm(k+1)-akm(k)
+ dbk=bkm(k+1)-bkm(k)
+ do i=0,nxmin1
+ do j=0,nymin1
+ wwh(i,j,k)=2*wwh(i,j,k)*ps(i,j,1,n)*(dak/ps(i,j,1,n)+dbk)/(dak/p00+dbk)
+ if (k.gt.1) then
+ wwh(i,j,k)=wwh(i,j,k)-wwh(i,j,k-1)
+ endif
+ end do
+ end do
+ end do
+ endif
+
! For global fields, assign the leftmost data column also to the rightmost
! data column; if required, shift whole grid by nxshift grid points
Index: /trunk/src/readwind_gfs.f90
===================================================================
--- /trunk/src/readwind_gfs.f90 (revision 19)
+++ /trunk/src/readwind_gfs.f90 (revision 20)
@@ -703,5 +703,4 @@
endif
-
if(iumax.ne.nuvz) stop 'READWIND: NUVZ NOT CONSISTENT'
if(iumax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
Index: /trunk/src/timemanager.f90
===================================================================
--- /trunk/src/timemanager.f90 (revision 19)
+++ /trunk/src/timemanager.f90 (revision 20)
@@ -128,5 +128,13 @@
!**********************************************************************
- !write (*,*) 'starting simulation'
+
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> starting simulation'
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+ endif
+ endif
+
do itime=0,ideltas,lsynctime
@@ -142,6 +150,10 @@
!********************************************************************
- if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) &
+ if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call wetdepo'
+ endif
call wetdepo(itime,lsynctime,loutnext)
+ endif
if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
@@ -156,23 +168,56 @@
!*************************************
- if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) &
- call convmix(itime)
+ if ((ldirect.eq.-1).and.(lconvection.eq.1).and.(itime.lt.0)) then
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call convmix -- backward'
+ endif
+ call convmix(itime)
+ if (verbosity.gt.1) then
+ !CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+ endif
+ endif
! Get necessary wind fields if not available
!*******************************************
-
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call getfields'
+ endif
call getfields(itime,nstop1)
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+ endif
if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
+
! Release particles
!******************
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> Release particles'
+ endif
+
if (mdomainfill.ge.1) then
if (itime.eq.0) then
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call init_domainfill'
+ endif
call init_domainfill
else
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call boundcond_domainfill'
+ endif
call boundcond_domainfill(itime,loutend)
endif
else
+ if (verbosity.gt.0) then
+ print*,'call releaseparticles'
+ endif
call releaseparticles(itime)
+ if (verbosity.gt.1) then
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
+ endif
endif
@@ -182,7 +227,10 @@
!**************************************************************
- if ((ldirect.eq.1).and.(lconvection.eq.1)) &
- call convmix(itime)
-
+ if ((ldirect.eq.1).and.(lconvection.eq.1)) then
+ if (verbosity.gt.0) then
+ write (*,*) 'timemanager> call convmix -- forward'
+ endif
+ call convmix(itime)
+ endif
! If middle of averaging period of output fields is reached, accumulated
@@ -299,7 +347,24 @@
if ((itime.eq.loutend).and.(outnum.gt.0.)) then
if ((iout.le.3.).or.(iout.eq.5)) then
+ if (surf_only.ne.1) then
call concoutput(itime,outnum,gridtotalunc, &
wetgridtotalunc,drygridtotalunc)
- if (nested_output.eq.1) call concoutput_nest(itime,outnum)
+ else
+ if (verbosity.eq.1) then
+ print*,'call concoutput_surf '
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+ call concoutput_surf(itime,outnum,gridtotalunc, &
+ wetgridtotalunc,drygridtotalunc)
+ if (verbosity.eq.1) then
+ print*,'called concoutput_surf '
+ CALL SYSTEM_CLOCK(count_clock)
+ WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0
+ endif
+ endif
+
+ if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
+ if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
outnum=0.
endif
Index: /trunk/src/verttransform.f90
===================================================================
--- /trunk/src/verttransform.f90 (revision 19)
+++ /trunk/src/verttransform.f90 (revision 20)
@@ -49,4 +49,6 @@
! Sabine Eckhardt, March 2007
! added the variable cloud for use with scavenging - descr. in com_mod
+ ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
+ ! note that also other subroutines are affected by the fix
!*****************************************************************************
! *
@@ -70,11 +72,13 @@
integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
- integer :: rain_cloud_above,kz_inv
+ integer :: rain_cloud_above,kz_inv !SE
+ integer icloudtop !PS
real :: f_qvsat,pressure
- real :: rh,lsp,convp
+ !real :: rh,lsp,convp
+ real :: rh,lsp,convp,prec,rhmin
real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
real :: xlon,ylat,xlonr,dzdx,dzdy
- real :: dzdx1,dzdx2,dzdy1,dzdy2
+ real :: dzdx1,dzdx2,dzdy1,dzdy2, precmin
real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
@@ -83,5 +87,7 @@
real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+ logical lconvectprec
real,parameter :: const=r_air/ga
+ parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics
logical :: init = .true.
@@ -545,38 +551,108 @@
! create a cloud and rainout/washout field, clouds occur where rh>80%
! total cloudheight is stored at level 0
+
+
+
do jy=0,nymin1
do ix=0,nxmin1
- rain_cloud_above=0
- lsp=lsprec(ix,jy,1,n)
- convp=convprec(ix,jy,1,n)
- cloudsh(ix,jy,n)=0
- do kz_inv=1,nz-1
- kz=nz-kz_inv+1
- pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
- rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
- clouds(ix,jy,kz,n)=0
- if (rh.gt.0.8) then ! in cloud
- if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
- rain_cloud_above=1
- cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
- height(kz)-height(kz-1)
- if (lsp.ge.convp) then
- clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
- else
- clouds(ix,jy,kz,n)=2 ! convp dominated rainout
- endif
- else ! no precipitation
- clouds(ix,jy,kz,n)=1 ! cloud
+
+
+
+ ! rain_cloud_above=0
+ ! lsp=lsprec(ix,jy,1,n)
+ ! convp=convprec(ix,jy,1,n)
+ ! cloudsh(ix,jy,n)=0
+ ! do kz_inv=1,nz-1
+ ! kz=nz-kz_inv+1
+ ! pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+ ! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+ ! clouds(ix,jy,kz,n)=0
+ ! if (rh.gt.0.8) then ! in cloud
+ ! if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+ ! rain_cloud_above=1
+ ! cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+ ! height(kz)-height(kz-1)
+ ! if (lsp.ge.convp) then
+ ! clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+ ! else
+ ! clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+ ! endif
+ ! else ! no precipitation
+ ! clouds(ix,jy,kz,n)=1 ! cloud
+ ! endif
+ ! else ! no cloud
+ ! if (rain_cloud_above.eq.1) then ! scavenging
+ ! if (lsp.ge.convp) then
+ ! clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+ ! else
+ ! clouds(ix,jy,kz,n)=4 ! convp dominated washout
+ ! endif
+ ! endif
+ ! endif
+ ! end do
+
+
+ ! PS 3012
+
+ lsp=lsprec(ix,jy,1,n)
+ convp=convprec(ix,jy,1,n)
+ prec=lsp+convp
+ if (lsp.gt.convp) then ! prectype='lsp'
+ lconvectprec = .false.
+ else ! prectype='cp '
+ lconvectprec = .true.
+ endif
+ rhmin = 0.90 ! standard condition for presence of clouds
+!PS note that original by Sabine Eckhart was 80%
+!PS however, for T<-20 C we consider saturation over ice
+!PS so I think 90% should be enough
+ icloudbot(ix,jy,n)=icmv
+ icloudtop=icmv ! this is just a local variable
+98 do kz=1,nz
+ pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+ rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+!ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
+ if (rh .gt. rhmin) then
+ if (icloudbot(ix,jy,n) .eq. icmv) then
+ icloudbot(ix,jy,n)=nint(height(kz))
+ endif
+ icloudtop=nint(height(kz)) ! use int to save memory
endif
- else ! no cloud
- if (rain_cloud_above.eq.1) then ! scavenging
- if (lsp.ge.convp) then
- clouds(ix,jy,kz,n)=5 ! lsp dominated washout
- else
- clouds(ix,jy,kz,n)=4 ! convp dominated washout
- endif
+ enddo
+
+!PS try to get a cloud thicker than 50 m
+!PS if there is at least .01 mm/h - changed to 0.002 and put into
+!PS parameter precpmin
+ if ((icloudbot(ix,jy,n) .eq. icmv .or. &
+ icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
+ prec .gt. precmin) then
+ rhmin = rhmin - 0.05
+ if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
+ endif
+!PS implement a rough fix for badly represented convection
+!PS is based on looking at a limited set of comparison data
+ if (lconvectprec .and. icloudtop .lt. 6000 .and. &
+ prec .gt. precmin) then
+
+ if (convp .lt. 0.1) then
+ icloudbot(ix,jy,n) = 500
+ icloudtop = 8000
+ else
+ icloudbot(ix,jy,n) = 0
+ icloudtop = 10000
endif
- endif
- end do
+ endif
+ if (icloudtop .ne. icmv) then
+ icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
+ else
+ icloudthck(ix,jy,n) = icmv
+ endif
+!PS get rid of too thin clouds
+ if (icloudthck(ix,jy,n) .lt. 50) then
+ icloudbot(ix,jy,n)=icmv
+ icloudthck(ix,jy,n)=icmv
+ endif
+
+
end do
end do
@@ -605,3 +681,5 @@
!104 continue
! close(4)
+
+
end subroutine verttransform
Index: /trunk/src/wetdepo.f90
===================================================================
--- /trunk/src/wetdepo.f90 (revision 19)
+++ /trunk/src/wetdepo.f90 (revision 20)
@@ -39,4 +39,7 @@
! Code may not be correct for decay of deposition! *
! *
+ ! Modification by Sabine Eckhart to introduce a new in-/below-cloud
+ ! scheme, not dated
+ ! Petra Seibert, 2011/2012: Fixing some deficiencies in this modification
!*****************************************************************************
! *
@@ -71,9 +74,11 @@
integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
- integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
- integer :: ks, kp
+ integer :: ngrid,itage,nage,hz,il,interp_time, n , clouds_v !NIK scheme
+ integer :: kz !PS scheme
+ integer :: ks, kp, n1,n2, icbot,ictop, indcloud
+ integer :: scheme_number ! NIK==1, PS ==2
real :: S_i, act_temp, cl, cle ! in cloud scavenging
real :: clouds_h ! cloud height for the specific grid point
- real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav
+ real :: xtn,ytn,lsp,convp,cc,grfraction,prec,wetscav, precsub,f
real :: wetdeposit(maxspec),restmass
real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
@@ -146,20 +151,38 @@
interp_time=nint(itime-0.5*ltsample)
- if (ngrid.eq.0) then
- call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
- 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
- memtime(1),memtime(2),interp_time,lsp,convp,cc)
- else
- call interpol_rain_nests(lsprecn,convprecn,tccn, &
- nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
- memtime(1),memtime(2),interp_time,lsp,convp,cc)
- endif
-
- if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+! PS nest case still needs to be implemented!!
+! if (ngrid.eq.0) then
+! call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
+! 1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
+! memtime(1),memtime(2),interp_time,lsp,convp,cc)
+ call interpol_rain(lsprec,convprec,tcc, &
+ icloudbot,icloudthck,nxmax,nymax,1,nx,ny, &
+ memind,sngl(xtra1(jpart)),sngl(ytra1(jpart)),1,memtime(1), &
+ memtime(2),interp_time,lsp,convp,cc,icbot,ictop,icmv)
+! else
+! call interpol_rain_nests(lsprecn,convprecn,tccn, &
+! nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
+! memtime(1),memtime(2),interp_time,lsp,convp,cc)
+! endif
+
+
+ ! if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
+ !PS 2012: subtract a small value, eg 0.01 mm/h, to remove spurious precip
+ prec = lsp+convp
+ precsub = 0.01
+ if (prec .lt. precsub) then
+ goto 20
+ else
+ f = (prec-precsub)/prec
+ lsp = f*lsp
+ convp = f*convp
+ endif
+
! get the level were the actual particle is in
do il=2,nz
if (height(il).gt.ztra1(jpart)) then
- hz=il-1
+ !hz=il-1
+ kz=il-1
goto 26
endif
@@ -174,15 +197,31 @@
! scavenging is done
- if (ngrid.eq.0) then
- clouds_v=clouds(ix,jy,hz,n)
- clouds_h=cloudsh(ix,jy,n)
- else
- clouds_v=cloudsn(ix,jy,hz,n,ngrid)
- clouds_h=cloudsnh(ix,jy,n,ngrid)
- endif
- !write(*,*) 'there is
- ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
- if (clouds_v.le.1) goto 20
- !write (*,*) 'there is scavenging'
+!old scheme
+! if (ngrid.eq.0) then
+! clouds_v=clouds(ix,jy,hz,n)
+! clouds_h=cloudsh(ix,jy,n)
+! else
+! clouds_v=cloudsn(ix,jy,hz,n,ngrid)
+! clouds_h=cloudsnh(ix,jy,n,ngrid)
+! endif
+! !write(*,*) 'there is
+! ! + precipitation',(clouds(ix,jy,ihz,n),ihz=1,20),lsp,convp,hz
+! if (clouds_v.le.1) goto 20
+! !write (*,*) 'there is scavenging'
+
+ ! PS: part of 2011/2012 fix
+
+ if (ztra1(jpart) .le. float(ictop)) then
+ if (ztra1(jpart) .gt. float(icbot)) then
+ indcloud = 2 ! in-cloud
+ else
+ indcloud = 1 ! below-cloud
+ endif
+ elseif (ictop .eq. icmv) then
+ indcloud = 0 ! no cloud found, use old scheme
+ else
+ goto 20 ! above cloud
+ endif
+
! 1) Parameterization of the the area fraction of the grid cell where the
@@ -228,6 +267,12 @@
!**********************************************************
- do ks=1,nspec ! loop over species
+ do ks=1,nspec ! loop over species
wetdeposit(ks)=0.
+
+
+ !conflicting changes to the same routine: 1=NIK 2 =PS
+ scheme_number=2
+ if (scheme_number.eq.1) then !NIK
+
if (weta(ks).gt.0.) then
if (clouds_v.ge.4) then
@@ -236,4 +281,5 @@
wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff.
! write(*,*) 'bel. wetscav: ',wetscav
+
else ! below_cloud clouds_v is lt 4 and gt 1 -> in cloud scavenging
! IN CLOUD SCAVENGING
@@ -245,12 +291,18 @@
act_temp=tt(ix,jy,hz,n)
endif
- cl=2E-7*prec**0.36
+
+! NIK 31.01.2013: SPECIES defined parameters for the in-cloud scavening
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+! wetd_in: Scaling factor for the total in-cloud scavenging (default 1.0-no scaling)
+ cl=weta_in(ks)*prec**wetb_in(ks)
if (dquer(ks).gt.0) then ! is particle
- S_i=0.9/cl
+ S_i=wetc_in(ks)/cl
else ! is gas
cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
S_i=1/cle
endif
- wetscav=S_i*prec/3.6E6/clouds_h
+ wetscav=S_i*prec/3.6E6/clouds_h/wetd_in(ks)
! write(*,*) 'in. wetscav:'
! + ,wetscav,cle,cl,act_temp,prec,clouds_h
@@ -289,4 +341,72 @@
wetdeposit(ks)=0.
endif ! weta(k)
+
+ elseif (scheme_number.eq.2) then ! PS
+
+!PS indcloud=0 ! Use this for FOR TESTING,
+!PS will skip the new in/below cloud method !!!
+
+ if (weta(ks).gt.0.) then
+ if (indcloud .eq. 1) then ! BELOW CLOUD SCAVENGING
+!C for aerosols and not highliy soluble substances weta=5E-6
+ wetscav=weta(ks)*prec**wetb(ks) ! scavenging coeff.
+!c write(*,*) 'bel. wetscav: ',wetscav
+ elseif (indcloud .eq. 2) then ! IN CLOUD SCAVENGING
+ if (ngrid.gt.0) then
+ act_temp=ttn(ix,jy,kz,n,ngrid)
+ else
+ act_temp=tt(ix,jy,kz,n)
+ endif
+
+! from NIK
+! weta_in=2.0E-07 (default)
+! wetb_in=0.36 (default)
+! wetc_in=0.9 (default)
+
+
+ cl=2E-7*prec**0.36
+ if (dquer(ks).gt.0) then ! is particle
+ S_i=0.9/cl
+ else ! is gas
+ cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
+ S_i=1/cle
+ endif
+ wetscav=S_i*prec/3.6E6/(ictop-icbot) ! 3.6e6 converts mm/h to m/s
+ else ! PS: no cloud diagnosed, old scheme,
+!CPS using with fixed a,b for simplicity, one may wish to change!!
+ wetscav = 1.e-4*prec**0.62
+ endif
+
+
+ wetdeposit(ks)=xmass1(jpart,ks)* &
+ ! (1.-exp(-wetscav*abs(ltsample)))*fraction ! wet deposition
+ (1.-exp(-wetscav*abs(ltsample)))*grfraction ! fraction = grfraction (PS)
+ restmass = xmass1(jpart,ks)-wetdeposit(ks)
+ if (ioutputforeachrelease.eq.1) then
+ kp=npoint(jpart)
+ else
+ kp=1
+ endif
+ if (restmass .gt. smallnum) then
+ xmass1(jpart,ks)=restmass
+!cccccccccccccccc depostatistic
+!c wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
+!cccccccccccccccc depostatistic
+ else
+ xmass1(jpart,ks)=0.
+ endif
+!C Correct deposited mass to the last time step when radioactive decay of
+!C gridded deposited mass was calculated
+ if (decay(ks).gt.0.) then
+ wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
+ endif
+ else ! weta(k)<0
+ wetdeposit(ks)=0.
+ endif
+
+ endif !on scheme
+
+
+
end do
@@ -295,11 +415,21 @@
!*****************************************************************************
- if (ldirect.eq.1) then
- call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
- real(ytra1(jpart)),nage,kp)
- if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
- wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
- nage,kp)
- endif
+! if (ldirect.eq.1) then
+! call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
+! real(ytra1(jpart)),nage,kp)
+! if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
+! wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)), &
+! nage,kp)
+! endif
+
+ !PS
+ if (ldirect.eq.1) then
+ call wetdepokernel(nclass(jpart),wetdeposit, &
+ sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+ if (nested_output.eq.1) &
+ call wetdepokernel_nest(nclass(jpart),wetdeposit, &
+ sngl(xtra1(jpart)),sngl(ytra1(jpart)),nage,kp)
+ endif
+
20 continue
Index: /trunk/src/writeheader.f90
===================================================================
--- /trunk/src/writeheader.f90 (revision 19)
+++ /trunk/src/writeheader.f90 (revision 20)
@@ -67,7 +67,7 @@
if (ldirect.eq.1) then
- write(unitheader) ibdate,ibtime,'FLEXPART V9.0'
+ write(unitheader) ibdate,ibtime, trim(flexversion)
else
- write(unitheader) iedate,ietime,'FLEXPART V9.0'
+ write(unitheader) iedate,ietime, trim(flexversion)
endif
Index: /trunk/src/writeheader_nest_surf.f90
===================================================================
--- /trunk/src/writeheader_nest_surf.f90 (revision 20)
+++ /trunk/src/writeheader_nest_surf.f90 (revision 20)
@@ -0,0 +1,156 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine writeheader_nest_surf
+
+ !*****************************************************************************
+ ! *
+ ! This routine produces a file header containing basic information on the *
+ ! settings of the FLEXPART run. *
+ ! The header file is essential and must be read in by any postprocessing *
+ ! program before reading in the output data. *
+ ! *
+ ! Author: A. Stohl *
+ ! *
+ ! 7 August 2002 *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! *
+ ! xlon longitude *
+ ! xl model x coordinate *
+ ! ylat latitude *
+ ! yl model y coordinate *
+ ! *
+ !*****************************************************************************
+
+ use point_mod
+ use outg_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+ real :: xp1,yp1,xp2,yp2
+
+
+ !************************
+ ! Open header output file
+ !************************
+
+ open(unitheader,file=path(2)(1:length(2))//'header_nest_grid_time', &
+ form='unformatted',err=998)
+
+
+ ! Write the header information
+ !*****************************
+
+ if (ldirect.eq.1) then
+ write(unitheader) ibdate,ibtime,trim(flexversion)
+ else
+ write(unitheader) iedate,ietime,trim(flexversion)
+ endif
+
+ ! Write info on output interval, averaging time, sampling time
+ !*************************************************************
+
+ write(unitheader) loutstep,loutaver,loutsample
+
+ ! Write information on output grid setup
+ !***************************************
+
+ write(unitheader) outlon0n,outlat0n,numxgridn,numygridn, &
+ dxoutn,dyoutn
+ write(unitheader) 1,(outheight(1),i=1,1)
+
+ call caldate(bdate,jjjjmmdd,ihmmss)
+ write(unitheader) jjjjmmdd,ihmmss
+
+ ! Write number of species, and name for each species (+extra name for depositions)
+ ! Indicate the dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+ ! concentration fields
+ !*****************************************************************************
+
+ write(unitheader) 3*nspec,maxpointspec_act
+ do i=1,nspec
+ write(unitheader) 1,'WD_'//species(i)(1:7)
+ write(unitheader) 1,'DD_'//species(i)(1:7)
+ write(unitheader) 1,species(i)
+ end do
+
+ ! Write information on release points: total number, then for each point:
+ ! start, end, coordinates, # of particles, name, mass
+ !************************************************************************
+
+ write(unitheader) numpoint
+ do i=1,numpoint
+ write(unitheader) ireleasestart(i),ireleaseend(i),kindz(i)
+ xp1=xpoint1(i)*dx+xlon0
+ yp1=ypoint1(i)*dy+ylat0
+ xp2=xpoint2(i)*dx+xlon0
+ yp2=ypoint2(i)*dy+ylat0
+ write(unitheader) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+ write(unitheader) npart(i),1
+ if (numpoint.le.1000) then
+ write(unitheader) compoint(i)
+ else
+ write(unitheader) compoint(1001)
+ endif
+ do j=1,nspec
+ write(unitheader) xmass(i,j)
+ write(unitheader) xmass(i,j)
+ write(unitheader) xmass(i,j)
+ end do
+ end do
+
+ ! Write information on some model switches
+ !*****************************************
+
+ write(unitheader) method,lsubgrid,lconvection, &
+ ind_source,ind_receptor
+
+ ! Write age class information
+ !****************************
+
+ write(unitheader) nageclass,(lage(i),i=1,nageclass)
+
+
+ ! Write topography to output file
+ !********************************
+
+ do ix=0,numxgridn-1
+ write(unitheader) (orooutn(ix,jy),jy=0,numygridn-1)
+ end do
+ close(unitheader)
+
+ return
+
+
+998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### '
+ write(*,*) ' #### '//path(2)(1:length(2))//'header'//' #### '
+ write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### '
+ write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+ write(*,*) ' #### THE PROGRAM AGAIN. #### '
+ stop
+
+end subroutine writeheader_nest_surf
Index: /trunk/src/writeheader_surf.f90
===================================================================
--- /trunk/src/writeheader_surf.f90 (revision 20)
+++ /trunk/src/writeheader_surf.f90 (revision 20)
@@ -0,0 +1,156 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine writeheader_surf
+
+ !*****************************************************************************
+ ! *
+ ! This routine produces a file header containing basic information on the *
+ ! settings of the FLEXPART run. *
+ ! The header file is essential and must be read in by any postprocessing *
+ ! program before reading in the output data. *
+ ! *
+ ! Author: A. Stohl *
+ ! *
+ ! 7 August 2002 *
+ ! *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! *
+ ! xlon longitude *
+ ! xl model x coordinate *
+ ! ylat latitude *
+ ! yl model y coordinate *
+ ! *
+ !*****************************************************************************
+
+ use point_mod
+ use outg_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+ real :: xp1,yp1,xp2,yp2
+
+
+ !************************
+ ! Open header output file
+ !************************
+
+ open(unitheader,file=path(2)(1:length(2))//'header_grid_time', &
+ form='unformatted',err=998)
+
+
+ ! Write the header information
+ !*****************************
+
+ if (ldirect.eq.1) then
+ write(unitheader) ibdate,ibtime, trim(flexversion)
+ else
+ write(unitheader) iedate,ietime, trim(flexversion)
+ endif
+
+ ! Write info on output interval, averaging time, sampling time
+ !*************************************************************
+
+ write(unitheader) loutstep,loutaver,loutsample
+
+ ! Write information on output grid setup
+ !***************************************
+
+ write(unitheader) outlon0,outlat0,numxgrid,numygrid, &
+ dxout,dyout
+ write(unitheader) 1,(outheight(1),i=1,1)
+
+ call caldate(bdate,jjjjmmdd,ihmmss)
+ write(unitheader) jjjjmmdd,ihmmss
+
+ ! Write number of species, and name for each species (+extra name for depositions)
+ ! Indicate the dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+ ! concentration fields
+ !*****************************************************************************
+
+ write(unitheader) 3*nspec,maxpointspec_act
+ do i=1,nspec
+ write(unitheader) 1,'WD_'//species(i)(1:7)
+ write(unitheader) 1,'DD_'//species(i)(1:7)
+ write(unitheader) 1,species(i)
+ end do
+
+ ! Write information on release points: total number, then for each point:
+ ! start, end, coordinates, # of particles, name, mass
+ !************************************************************************
+
+ write(unitheader) numpoint
+ do i=1,numpoint
+ write(unitheader) ireleasestart(i),ireleaseend(i),kindz(i)
+ xp1=xpoint1(i)*dx+xlon0
+ yp1=ypoint1(i)*dy+ylat0
+ xp2=xpoint2(i)*dx+xlon0
+ yp2=ypoint2(i)*dy+ylat0
+ write(unitheader) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+ write(unitheader) npart(i),1
+ if (numpoint.le.1000) then
+ write(unitheader) compoint(i)
+ else
+ write(unitheader) compoint(1001)
+ endif
+ do j=1,nspec
+ write(unitheader) xmass(i,j)
+ write(unitheader) xmass(i,j)
+ write(unitheader) xmass(i,j)
+ end do
+ end do
+
+ ! Write information on some model switches
+ !*****************************************
+
+ write(unitheader) method,lsubgrid,lconvection, &
+ ind_source,ind_receptor
+
+ ! Write age class information
+ !****************************
+
+ write(unitheader) nageclass,(lage(i),i=1,nageclass)
+
+
+ ! Write topography to output file
+ !********************************
+
+ do ix=0,numxgrid-1
+ write(unitheader) (oroout(ix,jy),jy=0,numygrid-1)
+ end do
+ close(unitheader)
+
+ return
+
+
+998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### '
+ write(*,*) ' #### '//path(2)(1:length(2))//'header'//' #### '
+ write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### '
+ write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+ write(*,*) ' #### THE PROGRAM AGAIN. #### '
+ stop
+
+end subroutine writeheader_surf
Index: /trunk/src/writeheader_txt.f90
===================================================================
--- /trunk/src/writeheader_txt.f90 (revision 20)
+++ /trunk/src/writeheader_txt.f90 (revision 20)
@@ -0,0 +1,190 @@
+!**********************************************************************
+! 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 . *
+!**********************************************************************
+
+subroutine writeheader_txt
+
+ !*****************************************************************************
+ ! *
+ ! This routine produces a file header containing basic information on the *
+ ! settings of the FLEXPART run. *
+ ! The header file is essential and must be read in by any postprocessing *
+ ! program before reading in the output data. *
+ ! *
+ ! Author: A. Stohl *
+ ! *
+ ! 7 August 2002 *
+ ! modified IP 2013 for text output *
+ !*****************************************************************************
+ ! *
+ ! Variables: *
+ ! *
+ ! xlon longitude *
+ ! xl model x coordinate *
+ ! ylat latitude *
+ ! yl model y coordinate *
+ ! *
+ !*****************************************************************************
+
+ use point_mod
+ use outg_mod
+ use par_mod
+ use com_mod
+
+ implicit none
+
+ integer :: jjjjmmdd,ihmmss,i,ix,jy,j
+ real :: xp1,yp1,xp2,yp2
+
+
+ !************************
+ ! Open header output file
+ !************************
+
+ open(unitheader,file=path(2)(1:length(2))//'header_txt', &
+ form='formatted',err=998)
+ open(unitheader_txt,file=path(2)(1:length(2))//'header_txt_releases', &
+ form='formatted',err=998)
+
+
+ ! Write the header information
+ !*****************************
+
+ write(unitheader,*) '# ibdate,ibtime, iedate, ietime, flexversion'
+ write(unitheader,*) ibdate, ibtime, iedate, ietime, trim(flexversion) ! 'FLEXPART V9.0'
+ !if (ldirect.eq.1) then
+ ! write(unitheader,*) ibdate,ibtime,trim(flexversion) ! 'FLEXPART V9.0'
+ !else
+ ! write(unitheader,*) iedate,ietime,trim(flexversion) ! 'FLEXPART V9.0'
+ !endif
+
+ ! Write info on output interval, averaging time, sampling time
+ !*************************************************************
+
+ write(unitheader,*) '# interval, averaging time, sampling time'
+ write(unitheader,*) loutstep,loutaver,loutsample
+
+ ! Write information on output grid setup
+ !***************************************
+
+ write(unitheader,*) '# information on grid setup '
+ write(unitheader,*) '#outlon0,outlat0,numxgrid,numygrid,dxout,dyout'
+ write(unitheader,*) outlon0,outlat0,numxgrid,numygrid, &
+ dxout,dyout
+ write(unitheader,*) '# numzgrid, outheight(.) '
+ write(unitheader,*) numzgrid,(outheight(i),i=1,numzgrid)
+
+ write(unitheader,*) '# jjjjmmdd,ihmmss'
+ call caldate(bdate,jjjjmmdd,ihmmss)
+ write(unitheader,*) jjjjmmdd,ihmmss
+
+ ! Write number of species, and name for each species (+extra name for depositions)
+ ! Indicate the vertical dimension of the fields (i.e., 1 for deposition fields, numzgrid for
+ ! concentration fields
+ !*****************************************************************************
+
+ write(unitheader,*) '# information on species'
+ write(unitheader,*) '# 3*nspec,maxpointspec_act'
+ write(unitheader,*) 3*nspec,maxpointspec_act
+ write(unitheader,*) '# for nspec:'
+ write(unitheader,*) '# 1, WD_ '
+ write(unitheader,*) '# 1, DD_ '
+ write(unitheader,*) '# numzgrid,species'
+ do i=1,nspec
+ write(unitheader,*) 1,'WD_'//species(i)(1:7)
+ write(unitheader,*) 1,'DD_'//species(i)(1:7)
+ write(unitheader,*) numzgrid,species(i)
+ end do
+
+ ! Write information on release points: total number, then for each point:
+ ! start, end, coordinates, # of particles, name, mass
+ !************************************************************************
+
+
+ write(unitheader_txt,*) '# information on release points'
+ write(unitheader_txt,*) '# numpoint'
+ write(unitheader_txt,*) numpoint
+ write(unitheader_txt,*) '# for numpoint:'
+ do i=1,numpoint
+ write(unitheader_txt,*) ireleasestart(i),ireleaseend(i),kindz(i)
+ xp1=xpoint1(i)*dx+xlon0
+ yp1=ypoint1(i)*dy+ylat0
+ xp2=xpoint2(i)*dx+xlon0
+ yp2=ypoint2(i)*dy+ylat0
+ write(unitheader_txt,*) xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i)
+ write(unitheader_txt,*) npart(i),1
+ if (numpoint.le.1000) then
+ write(unitheader_txt,*) compoint(i)
+ else
+ write(unitheader_txt,*) compoint(1001)
+ endif
+ do j=1,nspec
+ write(unitheader_txt,*) xmass(i,j)
+ write(unitheader_txt,*) xmass(i,j)
+ write(unitheader_txt,*) xmass(i,j)
+ end do
+ end do
+
+ ! Write information on model switches
+ !*****************************************
+
+ write(unitheader,*) '# information on model switches'
+ write(unitheader,*) '# method,lsubgrid,lconvection, ind_source,ind_receptor'
+ write(unitheader,*) method,lsubgrid,lconvection, &
+ ind_source,ind_receptor
+
+ ! Write age class information
+ !****************************
+
+ write(unitheader,*) '# information on age class '
+ write(unitheader,*) nageclass,(lage(i),i=1,nageclass)
+
+
+ !Do not write topography to text output file. Keep it on the binary one
+ !********************************
+
+ !do ix=0,numxgrid-1
+ ! write(unitheader,*) (oroout(ix,jy),jy=0,numygrid-1)
+ !end do
+
+
+
+
+
+ close(unitheader)
+ close(unitheader_txt)
+
+
+! open(unitheader,file=path(2)(1:length(2))//'header_nml', &
+! form='formatted',err=998)
+! write(unitheader,NML=COMMAND)
+! close(unitheader)
+
+ return
+
+
+998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### '
+ write(*,*) ' #### '//path(2)(1:length(2))//'header_txt'//' #### '
+ write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### '
+ write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
+ write(*,*) ' #### THE PROGRAM AGAIN. #### '
+ stop
+
+end subroutine writeheader_txt