Changes in / [9f59af6:37d6222] in flexpart.git


Ignore:
Location:
src
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • src/FLEXPART.f90

    r0c00f1f ra9cf4b1  
    5252  use com_mod
    5353  use conv_mod
    54   use netcdf_output_mod, only: writeheader_netcdf
     54
    5555  use random_mod, only: gasdev1
    5656  use class_gribfile
     57
     58#ifdef USE_NCF
     59  use netcdf_output_mod, only: writeheader_netcdf
     60#endif
    5761
    5862  implicit none
    … …  
    352356  !******************************************************************
    353357
     358#ifdef USE_NCF
    354359  if (lnetcdfout.eq.1) then
    355360    call writeheader_netcdf(lnest=.false.)
    … …  
    365370    endif
    366371  endif
     372#endif
    367373
    368374  if (verbosity.gt.0) then

  • src/FLEXPART_MPI.f90

    r0c00f1f r20963b1  
    5353  use conv_mod
    5454  use mpi_mod
    55   use netcdf_output_mod, only: writeheader_netcdf
    5655  use random_mod, only: gasdev1
    5756  use class_gribfile
     57
     58#ifdef USE_NCF
     59  use netcdf_output_mod, only: writeheader_netcdf
     60#endif
    5861
    5962  implicit none
    … …  
    6467  integer :: metdata_format = GRIBFILE_CENTRE_UNKNOWN
    6568  integer :: detectformat
    66 
     69  integer(selected_int_kind(16)), dimension(maxspec) :: tot_b=0, &
     70       & tot_i=0
    6771
    6872
    … …  
    203207
    204208  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
    205     print *,'ECMWF metdata detected'
     209    if (lroot) print *,'ECMWF metdata detected'
    206210  elseif (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
    207     print *,'NCEP metdata detected'
     211    if (lroot) print *,'NCEP metdata detected'
    208212  else
    209     print *,'Unknown metdata format'
     213    if (lroot) print *,'Unknown metdata format'
    210214    stop
    211215  endif
    … …  
    378382
    379383  if (mp_measure_time) call mpif_mtime('iotime',0)
     384
    380385  if (lroot) then ! MPI: this part root process only
    381 
    382   if (lnetcdfout.eq.1) then
    383     call writeheader_netcdf(lnest=.false.)
    384   else
    385     call writeheader
    386   end if
    387 
    388   if (nested_output.eq.1) then
    389     if (lnetcdfout.eq.1) then
    390       call writeheader_netcdf(lnest=.true.)
    391     else
    392       call writeheader_nest
    393     endif
    394   endif
    395 
    396 !
     386#ifdef USE_NCF
     387    if (lnetcdfout.eq.1) then
     388      call writeheader_netcdf(lnest=.false.)
     389    else
     390      call writeheader
     391    end if
     392   
     393    if (nested_output.eq.1) then
     394      if (lnetcdfout.eq.1) then
     395        call writeheader_netcdf(lnest=.true.)
     396      else
     397        call writeheader_nest
     398      endif
     399    endif
     400#endif
     401
    397402    if (verbosity.gt.0) then
    398403      print*,'call writeheader'
    … …  
    402407! FLEXPART 9.2 ticket ?? write header in ASCII format
    403408    call writeheader_txt
    404 !if (nested_output.eq.1) call writeheader_nest
     409
    405410    if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
    406411    if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
    … …  
    409414
    410415  if (mp_measure_time) call mpif_mtime('iotime',0)
    411 
    412   !open(unitdates,file=path(2)(1:length(2))//'dates')
    413416
    414417  if (verbosity.gt.0 .and. lroot) then
    … …  
    481484
    482485! NIK 16.02.2005
    483   if (lroot) then
    484     call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
     486  if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
     487    call MPI_Reduce(tot_blc_count, tot_b, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    485488         & mp_comm_used, mp_ierr)
    486     call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
     489    call MPI_Reduce(tot_inc_count, tot_i, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    487490         & mp_comm_used, mp_ierr)
    488   else
    489     if (mp_partgroup_pid.ge.0) then ! Skip for readwind process
    490       call MPI_Reduce(tot_blc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    491            & mp_comm_used, mp_ierr)
    492       call MPI_Reduce(tot_inc_count, 0, nspec, MPI_INTEGER8, MPI_SUM, id_root, &
    493            & mp_comm_used, mp_ierr)
    494     end if
    495491  end if
     492 
    496493
    497494  if (lroot) then
    … …  
    500497      write(*,*) 'Scavenging statistics for species ', species(i), ':'
    501498      write(*,*) 'Total number of occurences of below-cloud scavenging', &
    502            & tot_blc_count(i)
     499           & tot_b(i)
     500!           & tot_blc_count(i)
    503501      write(*,*) 'Total number of occurences of in-cloud    scavenging', &
    504            & tot_inc_count(i)
     502           & tot_i(i)
     503!           & tot_inc_count(i)
    505504      write(*,*) '**********************************************'
    506505    end do

  • src/concoutput.f90

    r08a38b5 r20963b1  
    140140  close(unitdates) 
    141141
     142  ! Overwrite existing dates file on first call, later append to it
     143  ! This fixes a bug where the dates file kept growing across multiple runs
    142144  IF (init) THEN
    143145    file_stat='OLD'
    … …  
    313315                   drygridsigma(ix,jy)* &
    314316                   sqrt(real(nclassunc))
    315 125           drygridsigmatotal=drygridsigmatotal+ &
     317              drygridsigmatotal=drygridsigmatotal+ &
    316318                   drygridsigma(ix,jy)
    317319            endif

  • src/concoutput_mpi.f90

    r6a678e3 r20963b1  
    104104  real,parameter :: weightair=28.97
    105105  logical :: sp_zer
    106   LOGICAL,save :: init=.true.
     106  logical,save :: init=.true.
    107107  character :: adate*8,atime*6
    108108  character(len=3) :: anspec
    … …  
    113113  integer :: ierr
    114114  character(LEN=100) :: dates_char
    115 !  character :: dates_char
    116115
    117116! Measure execution time
    … …  
    129128! This fixes a bug where the dates file kept growing across multiple runs
    130129
    131 ! If 'dates' file exists, make a backup
     130! If 'dates' file exists in output directory, make a backup
    132131  inquire(file=path(2)(1:length(2))//'dates', exist=ldates_file)
    133132  if (ldates_file.and.init) then
    … …  
    258257
    259258    write(anspec,'(i3.3)') ks
    260     if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
    261       if (ldirect.eq.1) then
    262         open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
     259
     260    if (DRYBKDEP.or.WETBKDEP) then !scavdep output
     261      if (DRYBKDEP) &
     262      open(unitoutgrid,file=path(2)(1:length(2))//'grid_drydep_'//adate// &
     263           atime//'_'//anspec,form='unformatted')
     264      if (WETBKDEP) &
     265      open(unitoutgrid,file=path(2)(1:length(2))//'grid_wetdep_'//adate// &
     266           atime//'_'//anspec,form='unformatted')
     267      write(unitoutgrid) itime
     268    else
     269      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     270        if (ldirect.eq.1) then
     271          open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
     272               atime//'_'//anspec,form='unformatted')
     273        else
     274          open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
     275               atime//'_'//anspec,form='unformatted')
     276        endif
     277        write(unitoutgrid) itime
     278      endif
     279
     280      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
     281        open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
    263282             atime//'_'//anspec,form='unformatted')
    264       else
    265         open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
    266              atime//'_'//anspec,form='unformatted')
     283
     284        write(unitoutgridppt) itime
    267285      endif
    268       write(unitoutgrid) itime
    269     endif
    270 
    271     if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
    272       open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
    273            atime//'_'//anspec,form='unformatted')
    274 
    275       write(unitoutgridppt) itime
    276     endif
     286    endif ! if deposition output
    277287
    278288    do kp=1,maxpointspec_act
    … …  
    354364! Concentration output
    355365!*********************
    356         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
     366        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
    357367
    358368! Wet deposition
    … …  
    449459                  endif
    450460                  sp_count_r=sp_count_r+1
     461                  if (lparticlecountoutput) then
     462                    sparse_dump_r(sp_count_r)= &
     463                         sp_fact* &
     464                         grid(ix,jy,kz)
     465                  else
    451466                  sparse_dump_r(sp_count_r)= &
    452467                       sp_fact* &
    453468                       grid(ix,jy,kz)* &
    454469                       factor3d(ix,jy,kz)/tot_mu(ks,kp)
     470                  end if
     471
     472
    455473!                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
    456474!    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
    … …  
    638656!*************************
    639657
    640   do ks=1,nspec
    641     do kp=1,maxpointspec_act
    642       do i=1,numreceptor
    643         creceptor(i,ks)=0.
    644       end do
    645       do jy=0,numygrid-1
    646         do ix=0,numxgrid-1
    647           do l=1,nclassunc
    648             do nage=1,nageclass
    649               do kz=1,numzgrid
    650                 gridunc(ix,jy,kz,ks,kp,l,nage)=0.
    651               end do
    652             end do
    653           end do
    654         end do
    655       end do
    656     end do
    657   end do
     658  creceptor(:,:)=0.
     659  gridunc(:,:,:,:,:,:,:)=0.
    658660
    659661  if (mp_measure_time) call mpif_mtime('rootonly',1)

  • src/getfields_mpi.f90

    r6ecb30a r20963b1  
    23323340  indmin=indj
    234234
     235   if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
     236        call writeprecip(itime,memind(1))
     237      endif
     238
    235239  else
    236240
    … …  
    28729160  indmin=indj
    288292
    289     mind3=memstat
     293!   if (WETBKDEP.and.lroot) then
     294    if (WETBKDEP.and.(lmpreader.or.(.not.lmp_use_reader.and.lroot))) then
     295      call writeprecip(itime,memind(1))
     296    endif
    290297
    291298  endif

  • src/makefile

    rccc9ec2 ra9cf4b1  
    3131#      make [-j] mpi-dbg
    3232#
     33#  NETCDF OUTPUT
     34#    To add support for output in netCDF format, append `ncf=yes` to the
     35#    `make` command
     36#
    3337################################################################################
    3438
    … …  
    5761        INCPATH2  = ${ROOT_DIR}/include
    5862        LIBPATH1 = ${ROOT_DIR}/lib
    59 
    60 else #ifeq ($(gcc), 5.4)
     63else
    6164# Compiled libraries under user ~flexpart, gfortran v5.4
    6265        ROOT_DIR = /homevip/flexpart/
    … …  
    6871        INCPATH2  = /usr/include
    6972        LIBPATH1 = ${ROOT_DIR}/gcc-5.4.0/lib
    70 
    71 #else
    72 # Default: System libraries at NILU, gfortran v4.6
    73 #       F90       = /usr/bin/gfortran
    74 #       MPIF90    = /usr/bin/mpif90.openmpi
    75 
    76 #       INCPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/include
    77 #       INCPATH2 = /usr/include
    78 #       LIBPATH1 = /xnilu_wrk/projects/FLEXPART/flex_wrk/bin64/grib_api/lib
    7973endif
     74
     75
     76### Enable netCDF output?
     77ifeq ($(ncf), yes)
     78        NCOPT = -DUSE_NCF -lnetcdff     
     79else
     80        NCOPT = -UUSE_NCF
     81endif
     82
    8083
    8184
    … …  
    8992
    9093## LIBRARIES
    91 LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
    92 
    93 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
    94 
    95 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
     94#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
     95LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
     96
     97FFLAGS   = -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
     98
     99DBGFLAGS = -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
    96100
    97101LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
    … …  
    139143        getfields_mpi.o \
    140144        readwind_ecmwf_mpi.o
     145
     146OBJECTS_NCF = netcdf_output_mod.o
    141147
    142148OBJECTS = \
    … …  
    196202dynamic_viscosity.o     get_settling.o  \
    197203initialize_cbl_vel.o    re_initialize_particle.o \
    198 cbl.o                   netcdf_output_mod.o
     204cbl.o
     205
     206ifeq ($(ncf), yes)
     207        OBJECTS := $(OBJECTS) $(OBJECTS_NCF)
     208endif
    199209
    200210%.o: %.mod

  • src/readOHfield.f90

    rd7aab4b r3f149cc  
    2222subroutine readOHfield
    2323
    24   !*****************************************************************************
    25   !                                                                            *
    26   ! Reads the OH field into memory                                             *
    27   !                                                                            *
    28   ! AUTHOR: R.L. Thompson, Nov 2014                                            *
    29   !                                                                            *
    30   !*****************************************************************************
    31   !                                                                            *
    32   ! Variables:                                                                 *
    33   !                                                                            *
    34   ! path(numpath)              contains the path names                         *
    35   ! lonOH(nxOH)                longitude of OH fields                          *
    36   ! latOH(nyOH)                latitude of OH fields                           *
    37   ! altOH(nzOH)                altitude of OH fields                           *
    38   ! etaOH(nzOH)                eta-levels of OH fields                         *
    39   ! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3)                *
    40   !                                                                            *
    41   !                                                                            *
    42   !*****************************************************************************
     24!*****************************************************************************
     25!                                                                            *
     26! Reads the OH field into memory                                             *
     27!                                                                            *
     28! AUTHOR: R.L. Thompson, Nov 2014                                            *
     29!                                                                            *
     30! UPDATES:                                                                   *
     31!   03/2018 SEC: Converted original netCDF files to binary format            *
     32!*****************************************************************************
     33!                                                                            *
     34! Variables:                                                                 *
     35!                                                                            *
     36! path(numpath)              contains the path names                         *
     37! lonOH(nxOH)                longitude of OH fields                          *
     38! latOH(nyOH)                latitude of OH fields                           *
     39! altOH(nzOH)                altitude of OH fields                           *
     40! etaOH(nzOH)                eta-levels of OH fields                         *
     41! OH_field(nxOH,nyOH,nzOH,m) OH concentration (molecules/cm3)                *
     42!                                                                            *
     43!                                                                            *
     44!*****************************************************************************
    4345
    4446  use oh_mod
    … …  
    4850  implicit none
    4951
    50   include 'netcdf.inc'
    51 
    52   character(len=150) :: thefile
    53   character(len=2) :: mm
    54   integer :: nid,ierr,xid,yid,zid,vid,m
     52  integer :: i,j,k,l,ierr
    5553  real, dimension(:), allocatable :: etaOH
    5654
    … …  
    6159  real, parameter :: scalehgt=7000. ! scale height in metres
    6260
    63   ! Read OH fields and level heights
    64   !********************************
    6561
    66   do m=1,12
    67  
    68     ! open netcdf file
    69     write(mm,fmt='(i2.2)') m
    70 !    thefile=trim(path(1))//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
    71     thefile=trim(ohfields_path)//'OH_FIELDS/'//'geos-chem.OH.2005'//mm//'01.nc'
    72     ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
    73     if(ierr.ne.0) then
    74       write (*,*) 'The OH field at: '//thefile// ' could not be opened'
    75       write (*,*) 'please copy the OH fields there, or change the path in the'
    76       write (*,*) 'COMMAND namelist!'
    77       write(*,*) nf_strerror(ierr)
    78       stop
    79     endif
     62  open(unitOH,file=trim(ohfields_path) &
     63       //'OH_FIELDS/OH_variables.bin',status='old', &
     64       form='UNFORMATTED', iostat=ierr, convert='little_endian')
    8065
    81     ! inquire about variables
    82     ierr=nf_inq_dimid(nid,'Lon-000',xid)
    83     if(ierr.ne.0) then
    84       write(*,*) nf_strerror(ierr)
    85       stop
    86     endif
    87     ierr=nf_inq_dimid(nid,'Lat-000',yid)
    88     if(ierr.ne.0) then
    89       write(*,*) nf_strerror(ierr)
    90       stop
    91     endif
    92     ierr=nf_inq_dimid(nid,'Alt-000',zid)
    93     if(ierr.ne.0) then
    94       write(*,*) nf_strerror(ierr)
    95       stop
    96     endif
    97 
    98     if(m.eq.1) then
    99 
    100       ! read dimension sizes
    101       ierr=nf_inq_dimlen(nid,xid,nxOH)
    102       if(ierr.ne.0) then
    103         write(*,*) nf_strerror(ierr)
    104         stop
    105       endif
    106       ierr=nf_inq_dimlen(nid,yid,nyOH)
    107       if(ierr.ne.0) then
    108         write(*,*) nf_strerror(ierr)
    109         stop
    110       endif
    111       ierr=nf_inq_dimlen(nid,zid,nzOH)
    112       if(ierr.ne.0) then
    113         write(*,*) nf_strerror(ierr)
    114         stop
    115       endif 
    116 
    117       ! allocate variables
    118       allocate(lonOH(nxOH))
    119       allocate(latOH(nyOH))
    120       allocate(etaOH(nzOH))
    121       allocate(altOH(nzOH))
    122       allocate(OH_field(nxOH,nyOH,nzOH,12))
    123       allocate(OH_hourly(nxOH,nyOH,nzOH,2))
    124 
    125       ! read dimension variables
    126       ierr=nf_inq_varid(nid,'LON',xid)
    127       ierr=nf_get_var_real(nid,xid,lonOH)
    128       if(ierr.ne.0) then
    129         write(*,*) nf_strerror(ierr)
    130         stop
    131       endif
    132       ierr=nf_inq_varid(nid,'LAT',yid)
    133       ierr=nf_get_var_real(nid,yid,latOH)
    134       if(ierr.ne.0) then
    135         write(*,*) nf_strerror(ierr)
    136         stop
    137       endif
    138       ierr=nf_inq_varid(nid,'ETAC',zid)
    139       ierr=nf_get_var_real(nid,zid,etaOH)
    140       if(ierr.ne.0) then
    141         write(*,*) nf_strerror(ierr)
    142         stop
    143       endif
    144 
    145       ! convert eta-level to altitude (assume surface pressure of 1010 hPa)
    146       altOH=log(1010./(etaOH*1010.))*scalehgt
    147 
    148     endif ! m.eq.1
    149 
    150     ! read OH_field
    151     ierr=nf_inq_varid(nid,'CHEM-L_S__OH',vid)
    152     ierr=nf_get_var_real(nid,vid,OH_field(:,:,:,m))
    153     if(ierr.ne.0) then
    154       write(*,*) nf_strerror(ierr)
    155       stop
    156     endif
    157 
    158     ierr=nf_close(nid)
    159 
    160   end do
    161  
    162   deallocate(etaOH)
    163 
    164   ! Read J(O1D) photolysis rates
    165   !******************************** 
    166 
    167   ! open netcdf file
    168 !  thefile=trim(path(1))//'OH_FIELDS/jrate_average.nc'
    169   thefile=trim(ohfields_path)//'OH_FIELDS/jrate_average.nc'
    170   ierr=nf_open(trim(thefile),NF_NOWRITE,nid)
    17166  if(ierr.ne.0) then
    172     write(*,*) nf_strerror(ierr)
     67    write(*,*) 'Cannot read binary OH fields in ',trim(ohfields_path)//'OH_FIELDS/OH_variables.bin'
    17368    stop
    17469  endif
    17570
    176   ! read dimension variables
    177   ierr=nf_inq_varid(nid,'longitude',xid)
    178   ierr=nf_get_var_real(nid,xid,lonjr)
    179   if(ierr.ne.0) then
    180     write(*,*) nf_strerror(ierr)
    181     stop
    182   endif
    183   ierr=nf_inq_varid(nid,'latitude',yid)
    184   ierr=nf_get_var_real(nid,yid,latjr)
    185   if(ierr.ne.0) then
    186     write(*,*) nf_strerror(ierr)
    187     stop
    188   endif
     71  read(unitOH) nxOH
     72  read(unitOH) nyOH
     73  read(unitOH) nzOH
     74  write(*,*) nxOH,nyOH,nzOH
    18975
    190   ! read jrate_average
    191   ierr=nf_inq_varid(nid,'jrate',vid)
    192   ierr=nf_get_var_real(nid,vid,jrate_average)
    193   if(ierr.ne.0) then
    194     write(*,*) nf_strerror(ierr)
    195     stop
    196   endif
     76! allocate variables
     77  allocate(lonOH(nxOH))
     78  allocate(latOH(nyOH))
     79  allocate(etaOH(nzOH))
     80  allocate(altOH(nzOH))
     81  allocate(OH_field(nxOH,nyOH,nzOH,12))
     82  allocate(OH_hourly(nxOH,nyOH,nzOH,2))
    19783
    198   ierr=nf_close(nid)
     84  read(unitOH) (lonjr(i),i=1,360)
     85  read(unitOH) (latjr(i),i=1,180)
     86  read(unitOH) (((jrate_average(i,j,k),i=1,360),j=1,180),k=1,12)
     87  read(unitOH) (lonOH(i),i=1,nxOH)
     88  read(unitOH) (latOH(i),i=1,nyOH)
     89  read(unitOH) (lonOH(i),i=1,nxOH)
    19990
    200   return
     91  read(unitOH) (altOH(i),i=1,nzOH)
     92  read(unitOH) ((((OH_field(i,j,k,l),i=1,nxOH),j=1,nyOH),k=1,nzOH),l=1,12)
     93  read(unitOH) ((((OH_hourly(i,j,k,l),i=1,nxOH),j=1,nyOH),k=1,nzOH),l=1,2)
    20194
    20295end subroutine readOHfield

  • src/readcommand.f90

    r01e6052 r20963b1  
    332332     case (3)  ! 3 .. wet deposition in outputfield
    333333        ind_rel = 3
    334          write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
    335          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
    336          write(*,*) ' #### Release is performed above ground lev    #### '
     334        if (lroot) then
     335          write(*,*) ' #### FLEXPART WET DEPOSITION BACKWARD MODE    #### '
     336          write(*,*) ' #### Releaseheight is forced to 0 - 20km      #### '
     337          write(*,*) ' #### Release is performed above ground lev    #### '
     338        end if
    337339         WETBKDEP=.true.
    338340         allocate(xscav_frac1(maxpart,maxspec))
    339341     case (4)  ! 4 .. dry deposition in outputfield
    340342         ind_rel = 4
    341          write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
    342          write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
    343          write(*,*) ' #### Release is performed above ground lev    #### '
     343         if (lroot) then
     344           write(*,*) ' #### FLEXPART DRY DEPOSITION BACKWARD MODE    #### '
     345           write(*,*) ' #### Releaseheight is forced to 0 - 2*href    #### '
     346           write(*,*) ' #### Release is performed above ground lev    #### '
     347         end if
    344348         DRYBKDEP=.true.
    345349         allocate(xscav_frac1(maxpart,maxspec))
    … …  
    392396  endif
    393397
    394 !  check for netcdf output switch (use for non-namelist input only!)
     398! Check for netcdf output switch
     399!*******************************
    395400  if (iout.ge.8) then
    396401     lnetcdfout = 1
    397402     iout = iout - 8
    398 ! #ifndef NETCDF_OUTPUT
    399 !      print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
    400 !      print*,'Please recompile with netcdf library or use standard output format.'
    401 !      stop
    402 ! #endif
     403#ifndef USE_NCF
     404     write(*,*) 'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
     405     write(*,*) 'Please recompile with netcdf library (`make [...] ncf=yes`) or use standard output format.'
     406     stop
     407#endif
    403408  endif
    404409

  • src/readwind_ecmwf_mpi.f90

    r61e07ba r20963b1  
    103103
    104104  integer :: isec1(56),isec2(22+nxmax+nymax)
    105   real(kind=4) :: zsec4(jpunp)
    106   real(kind=4) :: xaux,yaux
    107   real(kind=8) :: xauxin,yauxin
    108   real,parameter :: eps=1.e-4
    109   real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
    110   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
     105  real(sp) :: zsec4(jpunp)
     106  real(sp) :: xaux,yaux
     107  real(dp) :: xauxin,yauxin
     108  real(sp),parameter :: eps=1.e-4
     109  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
     110  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    111111
    112112  logical :: hflswitch,strswitch !,readcloud
    … …  
    123123!ZHG test the grib fields that have lcwc without using them
    124124!  readcloud=.false.
    125 !hg end
     125
    126126  levdiff2=nlev_ec-nwz+1
    127127  iumax=0

  • src/releaseparticles_mpi.f90

    r16b61a5 r20963b1  
    216216              xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
    217217                   *timecorrect(k)/average_timecorrect
    218 !             write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
     218              if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
     219!              if ( henry(k).gt.0 .or. &
     220!                   crain_aero(k).gt.0. .or. csnow_aero(k).gt.0. .or. &
     221!                   ccn_aero(k).gt.0. .or. in_aero(k).gt.0. )  then
     222                xscav_frac1(ipart,k)=-1.
     223               endif
    219224! Assign certain properties to particle
    220225!**************************************
    … …  
    366371!Af ind_rel is defined in readcommand.f
    367372
    368             if (ind_rel .eq. 1) then
     373            if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
    369374
    370375! Interpolate the air density

  • src/timemanager.f90

    rfe32dca ra9cf4b1  
    101101  use par_mod
    102102  use com_mod
     103#ifdef USE_NCF
    103104  use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
    104105       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
     106#endif
    105107
    106108  implicit none
    … …  
    389391          if (surf_only.ne.1) then
    390392            if (lnetcdfout.eq.1) then
     393#ifdef USE_NCF
    391394              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     395#endif
    392396            else
    393397              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
    … …  
    399403             write(*,*) 'system clock',count_clock - count_clock0   
    400404            endif
    401             if (lnetcdfout.eq.1) then
     405            if (lnetcdfout.eq.1) then
     406#ifdef USE_NCF
    402407              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     408#endif
    403409            else
    404410              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
    … …  
    419425              endif
    420426            else
     427#ifdef USE_NCF
    421428              if (surf_only.ne.1) then
    422429                call concoutput_nest_netcdf(itime,outnum)
    … …  
    424431                call concoutput_surf_nest_netcdf(itime,outnum)
    425432              endif
     433#endif
    426434            endif
    427435          endif

  • src/timemanager_mpi.f90

    rd8eed02 r20963b1  
    102102  use com_mod
    103103  use mpi_mod
     104#ifdef USE_NCF
    104105  use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
    105106       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
     107#endif
    106108
    107109  implicit none
    … …  
    113115  integer :: ip
    114116  integer :: loutnext,loutstart,loutend
    115   integer :: ix,jy,ldeltat,itage,nage
     117  integer :: ix,jy,ldeltat,itage,nage,idummy
    116118  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
    117119  integer :: numpart_tot_mpi ! for summing particles on all processes
    118   real :: outnum,weight,prob(maxspec)
    119   real :: decfact
     120  real :: outnum,weight,prob(maxspec), prob_rec(maxspec), decfact,wetscav
    120121
    121122  real(sp) :: gridtotalunc
    … …  
    123124       & drydeposit(maxspec)=0_dep_prec
    124125  real :: xold,yold,zold,xmassfract
     126  real :: grfraction(3)
    125127  real, parameter :: e_inv = 1.0/exp(1.0)
    126128
    … …  
    157159!CGZ-lifetime: set lifetime to 0
    158160
     161  if (.not.lusekerneloutput) write(*,*) 'Not using the kernel'
     162  if (turboff) write(*,*) 'Turbulence switched off'
    159163
    160164
    … …  
    481485            if (lroot) then
    482486              if (lnetcdfout.eq.1) then
     487#ifdef USE_NCF
    483488                call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
    484489                     &drygridtotalunc)
     490#endif
    485491              else
    486492                call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
    … …  
    494500            if (lroot) then
    495501              if (lnetcdfout.eq.1) then
     502#ifdef USE_NCF
    496503                call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
    497504                     &drygridtotalunc)
     505#endif
    498506              else
    499507                call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
    … …  
    513521            call mpif_tm_reduce_grid_nest
    514522 
    515            if (mp_measure_time) call mpif_mtime('iotime',0)
     523            if (mp_measure_time) call mpif_mtime('iotime',0)
    516524
    517525            if (lnetcdfout.eq.0) then
    … …  
    526534              else  ! :TODO: check for zeroing in the netcdf module
    527535                call concoutput_surf_nest(itime,outnum)
    528 
    529536              end if
    530 
    531537            else
    532 
     538#ifdef USE_NCF
    533539              if (surf_only.ne.1) then
    534540                if (lroot) then             
    … …  
    544550                end if
    545551              endif
    546 
    547 
     552#endif
    548553            end if
    549554          end if
    550          
    551 
    552555          outnum=0.
    553556        endif
    … …  
    707710        zold=ztra1(j)
    708711
     712   
     713  ! RECEPTOR: dry/wet depovel
     714  !****************************
     715  ! Before the particle is moved
     716  ! the calculation of the scavenged mass shall only be done once after release
     717  ! xscav_frac1 was initialised with a negative value
     718
     719      if  (DRYBKDEP) then
     720       do ks=1,nspec
     721         if  ((xscav_frac1(j,ks).lt.0)) then
     722            call get_vdep_prob(itime,xtra1(j),ytra1(j),ztra1(j),prob_rec)
     723            if (DRYDEPSPEC(ks)) then        ! dry deposition
     724               xscav_frac1(j,ks)=prob_rec(ks)
     725             else
     726                xmass1(j,ks)=0.
     727                xscav_frac1(j,ks)=0.
     728             endif
     729         endif
     730        enddo
     731       endif
     732
     733       if (WETBKDEP) then
     734       do ks=1,nspec
     735         if  ((xscav_frac1(j,ks).lt.0)) then
     736            call get_wetscav(itime,lsynctime,loutnext,j,ks,grfraction,idummy,idummy,wetscav)
     737            if (wetscav.gt.0) then
     738                xscav_frac1(j,ks)=wetscav* &
     739                       (zpoint2(npoint(j))-zpoint1(npoint(j)))*grfraction(1)
     740            else
     741                xmass1(j,ks)=0.
     742                xscav_frac1(j,ks)=0.
     743            endif
     744         endif
     745        enddo
     746       endif
     747
    709748! Integrate Lagevin equation for lsynctime seconds
    710749!*************************************************
    … …  
    808847                 call initial_cond_calc(itime+lsynctime,j)
    809848            itra1(j)=-999999999
    810             !print*, 'terminated particle ',j,'for age'
     849            if (verbosity.gt.0) then
     850              print*, 'terminated particle ',j,'for age'
     851            endif
    811852          endif
    812853        endif
    … …  
    819860
    820861
    821 ! Added by mc: counter of "unstable" particle velocity during a time scale
    822 !   of maximumtl=20 minutes (defined in com_mod)
    823 
     862! Counter of "unstable" particle velocity during a time scale
     863! of maximumtl=20 minutes (defined in com_mod)
     864!************************************************************
    824865    total_nan_intl=0
    825866    i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG