Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
flexpart
flexpart
Commits
37d62221
Commit
37d62221
authored
Apr 13, 2018
by
Espen Sollum
Browse files
Merge branch 'dev' into release-10
parents
9f59af6c
20963b15
Changes
12
Hide whitespace changes
Inline
Side-by-side
src/FLEXPART.f90
View file @
37d62221
...
...
@@ -51,10 +51,14 @@ program flexpart
use
par_mod
use
com_mod
use
conv_mod
use
netcdf_output_mod
,
only
:
writeheader_netcdf
use
random_mod
,
only
:
gasdev1
use
class_gribfile
#ifdef USE_NCF
use
netcdf_output_mod
,
only
:
writeheader_netcdf
#endif
implicit
none
integer
::
i
,
j
,
ix
,
jy
,
inest
...
...
@@ -351,6 +355,7 @@ program flexpart
! and open files that are to be kept open throughout the simulation
!******************************************************************
#ifdef USE_NCF
if
(
lnetcdfout
.eq.
1
)
then
call
writeheader_netcdf
(
lnest
=
.false.
)
else
...
...
@@ -364,6 +369,7 @@ program flexpart
call
writeheader_nest
endif
endif
#endif
if
(
verbosity
.gt.
0
)
then
print
*
,
'call writeheader'
...
...
src/FLEXPART_MPI.f90
View file @
37d62221
...
...
@@ -52,10 +52,13 @@ program flexpart
use
com_mod
use
conv_mod
use
mpi_mod
use
netcdf_output_mod
,
only
:
writeheader_netcdf
use
random_mod
,
only
:
gasdev1
use
class_gribfile
#ifdef USE_NCF
use
netcdf_output_mod
,
only
:
writeheader_netcdf
#endif
implicit
none
integer
::
i
,
j
,
ix
,
jy
,
inest
...
...
@@ -63,7 +66,8 @@ program flexpart
character
(
len
=
256
)
::
inline_options
!pathfile, flexversion, arg2
integer
::
metdata_format
=
GRIBFILE_CENTRE_UNKNOWN
integer
::
detectformat
integer
(
selected_int_kind
(
16
)),
dimension
(
maxspec
)
::
tot_b
=
0
,
&
&
tot_i
=
0
! Initialize mpi
...
...
@@ -202,11 +206,11 @@ program flexpart
metdata_format
=
detectformat
()
if
(
metdata_format
.eq.
GRIBFILE_CENTRE_ECMWF
)
then
print
*
,
'ECMWF metdata detected'
if
(
lroot
)
print
*
,
'ECMWF metdata detected'
elseif
(
metdata_format
.eq.
GRIBFILE_CENTRE_NCEP
)
then
print
*
,
'NCEP metdata detected'
if
(
lroot
)
print
*
,
'NCEP metdata detected'
else
print
*
,
'Unknown metdata format'
if
(
lroot
)
print
*
,
'Unknown metdata format'
stop
endif
...
...
@@ -377,23 +381,24 @@ program flexpart
!******************************************************************
if
(
mp_measure_time
)
call
mpif_mtime
(
'iotime'
,
0
)
if
(
lroot
)
then
! MPI: this part root process only
if
(
lnetcdfout
.eq.
1
)
then
call
writeheader_netcdf
(
lnest
=
.false.
)
else
call
writeheader
end
if
if
(
nested_output
.eq.
1
)
then
if
(
lnetcdfout
.eq.
1
)
then
call
writeheader_netcdf
(
lnest
=
.true.
)
else
call
writeheader_nest
if
(
lroot
)
then
! MPI: this part root process only
#ifdef USE_NCF
if
(
lnetcdfout
.eq.
1
)
then
call
writeheader_netcdf
(
lnest
=
.false.
)
else
call
writeheader
end
if
if
(
nested_output
.eq.
1
)
then
if
(
lnetcdfout
.eq.
1
)
then
call
writeheader_netcdf
(
lnest
=
.true.
)
else
call
writeheader_nest
endif
endif
endif
#
endif
!
if
(
verbosity
.gt.
0
)
then
print
*
,
'call writeheader'
endif
...
...
@@ -401,7 +406,7 @@ program flexpart
call
writeheader
! 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
...
...
@@ -409,8 +414,6 @@ program flexpart
if
(
mp_measure_time
)
call
mpif_mtime
(
'iotime'
,
0
)
!open(unitdates,file=path(2)(1:length(2))//'dates')
if
(
verbosity
.gt.
0
.and.
lroot
)
then
print
*
,
'call openreceptors'
endif
...
...
@@ -480,28 +483,24 @@ program flexpart
! NIK 16.02.2005
if
(
lroot
)
then
call
MPI_Reduce
(
MPI_IN_PLACE
,
tot_blc_count
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
if
(
mp_partgroup_pid
.ge.
0
)
then
! Skip for readwind process
call
MPI_Reduce
(
tot_blc_count
,
tot_b
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
&
mp_comm_used
,
mp_ierr
)
call
MPI_Reduce
(
MPI_IN_PLACE
,
tot_inc_count
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
call
MPI_Reduce
(
tot_inc_count
,
tot_i
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
&
mp_comm_used
,
mp_ierr
)
else
if
(
mp_partgroup_pid
.ge.
0
)
then
! Skip for readwind process
call
MPI_Reduce
(
tot_blc_count
,
0
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
&
mp_comm_used
,
mp_ierr
)
call
MPI_Reduce
(
tot_inc_count
,
0
,
nspec
,
MPI_INTEGER8
,
MPI_SUM
,
id_root
,
&
&
mp_comm_used
,
mp_ierr
)
end
if
end
if
if
(
lroot
)
then
do
i
=
1
,
nspec
write
(
*
,
*
)
'**********************************************'
write
(
*
,
*
)
'Scavenging statistics for species '
,
species
(
i
),
':'
write
(
*
,
*
)
'Total number of occurences of below-cloud scavenging'
,
&
&
tot_blc_count
(
i
)
&
tot_b
(
i
)
! & tot_blc_count(i)
write
(
*
,
*
)
'Total number of occurences of in-cloud scavenging'
,
&
&
tot_inc_count
(
i
)
&
tot_i
(
i
)
! & tot_inc_count(i)
write
(
*
,
*
)
'**********************************************'
end
do
...
...
src/concoutput.f90
View file @
37d62221
...
...
@@ -139,6 +139,8 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
write
(
unitdates
,
'(a)'
)
adate
//
atime
close
(
unitdates
)
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
IF
(
init
)
THEN
file_stat
=
'OLD'
init
=
.false.
...
...
@@ -312,7 +314,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
drygridsigma
(
ix
,
jy
)
=
&
drygridsigma
(
ix
,
jy
)
*
&
sqrt
(
real
(
nclassunc
))
125
drygridsigmatotal
=
drygridsigmatotal
+
&
drygridsigmatotal
=
drygridsigmatotal
+
&
drygridsigma
(
ix
,
jy
)
endif
...
...
src/concoutput_mpi.f90
View file @
37d62221
...
...
@@ -103,7 +103,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
real
,
parameter
::
smallnum
=
tiny
(
0.0
)
! smallest number that can be handled
real
,
parameter
::
weightair
=
28.97
logical
::
sp_zer
LOGICAL
,
save
::
init
=
.true.
logical
,
save
::
init
=
.true.
character
::
adate
*
8
,
atime
*
6
character
(
len
=
3
)
::
anspec
integer
::
mind
...
...
@@ -112,7 +112,6 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
logical
::
ldates_file
integer
::
ierr
character
(
LEN
=
100
)
::
dates_char
! character :: dates_char
! Measure execution time
if
(
mp_measure_time
)
call
mpif_mtime
(
'rootonly'
,
0
)
...
...
@@ -128,7 +127,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Overwrite existing dates file on first call, later append to it
! This fixes a bug where the dates file kept growing across multiple runs
! If 'dates' file exists, make a backup
! If 'dates' file exists
in output directory
, make a backup
inquire
(
file
=
path
(
2
)(
1
:
length
(
2
))//
'dates'
,
exist
=
ldates_file
)
if
(
ldates_file
.and.
init
)
then
open
(
unit
=
unitdates
,
file
=
path
(
2
)(
1
:
length
(
2
))//
'dates'
,
form
=
'formatted'
,
&
...
...
@@ -257,23 +256,34 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
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
//
&
if
(
DRYBKDEP
.or.
WETBKDEP
)
then
!scavdep output
if
(
DRYBKDEP
)
&
open
(
unitoutgrid
,
file
=
path
(
2
)(
1
:
length
(
2
))//
'grid_drydep_'
//
adate
//
&
atime
//
'_'
//
anspec
,
form
=
'unformatted'
)
if
(
WETBKDEP
)
&
open
(
unitoutgrid
,
file
=
path
(
2
)(
1
:
length
(
2
))//
'grid_wetdep_'
//
adate
//
&
atime
//
'_'
//
anspec
,
form
=
'unformatted'
)
write
(
unitoutgrid
)
itime
else
if
((
iout
.eq.
1
)
.or.
(
iout
.eq.
3
)
.or.
(
iout
.eq.
5
))
then
if
(
ldirect
.eq.
1
)
then
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
write
(
unitoutgridppt
)
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
endif
! if deposition output
do
kp
=
1
,
maxpointspec_act
do
nage
=
1
,
nageclass
...
...
@@ -353,7 +363,7 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! Concentration output
!*********************
if
((
iout
.eq.
1
)
.or.
(
iout
.eq.
3
)
.or.
(
iout
.eq.
5
))
then
if
((
iout
.eq.
1
)
.or.
(
iout
.eq.
3
)
.or.
(
iout
.eq.
5
)
.or.
(
iout
.eq.
6
)
)
then
! Wet deposition
sp_count_i
=
0
...
...
@@ -448,10 +458,18 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
sp_fact
=
sp_fact
*
(
-1.
)
endif
sp_count_r
=
sp_count_r
+1
if
(
lparticlecountoutput
)
then
sparse_dump_r
(
sp_count_r
)
=
&
sp_fact
*
&
grid
(
ix
,
jy
,
kz
)
else
sparse_dump_r
(
sp_count_r
)
=
&
sp_fact
*
&
grid
(
ix
,
jy
,
kz
)
*
&
factor3d
(
ix
,
jy
,
kz
)/
tot_mu
(
ks
,
kp
)
end
if
! 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)=
...
...
@@ -637,24 +655,8 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
! 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
creceptor
(:,:)
=
0.
gridunc
(:,:,:,:,:,:,:)
=
0.
if
(
mp_measure_time
)
call
mpif_mtime
(
'rootonly'
,
1
)
...
...
src/getfields_mpi.f90
View file @
37d62221
...
...
@@ -232,6 +232,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end
do
40
indmin
=
indj
if
(
WETBKDEP
.and.
(
lmpreader
.or.
(
.not.
lmp_use_reader
.and.
lroot
)))
then
call
writeprecip
(
itime
,
memind
(
1
))
endif
else
! No wind fields, which can be used, are currently in memory
...
...
@@ -286,7 +290,10 @@ subroutine getfields(itime,nstop,memstat,metdata_format)
end
do
60
indmin
=
indj
mind3
=
memstat
! if (WETBKDEP.and.lroot) then
if
(
WETBKDEP
.and.
(
lmpreader
.or.
(
.not.
lmp_use_reader
.and.
lroot
)))
then
call
writeprecip
(
itime
,
memind
(
1
))
endif
endif
...
...
src/makefile
View file @
37d62221
...
...
@@ -30,6 +30,10 @@ SHELL = /bin/bash
# Compile for debugging parallel FLEXPART
# make [-j] mpi-dbg
#
# NETCDF OUTPUT
# To add support for output in netCDF format, append `ncf=yes` to the
# `make` command
#
################################################################################
## PROGRAMS
...
...
@@ -56,8 +60,7 @@ ifeq ($(gcc), 4.9)
INCPATH1
=
${ROOT_DIR}
/gcc-4.9.1/include
INCPATH2
=
${ROOT_DIR}
/include
LIBPATH1
=
${ROOT_DIR}
/lib
else
#ifeq ($(gcc), 5.4)
else
# Compiled libraries under user ~flexpart, gfortran v5.4
ROOT_DIR
=
/homevip/flexpart/
...
...
@@ -67,18 +70,18 @@ else #ifeq ($(gcc), 5.4)
INCPATH1
=
${ROOT_DIR}
/gcc-5.4.0/include
INCPATH2
=
/usr/include
LIBPATH1
=
${ROOT_DIR}
/gcc-5.4.0/lib
endif
#else
# Default: System libraries at NILU, gfortran v4.6
# F90 = /usr/bin/gfortran
# MPIF90 = /usr/bin/mpif90.openmpi
# INCPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/include
# INCPATH2 = /usr/include
# LIBPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/lib
### Enable netCDF output?
ifeq
($(ncf), yes)
NCOPT
=
-DUSE_NCF
-lnetcdff
else
NCOPT
=
-UUSE_NCF
endif
# path to gributils used to detect meteodata format
VPATH
=
gributils/
...
...
@@ -88,11 +91,12 @@ O_LEV = 2 # [0,1,2,3,g,s,fast]
O_LEV_DBG
=
g
# [0,g]
## LIBRARIES
LIBS
=
-lgrib_api_f90
-lgrib_api
-lm
-ljasper
-lnetcdff
# -fopenmp
#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
LIBS
=
-lgrib_api_f90
-lgrib_api
-lm
-ljasper
$(NCOPT)
FFLAGS
=
-I
$(INCPATH1)
-I
$(INCPATH2)
-O
$(O_LEV)
-g
-cpp
-m64
-mcmodel
=
medium
-fconvert
=
little-endian
-frecord-marker
=
4
-fmessage-length
=
0
-flto
=
jobserver
-O
$(O_LEV)
$(FUSER)
#-Warray-bounds -fcheck=all # -march=native
FFLAGS
=
-I
$(INCPATH1)
-I
$(INCPATH2)
-O
$(O_LEV)
-g
-cpp
-m64
-mcmodel
=
medium
-fconvert
=
little-endian
-frecord-marker
=
4
-fmessage-length
=
0
-flto
=
jobserver
-O
$(O_LEV)
$(NCOPT)
$(FUSER)
#-Warray-bounds -fcheck=all # -march=native
DBGFLAGS
=
-I
$(INCPATH1)
-I
$(INCPATH2)
-O
$(O_LEV_DBG)
-g3
-ggdb3
-cpp
-m64
-mcmodel
=
medium
-fconvert
=
little-endian
-frecord-marker
=
4
-fmessage-length
=
0
-flto
=
jobserver
-O
$(O_LEV_DBG)
-fbacktrace
-Wall
-fdump-core
$(FUSER)
# -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
DBGFLAGS
=
-I
$(INCPATH1)
-I
$(INCPATH2)
-O
$(O_LEV_DBG)
-g3
-ggdb3
-cpp
-m64
-mcmodel
=
medium
-fconvert
=
little-endian
-frecord-marker
=
4
-fmessage-length
=
0
-flto
=
jobserver
-O
$(O_LEV_DBG)
$(NCOPT)
-fbacktrace
-Wall
-fdump-core
$(FUSER)
# -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
LDFLAGS
=
$(FFLAGS)
-L
$(LIBPATH1)
-Wl
,-rpath,
$(LIBPATH1)
$(LIBS)
#-L
$(LIBPATH2)
LDDEBUG
=
$(DBGFLAGS)
-L
$(LIBPATH1)
$(LIBS)
#-L
$(LIBPATH2)
...
...
@@ -139,6 +143,8 @@ OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
getfields_mpi.o
\
readwind_ecmwf_mpi.o
OBJECTS_NCF
=
netcdf_output_mod.o
OBJECTS
=
\
advance.o initialize.o
\
writeheader.o writeheader_txt.o
\
...
...
@@ -195,7 +201,11 @@ ohreaction.o getvdep_nests.o \
initial_cond_calc.o initial_cond_output.o
\
dynamic_viscosity.o get_settling.o
\
initialize_cbl_vel.o re_initialize_particle.o
\
cbl.o netcdf_output_mod.o
cbl.o
ifeq
($(ncf), yes)
OBJECTS
:=
$(OBJECTS)
$(OBJECTS_NCF)
endif
%.o
:
%.mod
...
...
src/readOHfield.f90
View file @
37d62221
...
...
@@ -21,25 +21,27 @@
subroutine
readOHfield
!*****************************************************************************
! *
! Reads the OH field into memory *
! *
! AUTHOR: R.L. Thompson, Nov 2014 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonOH(nxOH) longitude of OH fields *
! latOH(nyOH) latitude of OH fields *
! altOH(nzOH) altitude of OH fields *
! etaOH(nzOH) eta-levels of OH fields *
! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3) *
! *
! *
!*****************************************************************************
!*****************************************************************************
! *
! Reads the OH field into memory *
! *
! AUTHOR: R.L. Thompson, Nov 2014 *
! *
! UPDATES: *
! 03/2018 SEC: Converted original netCDF files to binary format *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonOH(nxOH) longitude of OH fields *
! latOH(nyOH) latitude of OH fields *
! altOH(nzOH) altitude of OH fields *
! etaOH(nzOH) eta-levels of OH fields *
! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3) *
! *
! *
!*****************************************************************************
use
oh_mod
use
par_mod
...
...
@@ -47,11 +49,7 @@ subroutine readOHfield
implicit
none
include
'netcdf.inc'
character
(
len
=
150
)
::
thefile
character
(
len
=
2
)
::
mm
integer
::
nid
,
ierr
,
xid
,
yid
,
zid
,
vid
,
m
integer
::
i
,
j
,
k
,
l
,
ierr
real
,
dimension
(:),
allocatable
::
etaOH
! real, parameter :: gasct=8.314 ! gas constant
...
...
@@ -60,144 +58,39 @@ subroutine readOHfield
! real, parameter :: lrate=0.0065 ! K m-1
real
,
parameter
::
scalehgt
=
7000.
! scale height in metres
! Read OH fields and level heights
!********************************
do
m
=
1
,
12
! open netcdf file
write
(
mm
,
fmt
=
'(i2.2)'
)
m
! thefile=trim(path(1))//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
thefile
=
trim
(
ohfields_path
)//
'OH_FIELDS/'
//
'geos-chem.OH.2005'
//
mm
//
'01.nc'
ierr
=
nf_open
(
trim
(
thefile
),
NF_NOWRITE
,
nid
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
'The OH field at: '
//
thefile
//
' could not be opened'
write
(
*
,
*
)
'please copy the OH fields there, or change the path in the'
write
(
*
,
*
)
'COMMAND namelist!'
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
! inquire about variables
ierr
=
nf_inq_dimid
(
nid
,
'Lon-000'
,
xid
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_dimid
(
nid
,
'Lat-000'
,
yid
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_dimid
(
nid
,
'Alt-000'
,
zid
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
if
(
m
.eq.
1
)
then
! read dimension sizes
ierr
=
nf_inq_dimlen
(
nid
,
xid
,
nxOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_dimlen
(
nid
,
yid
,
nyOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_dimlen
(
nid
,
zid
,
nzOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
! allocate variables
allocate
(
lonOH
(
nxOH
))
allocate
(
latOH
(
nyOH
))
allocate
(
etaOH
(
nzOH
))
allocate
(
altOH
(
nzOH
))
allocate
(
OH_field
(
nxOH
,
nyOH
,
nzOH
,
12
))
allocate
(
OH_hourly
(
nxOH
,
nyOH
,
nzOH
,
2
))
! read dimension variables
ierr
=
nf_inq_varid
(
nid
,
'LON'
,
xid
)
ierr
=
nf_get_var_real
(
nid
,
xid
,
lonOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_varid
(
nid
,
'LAT'
,
yid
)
ierr
=
nf_get_var_real
(
nid
,
yid
,
latOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_varid
(
nid
,
'ETAC'
,
zid
)
ierr
=
nf_get_var_real
(
nid
,
zid
,
etaOH
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
! convert eta-level to altitude (assume surface pressure of 1010 hPa)
altOH
=
log
(
1010.
/(
etaOH
*
1010.
))
*
scalehgt
endif
! m.eq.1
! read OH_field
ierr
=
nf_inq_varid
(
nid
,
'CHEM-L_S__OH'
,
vid
)
ierr
=
nf_get_var_real
(
nid
,
vid
,
OH_field
(:,:,:,
m
))
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_close
(
nid
)
end
do
deallocate
(
etaOH
)
! Read J(O1D) photolysis rates
!********************************
open
(
unitOH
,
file
=
trim
(
ohfields_path
)
&
//
'OH_FIELDS/OH_variables.bin'
,
status
=
'old'
,
&
form
=
'UNFORMATTED'
,
iostat
=
ierr
,
convert
=
'little_endian'
)
! open netcdf file
! thefile=trim(path(1))//'OH_FIELDS/jrate_average.nc'
thefile
=
trim
(
ohfields_path
)//
'OH_FIELDS/jrate_average.nc'
ierr
=
nf_open
(
trim
(
thefile
),
NF_NOWRITE
,
nid
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
write
(
*
,
*
)
'Cannot read binary OH fields in '
,
trim
(
ohfields_path
)//
'OH_FIELDS/OH_variables.bin'
stop
endif
! read dimension variables
ierr
=
nf_inq_varid
(
nid
,
'longitude'
,
xid
)
ierr
=
nf_get_var_real
(
nid
,
xid
,
lonjr
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
ierr
=
nf_inq_varid
(
nid
,
'latitude'
,
yid
)
ierr
=
nf_get_var_real
(
nid
,
yid
,
latjr
)
if
(
ierr
.ne.
0
)
then
write
(
*
,
*
)
nf_strerror
(
ierr
)
stop
endif
! read jrate_average
ierr
=
nf_inq_varid
(
nid
,
'jrate'
,
vid
)