Changeset f203036 in flexpart.git


Ignore:
Timestamp:
Jun 18, 2019, 1:59:41 PM (5 years ago)
Author:
Ignacio Pisso <Ignacio.Pisso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug
Children:
79d7d08
Parents:
6741557 (diff), 941db73 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'release-10' into dev

Files:
6 added
29 edited

Legend:

Unmodified
Added
Removed
  • create_tarball.sh

    ra2e9de4 r941db73  
    11#!/bin/bash
    22
     3echo CREATE A NEW FLEXPART DISTRIBUTION
     4
    35#define version number
    4 
    56githash=$(git rev-parse --short --verify HEAD)
    6 
    7 
    8 version=10.3beta5_$githash
     7echo githash $githash
     8
     9version=10.3.1_$githash
     10echo version $version 
    911
    1012# define tarball name
    1113targetdir=../flexpart_distribution/
    12 tarball_tmp=${targetdir}flexpart_v$version
     14echo targetdir $targetdir
     15
     16distribution_name=flexpart_v$version
     17
     18tarball_tmp=${targetdir}flexpart_v$version
     19echo tarball_tmp $tarball_tmp
     20
    1321#tarball=${targetdir}flexpart_v$version.tar
    1422tarball=${tarball_tmp}.tar
     23echo tarball $tarball
    1524
    1625# clean old package
    1726if [ -d $tarball_tmp ]; then
    18   echo $tarball_tmp exists: move to $tarball_tmp.bk and exit 
    19   mkdir $tarball_tmp.bk
    20   mv $tarball_tmp ${tarball_tmp}.bk/
    21   mv $tarball ${tarball_tmp}.bk/
    22   exit
     27  echo
     28  echo clean old tarball
     29  hora=$(date +"%Y-%m-%d_%H%M%S")
     30  tarball_tmp_bk=$tarball_tmp$tarball_tmp_$hora
     31  echo tarball_tmp=$tarball_tmp exists: move to tarball_tmp_bk=$tarball_tmp_bk #and exit 
     32  mkdir $tarball_tmp_bk
     33  mv $tarball_tmp $tarball_tmp_bk/
     34  mv $tarball $tarball_tmp_bk/
     35  #exit
     36  echo old files moved to tarball_tmp_bk=$tarball_tmp_bk
     37  echo
    2338fi
    2439
    2540echo ---------------------------------------------------------
    26 echo ')' create basic dir structure
     41echo ')' create basis dir $tarball_tmp
    2742mkdir $tarball_tmp
    2843echo ---------------------------------------------------------
    29 echo ---------------------------------------------------------
    30 ##############################################################
    31 echo ')' pathnames
     44
     45echo
     46
     47echo ---------------------------------------------------------
     48echo ')' copy pathnames
    3249#cp pathnames_distribution $tarball_tmp/pathnames
    3350cp pathnames $tarball_tmp/pathnames
    3451echo ---------------------------------------------------------
    35 ##############################################################
    36 echo ')' src/
     52
     53echo
     54
     55echo ---------------------------------------------------------
     56echo ')' copy src/
    3757mkdir $tarball_tmp/src
    3858cp src/*.f90 $tarball_tmp/src
     
    4363echo ---------------------------------------------------------
    4464################################################################
    45 echo ')' options
     65
     66echo
     67
     68echo ---------------------------------------------------------
     69echo ')' copy options/
    4670# (for the distribution they work with the defult flex_ecmwf test winds)
    4771#cp -r options_flex_ecmwf_EA $tarball_tmp/options
     
    5478  echo $i
    5579  cp -r options/$i $tarball_tmp/options
     80  #echo copy $i to $tarball_tmp/options
    5681done
    5782
     
    6085cp options/SPECIES/SPECIES* $tarball_tmp/options/SPECIES/
    6186cp options/SPECIES/specoverview.f90 $tarball_tmp/options/SPECIES/
    62 echo ---------------------------------------------------------
    63 ################################################################
    64 echo ')' AVAILABLE
     87echo copy options/SPECIES/ to $tarball_tmp/options/SPECIES/
     88
     89echo ---------------------------------------------------------
     90
     91echo
     92
     93echo ---------------------------------------------------------
     94echo ')' copy AVAILABLE
    6595#cp AVAILABLE_flex_ecmwf_EA $tarball_tmp/AVAILABLE
    6696cp AVAILABLE $tarball_tmp/AVAILABLE
    67 
     97echo ---------------------------------------------------------
     98
     99echo
     100
     101echo ---------------------------------------------------------
     102echo  ')' create output/ #  mkdir $tarball_tmp/output
     103mkdir $tarball_tmp/output
    68104echo ---------------------------------------------------------
    69105################################################################
    70 echo  ')' output / #  mkdir $tarball_tmp/output
    71 mkdir $tarball_tmp/output
    72 echo ---------------------------------------------------------
    73 ################################################################
     106
     107echo
     108
     109echo ---------------------------------------------------------
    74110echo ')' preprocess/
    75111mkdir $tarball_tmp/preprocess
    76 #############################
     112
     113echo
     114
    77115echo -----------------flex_extract-------------------
    78116#echo '6)'  mkdir $tarball_tmp/flex_extract [a separate repository]
     
    88126## cp -r flex_extract/work/EA* $tarball_tmp/preprocess/flex_extract/work   
    89127
     128flex_extract=../flex_extract_v7.0.4/
    90129echo include flex_extract v7.0.4 b7c1c04a204c91e53759ef590504bf52dfaece64
    91 flex_extract=../flex_extract_v7.0.4/
     130echo from $flex_extract [use git modules?] IP 3/2018
     131
    92132cp $flex_extract/README.md $tarball_tmp/preprocess/flex_extract
    93133cp -r $flex_extract/docs $tarball_tmp/preprocess/flex_extract
     
    95135cp -r $flex_extract/python $tarball_tmp/preprocess/flex_extract
    96136cp -r $flex_extract/src $tarball_tmp/preprocess/flex_extract
    97 
    98 
    99 
     137echo flex_extract copied
     138echo ---------------------------------------------------------
    100139
    101140
     
    106145#cp -r examples/*.sh $tarball_tmp/examples/
    107146#cp -r examples/Makefile $tarball_tmp/examples/
     147
     148echo
     149
    108150echo ---------------------------------------------------------
    109151################################################################
    110 echo postprocess/
     152echo ')' postprocess/
    111153
    112154postprocess=postprocess
     
    119161cp $postprocess/flex_read_fortran/*.f90 $tarball_tmp/$postprocess/flex_read_fortran
    120162cp $postprocess/flex_read_fortran/makefile $tarball_tmp/$postprocess/flex_read_fortran
     163echo flex_read_fortran copied
    121164
    122165echo -----------------flex_read_matlab-------------------
     166echo flex_read_fortran NOT copied
    123167
    124168# add matlab reading routines
    125169#mkdir $tarball_tmp/postprocess/flex_read_matlab
    126170#cp postprocess/flex_read_matlab/*.m $tarball_tmp/postprocess/flex_read_matlab
    127 
    128 ###############################################################
    129 
    130 echo ---------------------------------------------------------
    131 echo tests/
    132 
     171echo ---------------------------------------------------------
     172
     173echo
     174
     175echo ---------------------------------------------------------
     176echo ')' tests/
     177###############################################################
    133178#echo '13) tests'
    134179mkdir $tarball_tmp/tests
    135 
    136 ###############################################################
    137180echo -----------------flex_read_fortran-------------------
    138 
    139181#echo 'b) ./tests/flex_read_fortran/'
    140182echo fixme
    141183#mkdir $tarball_tmp/tests/flex_read_fortran
    142184#cp tests/flex_read_fortran/test_read_default.sh  $tarball_tmp/tests/flex_read_fortran
     185
    143186
    144187###############################################################
     
    176219cp tests/declare_examples $tarball_tmp/tests/
    177220
    178 
     221echo
    179222
    180223# ~/repos/flexpart/tests$./compare_grids.sh
     
    182225#echo mkdir $tarball_tmp/tests/examples2/
    183226#echo cp tests/examples2/setup.sh $tarball_tmp/tests/examples2/
    184 echo --repeat examples-------------------
    185 #echo FIXME
    186 
    187 ###############################################################
    188 echo -----------------ctbto-------------------
    189 mkdir $tarball_tmp/tests/ctbto
     227# echo --repeat examples-------------------
     228# echo FIXME
     229
     230###############################################################
     231#echo -----------------ctbto-------------------
     232# mkdir $tarball_tmp/tests/ctbto
    190233
    191234# cp -r tests/NILU/test_1 $tarball_tmp/tests/
    192235# cp -r tests/default_cases $tarball_tmp/tests/
    193236
    194 tar cvf  $tarball  $tarball_tmp 
    195 
    196 echo  $tarball complete
     237echo ---------------------------------------------------------
     238echo create tarball
     239#tar cvf $tarball  $tarball_tmp 
     240#tar cf $tarball  $tarball_tmp 
     241#cd
     242
     243cd $targetdir
     244tar cf $distribution_name.tar $distribution_name
     245
     246pwd
     247
     248
     249echo  tarball $tarball complete
    197250echo exported untarred files in $tarball_tmp
    198251exit
     
    200253###############################################################
    201254
    202 
    203 
     255# obtain $FLEXHOME (and set)
     256#1 cd $FLEXHOME/src
     257
     258#2 compile
     259#
     260#[laptop] source /Users/ignacio/repos/flexpart/src/make_in_laptop.sh
     261# [njord] make
     262# ->created executable (FLEXPART)
     263
     264#3 execute in src (absolute paths)
     265#
     266#[laptop] cp  /Users/ignacio/repos/flexpart/src/pathnames .
     267#[njord] FIXME
     268#
     269# mkdir output
     270# ./FLEXPART
     271# ->created output in output/
     272
     273#4 read output
     274# cd  $FLEXHOME/postprocess/flex_read_fortran/
     275# make
     276# -> printheader* printgrid* flex_read_compare2*
     277#/postprocess/flex_read_fortran$./printheader ../../src/output/
     278#/postprocess/flex_read_fortran$./printgrid ../../src/output/ conc
     279# -> output in stdout (max:   11122924.0     sum:   90330784.0)
     280
     281#5 execute in $FLEXHOME
     282# cd $FLEXHOME
     283# get winds
     284#[laptop] cp -r ~/repos/flex_winds/work/ ./preprocess/flex_extract/
     285#[njord] curl https://folk.nilu.no/~ignacio/FLEXPART/EA120101.tar --output EA120101.tar ; tar -xvf EA120101.tar ; mv flex_extract/work preprocess/flex_extract/ ; rmdir flex_extract
     286
     287# src/FLEXPART
     288# -> output in $FLEXHOME/output/
     289
     290#6 read output
     291# postprocess/flex_read_fortran/printheader output/
     292# postprocess/flex_read_fortran/printgrid output/ conc
     293# -> output in stdout ( max:   11578738.0     sum:   104058720.)
     294
     295#7 gnererate examples
     296# cd $FLEXHOME/tests/examples
     297
     298#make run
     299
     300#make examples
     301#make batch
     302#./run_batch_cl.sh
     303
     304#make (set_default_example.sh)
     305#tests/examples$../../src/FLEXPART
     306#output
     307
     308#8 read examples:
     309#cd $FLEXHOME/tests/read_examples 
     310# ./read_headers.sh
     311# ./read_grids.sh
     312
     313#9 compare examples with reference
  • options/COMMAND

    r2753a5c r0a94e13  
    1919 IFINE=                 4, ! Reduction for time step in vertical transport, used only if CTL>1
    2020 IOUT=                  1, ! Output type: [1]mass 2]pptv 3]1&2 4]plume 5]1&4, +8 for NetCDF output     
    21  IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end  
     21 IPOUT=                 0, ! Particle position output: 0]no 1]every output 2]only at end 3]time averaged
    2222 LSUBGRID=              0, ! Increase of ABL heights due to sub-grid scale orographic variations;[0]off 1]on
    2323 LCONVECTION=           1, ! Switch for convection parameterization;0]off [1]on   
  • src/FLEXPART.f90

    r2753a5c r2eefa58  
    6868  integer :: detectformat
    6969
    70 
    71 
    72   ! Initialize arrays in com_mod
    73   !*****************************
    74   call com_mod_allocate_part(maxpart)
    7570 
    76 
    7771  ! Generate a large number of random numbers
    7872  !******************************************
     
    143137    write(*,*) 'call readpaths'
    144138  endif
    145   call readpaths(pathfile)
     139  call readpaths
    146140 
    147141  if (verbosity.gt.1) then !show clock info
     
    171165    endif     
    172166  endif
     167
     168  ! Initialize arrays in com_mod
     169  !*****************************
     170  call com_mod_allocate_part(maxpart)
     171
    173172
    174173  ! Read the age classes to be used
     
    453452  endif
    454453
     454  if (verbosity.gt.0) write (*,*) 'timemanager> call wetdepo'
    455455  call timemanager(metdata_format)
     456 
    456457
    457458  if (verbosity.gt.0) then
     
    468469      endif
    469470    end do
    470     write (*,*) 'timemanager> call wetdepo'
    471471  endif
    472472 
  • src/FLEXPART_MPI.f90

    r20963b1 r2eefa58  
    7777  if (mp_measure_time) call mpif_mtime('flexpart',0)
    7878
    79   ! Initialize arrays in com_mod
    80   !*****************************
    81 
    82   if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
    83 
    84 
     79 
    8580  ! Generate a large number of random numbers
    8681  !******************************************
     
    151146    write(*,*) 'call readpaths'
    152147  endif
    153   call readpaths(pathfile)
     148  call readpaths
    154149 
    155150  if (verbosity.gt.1) then !show clock info
     
    179174    endif
    180175  endif
     176
     177  ! Initialize arrays in com_mod
     178  !*****************************
     179
     180  if(.not.(lmpreader.and.lmp_use_reader)) call com_mod_allocate_part(maxpart_mpi)
    181181
    182182
     
    413413  end if ! (mpif_pid == 0)
    414414
    415   if (mp_measure_time) call mpif_mtime('iotime',0)
     415  if (mp_measure_time) call mpif_mtime('iotime',1)
    416416
    417417  if (verbosity.gt.0 .and. lroot) then
  • src/com_mod.f90

    re9e0f06 r2eefa58  
    1818
    1919  implicit none
     20
     21
    2022
    2123  !****************************************************************
     
    6971
    7072  real :: ctl,fine
    71   integer :: ifine,iout,ipout,ipin,iflux,mdomainfill
     73  integer :: ifine,iout,ipout,ipin,iflux,mdomainfill,ipoutfac
    7274  integer :: mquasilag,nested_output,ind_source,ind_receptor
    7375  integer :: ind_rel,ind_samp,ioutputforeachrelease,linit_cond,surf_only
     
    8284  ! iout     output options: 1 conc. output (ng/m3), 2 mixing ratio (pptv), 3 both
    8385  ! ipout    particle dump options: 0 no, 1 every output interval, 2 only at end
     86  ! ipoutfac increase particle dump interval by factor (default 1)
    8487  ! ipin     read in particle positions from dumped file from a previous run
    8588  ! fine     real(ifine)
     
    121124  ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
    122125
     126  integer :: linversionout
     127  ! linversionout 1 for one grid_time output file for each release containing all timesteps
     128
    123129  integer :: nageclass,lage(maxageclass)
    124130
     
    128134
    129135  logical :: gdomainfill
    130 
    131136  ! gdomainfill             .T., if domain-filling is global, .F. if not
    132137
     
    174179  real :: ri(5,numclass),rac(5,numclass),rcl(maxspec,5,numclass)
    175180  real :: rgs(maxspec,5,numclass),rlu(maxspec,5,numclass)
    176   real :: rm(maxspec),dryvel(maxspec)
     181  real :: rm(maxspec),dryvel(maxspec),kao(maxspec)
    177182  real :: ohcconst(maxspec),ohdconst(maxspec),ohnconst(maxspec)
    178183
     
    359364  real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice      [kg/kg]
    360365  real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0  !combined [m3/m3]
    361 
     366! RLT add pressure and dry air density
     367  real :: prs(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
     368  real :: rho_dry(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
    362369  real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
    363370  real :: rho(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
     
    381388  ! uupol,vvpol [m/s]    wind components in polar stereographic projection
    382389  ! tt [K]               temperature data
     390  ! prs                  air pressure
    383391  ! qv                   specific humidity data
    384392  ! pv (pvu)             potential vorticity
     
    651659  real :: receptorarea(maxreceptor)
    652660  real :: creceptor(maxreceptor,maxspec)
     661  real, allocatable, dimension(:,:) :: creceptor0
    653662  character(len=16) :: receptorname(maxreceptor)
    654663  integer :: numreceptor
     
    673682  real, allocatable, dimension(:,:) :: xmass1
    674683  real, allocatable, dimension(:,:) :: xscav_frac1
     684
     685! Variables used for writing out interval averages for partoutput
     686!****************************************************************
     687
     688  integer, allocatable, dimension(:) :: npart_av
     689  real, allocatable, dimension(:) :: part_av_cartx,part_av_carty,part_av_cartz,part_av_z,part_av_topo
     690  real, allocatable, dimension(:) :: part_av_pv,part_av_qv,part_av_tt,part_av_rho,part_av_tro,part_av_hmix
     691  real, allocatable, dimension(:) :: part_av_uu,part_av_vv,part_av_energy
    675692
    676693  ! eso: Moved from timemanager
     
    780797         & idt(nmpart),itramem(nmpart),itrasplit(nmpart),&
    781798         & xtra1(nmpart),ytra1(nmpart),ztra1(nmpart),&
    782          & xmass1(nmpart, maxspec),&
    783          & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
     799         & xmass1(nmpart, maxspec))  ! ,&
     800!         & checklifetime(nmpart,maxspec), species_lifetime(maxspec,2))!CGZ-lifetime
     801
     802    if (ipout.eq.3) then
     803      allocate(npart_av(nmpart),part_av_cartx(nmpart),part_av_carty(nmpart),&
     804           & part_av_cartz(nmpart),part_av_z(nmpart),part_av_topo(nmpart))
     805      allocate(part_av_pv(nmpart),part_av_qv(nmpart),part_av_tt(nmpart),&
     806           & part_av_rho(nmpart),part_av_tro(nmpart),part_av_hmix(nmpart))
     807      allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart))
     808    end if
    784809
    785810
    786811    allocate(uap(nmpart),ucp(nmpart),uzp(nmpart),us(nmpart),&
    787812         & vs(nmpart),ws(nmpart),cbt(nmpart))
    788    
     813
    789814  end subroutine com_mod_allocate_part
    790815
  • src/concoutput.f90

    r20963b1 r2eefa58  
    7272  real :: sp_fact
    7373  real :: outnum,densityoutrecept(maxreceptor),xl,yl
     74! RLT
     75  real :: densitydryrecept(maxreceptor)
     76  real :: factor_dryrecept(maxreceptor)
    7477
    7578!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
     
    106109  character(LEN=8),save :: file_stat='REPLACE'
    107110  logical :: ldates_file
     111  logical :: lexist
    108112  integer :: ierr
    109113  character(LEN=100) :: dates_char
     
    203207        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
    204208             rho(iix,jjy,kzz-1,mind)*dz2)/dz
     209! RLT
     210        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
     211             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz 
    205212      end do
    206213    end do
     
    214221!densityoutrecept(i)=rho(iix,jjy,1,2)
    215222    densityoutrecept(i)=rho(iix,jjy,1,mind)
     223! RLT
     224    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
    216225  end do
    217226
     227! RLT
     228! conversion factor for output relative to dry air
     229  factor_drygrid=densityoutgrid/densitydrygrid
     230  factor_dryrecept=densityoutrecept/densitydryrecept
    218231
    219232! Output is different for forward and backward simulations
     
    353366! Concentration output
    354367!*********************
    355         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
     368        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
    356369
    357370! Wet deposition
     
    614627  end do
    615628
     629! RLT Aug 2017
     630! Write out conversion factor for dry air
     631  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
     632  if (lexist) then
     633    ! open and append
     634    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
     635            status='old',action='write',access='append')
     636  else
     637    ! create new
     638    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
     639            status='new',action='write')
     640  endif
     641  sp_count_i=0
     642  sp_count_r=0
     643  sp_fact=-1.
     644  sp_zer=.true.
     645  do kz=1,numzgrid
     646    do jy=0,numygrid-1
     647      do ix=0,numxgrid-1
     648        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
     649          if (sp_zer.eqv..true.) then ! first value not equal to one
     650            sp_count_i=sp_count_i+1
     651            sparse_dump_i(sp_count_i)= &
     652                  ix+jy*numxgrid+kz*numxgrid*numygrid
     653            sp_zer=.false.
     654            sp_fact=sp_fact*(-1.)
     655          endif
     656          sp_count_r=sp_count_r+1
     657          sparse_dump_r(sp_count_r)= &
     658               sp_fact*factor_drygrid(ix,jy,kz)
     659        else ! factor is one
     660          sp_zer=.true.
     661        endif
     662      end do
     663    end do
     664  end do
     665  write(unitoutfactor) sp_count_i
     666  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
     667  write(unitoutfactor) sp_count_r
     668  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
     669  close(unitoutfactor)
     670
     671
    616672  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
    617673  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
     
    640696  endif
    641697
    642 
     698! RLT Aug 2017
     699! Write out conversion factor for dry air
     700  if (numreceptor.gt.0) then
     701    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
     702     if (lexist) then
     703     ! open and append
     704      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
     705              status='old',action='write',access='append')
     706    else
     707      ! create new
     708      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
     709              status='new',action='write')
     710    endif
     711    write(unitoutfactor) itime
     712    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
     713    close(unitoutfactor)
     714  endif
    643715
    644716! Reinitialization of grid
  • src/concoutput_mpi.f90

    r20963b1 r6741557  
    364364! Concentration output
    365365!*********************
    366         if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5).or.(iout.eq.6)) then
     366        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
    367367
    368368! Wet deposition
  • src/concoutput_nest.f90

    r7d02c2f r2eefa58  
    7070  real :: sp_fact
    7171  real :: outnum,densityoutrecept(maxreceptor),xl,yl
     72! RLT
     73  real :: densitydryrecept(maxreceptor)
     74  real :: factor_dryrecept(maxreceptor)
    7275
    7376  !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
     
    9699  character :: adate*8,atime*6
    97100  character(len=3) :: anspec
     101  logical :: lexist
    98102  integer :: mind
    99103! mind        eso:added to ensure identical results between 2&3-fields versions
     
    164168        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
    165169             rho(iix,jjy,kzz-1,mind)*dz2)/dz
     170! RLT
     171        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,mind)*dz1+ &
     172             rho_dry(iix,jjy,kzz-1,mind)*dz2)/dz
    166173      end do
    167174    end do
     
    175182    !densityoutrecept(i)=rho(iix,jjy,1,2)
    176183    densityoutrecept(i)=rho(iix,jjy,1,mind)
     184! RLT
     185    densitydryrecept(i)=rho_dry(iix,jjy,1,mind)
    177186  end do
    178187
     188! RLT
     189! conversion factor for output relative to dry air
     190  factor_drygrid=densityoutgrid/densitydrygrid
     191  factor_dryrecept=densityoutrecept/densitydryrecept
    179192
    180193  ! Output is different for forward and backward simulations
     
    551564  end do
    552565
     566! RLT Aug 2017
     567! Write out conversion factor for dry air
     568  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
     569  if (lexist) then
     570    ! open and append
     571    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
     572            status='old',action='write',access='append')
     573  else
     574    ! create new
     575    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
     576            status='new',action='write')
     577  endif
     578  sp_count_i=0
     579  sp_count_r=0
     580  sp_fact=-1.
     581  sp_zer=.true.
     582  do kz=1,numzgrid
     583    do jy=0,numygridn-1
     584      do ix=0,numxgridn-1
     585        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
     586          if (sp_zer.eqv..true.) then ! first value not equal to one
     587            sp_count_i=sp_count_i+1
     588            sparse_dump_i(sp_count_i)= &
     589                  ix+jy*numxgridn+kz*numxgridn*numygridn
     590            sp_zer=.false.
     591            sp_fact=sp_fact*(-1.)
     592          endif
     593          sp_count_r=sp_count_r+1
     594          sparse_dump_r(sp_count_r)= &
     595               sp_fact*factor_drygrid(ix,jy,kz)
     596        else ! factor is one
     597          sp_zer=.true.
     598        endif
     599      end do
     600    end do
     601  end do
     602  write(unitoutfactor) sp_count_i
     603  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
     604  write(unitoutfactor) sp_count_r
     605  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
     606  close(unitoutfactor)
    553607
    554608
  • src/concoutput_surf.f90

    r16b61a5 r2eefa58  
    7272  real :: sp_fact
    7373  real :: outnum,densityoutrecept(maxreceptor),xl,yl
     74! RLT
     75  real :: densitydryrecept(maxreceptor)
     76  real :: factor_dryrecept(maxreceptor)
    7477
    7578!real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
     
    101104  character :: adate*8,atime*6
    102105  character(len=3) :: anspec
     106  logical :: lexist
    103107
    104108
     
    180184        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    181185             rho(iix,jjy,kzz-1,2)*dz2)/dz
     186! RLT
     187        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
     188             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
    182189      end do
    183190    end do
     
    190197    jjy=max(min(nint(yl),nymin1),0)
    191198    densityoutrecept(i)=rho(iix,jjy,1,2)
     199! RLT
     200    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
    192201  end do
    193202
     203! RLT
     204! conversion factor for output relative to dry air
     205  factor_drygrid=densityoutgrid/densitydrygrid
     206  factor_dryrecept=densityoutrecept/densitydryrecept
    194207
    195208! Output is different for forward and backward simulations
     
    596609  end do
    597610
     611! RLT Aug 2017
     612! Write out conversion factor for dry air
     613  inquire(file=path(2)(1:length(2))//'factor_drygrid',exist=lexist)
     614  if (lexist) then
     615    ! open and append
     616    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
     617            status='old',action='write',access='append')
     618  else
     619    ! create new
     620    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid',form='unformatted',&
     621            status='new',action='write')
     622  endif
     623  sp_count_i=0
     624  sp_count_r=0
     625  sp_fact=-1.
     626  sp_zer=.true.
     627  do kz=1,1
     628    do jy=0,numygrid-1
     629      do ix=0,numxgrid-1
     630        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
     631          if (sp_zer.eqv..true.) then ! first value not equal to one
     632            sp_count_i=sp_count_i+1
     633            sparse_dump_i(sp_count_i)= &
     634                  ix+jy*numxgrid+kz*numxgrid*numygrid
     635            sp_zer=.false.
     636            sp_fact=sp_fact*(-1.)
     637          endif
     638          sp_count_r=sp_count_r+1
     639          sparse_dump_r(sp_count_r)= &
     640               sp_fact*factor_drygrid(ix,jy,kz)
     641        else ! factor is one
     642          sp_zer=.true.
     643        endif
     644      end do
     645    end do
     646  end do
     647  write(unitoutfactor) sp_count_i
     648  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
     649  write(unitoutfactor) sp_count_r
     650  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
     651  close(unitoutfactor)
     652
     653
    598654  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
    599655  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
     
    622678  endif
    623679
    624 
     680! RLT Aug 2017
     681! Write out conversion factor for dry air
     682  if (numreceptor.gt.0) then
     683    inquire(file=path(2)(1:length(2))//'factor_dryreceptor',exist=lexist)
     684     if (lexist) then
     685     ! open and append
     686      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
     687              status='old',action='write',access='append')
     688    else
     689      ! create new
     690      open(unitoutfactor,file=path(2)(1:length(2))//'factor_dryreceptor',form='unformatted',&
     691              status='new',action='write')
     692    endif
     693    write(unitoutfactor) itime
     694    write(unitoutfactor) (factor_dryrecept(i),i=1,numreceptor)
     695    close(unitoutfactor)
     696  endif
    625697
    626698! Reinitialization of grid
  • src/concoutput_surf_nest.f90

    r6a678e3 r2eefa58  
    7070  real :: sp_fact
    7171  real :: outnum,densityoutrecept(maxreceptor),xl,yl
     72! RLT
     73  real :: densitydryrecept(maxreceptor)
     74  real :: factor_dryrecept(maxreceptor)
    7275
    7376  !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
     
    9699  character :: adate*8,atime*6
    97100  character(len=3) :: anspec
    98 
     101  logical :: lexist
    99102
    100103  ! Determine current calendar date, needed for the file name
     
    159162        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    160163             rho(iix,jjy,kzz-1,2)*dz2)/dz
     164! RLT
     165        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
     166             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
    161167      end do
    162168    end do
     
    169175      jjy=max(min(nint(yl),nymin1),0)
    170176      densityoutrecept(i)=rho(iix,jjy,1,2)
     177! RLT
     178    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
    171179    end do
    172180
     181! RLT
     182! conversion factor for output relative to dry air
     183  factor_drygrid=densityoutgrid/densitydrygrid
     184  factor_dryrecept=densityoutrecept/densitydryrecept
    173185
    174186  ! Output is different for forward and backward simulations
     
    317329         write(unitoutgrid) sp_count_r
    318330         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
    319          write(unitoutgrid) sp_count_r
    320          write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
     331!         write(unitoutgrid) sp_count_r
     332!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
    321333
    322334  ! Dry deposition
     
    354366         write(unitoutgrid) sp_count_r
    355367         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
    356          write(unitoutgrid) sp_count_r
    357          write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
     368!         write(unitoutgrid) sp_count_r
     369!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
    358370
    359371
     
    399411         write(unitoutgrid) sp_count_r
    400412         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
    401          write(unitoutgrid) sp_count_r
    402          write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
     413!         write(unitoutgrid) sp_count_r
     414!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
    403415         else
    404416
     
    440452         write(unitoutgrid) sp_count_r
    441453         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
    442          write(unitoutgrid) sp_count_r
    443          write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
     454!         write(unitoutgrid) sp_count_r
     455!         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
    444456         endif ! surf_only
    445457
     
    487499         write(unitoutgridppt) sp_count_r
    488500         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
    489          write(unitoutgridppt) sp_count_r
    490          write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
     501!         write(unitoutgridppt) sp_count_r
     502!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
    491503
    492504
     
    526538         write(unitoutgridppt) sp_count_r
    527539         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
    528          write(unitoutgridppt) sp_count_r
    529          write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
     540!         write(unitoutgridppt) sp_count_r
     541!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
    530542
    531543
     
    570582         write(unitoutgridppt) sp_count_r
    571583         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
    572          write(unitoutgridppt) sp_count_r
    573          write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
     584!         write(unitoutgridppt) sp_count_r
     585!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
    574586         else
    575587
     
    611623         write(unitoutgridppt) sp_count_r
    612624         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
    613          write(unitoutgridppt) sp_count_r
    614          write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
     625!         write(unitoutgridppt) sp_count_r
     626!         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
    615627         endif ! surf_only
    616628
     
    624636
    625637  end do
     638
     639! RLT Aug 2017
     640! Write out conversion factor for dry air
     641  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
     642  if (lexist) then
     643    ! open and append
     644    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
     645            status='old',action='write',access='append')
     646  else
     647    ! create new
     648    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
     649            status='new',action='write')
     650  endif
     651  sp_count_i=0
     652  sp_count_r=0
     653  sp_fact=-1.
     654  sp_zer=.true.
     655  do kz=1,1
     656    do jy=0,numygridn-1
     657      do ix=0,numxgridn-1
     658        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
     659          if (sp_zer.eqv..true.) then ! first value not equal to one
     660            sp_count_i=sp_count_i+1
     661            sparse_dump_i(sp_count_i)= &
     662                  ix+jy*numxgridn+kz*numxgridn*numygridn
     663            sp_zer=.false.
     664            sp_fact=sp_fact*(-1.)
     665          endif
     666          sp_count_r=sp_count_r+1
     667          sparse_dump_r(sp_count_r)= &
     668               sp_fact*factor_drygrid(ix,jy,kz)
     669        else ! factor is one
     670          sp_zer=.true.
     671        endif
     672      end do
     673    end do
     674  end do
     675  write(unitoutfactor) sp_count_i
     676  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
     677  write(unitoutfactor) sp_count_r
     678  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
     679  close(unitoutfactor)
    626680
    627681
  • src/getfields.f90

    r6ecb30a r2eefa58  
    8787  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
    8888  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
     89! RLT added partial pressure water vapor
     90  real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
     91  integer :: kz, ix
     92  character(len=100) :: rowfmt
    8993
    9094  integer :: indmin = 1
     
    15315740  indmin=indj
    154158
    155    if (WETBKDEP) then
    156         call writeprecip(itime,memind(1))
    157    endif
     159    if (WETBKDEP) then
     160      call writeprecip(itime,memind(1))
     161    endif
    158162
    159163  else
     
    20420860  indmin=indj
    205209
    206    if (WETBKDEP) then
    207         call writeprecip(itime,memind(1))
    208    endif
    209 
    210   endif
     210    if (WETBKDEP) then
     211      call writeprecip(itime,memind(1))
     212    endif
     213
     214  end if
     215
     216  ! RLT calculate dry air density
     217  pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
     218  rho_dry=(prs-pwater)/(r_air*tt)
     219
     220  ! test density
     221!  write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
     222!  if(itime.eq.0) then
     223!    open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
     224!    do kz=1,nzmax
     225!      do ix=1,nxmax
     226!        write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
     227!      end do
     228!    end do
     229!    close(500)
     230!    open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
     231!    do kz=1,nzmax
     232!      do ix=1,nxmax
     233!        write(500,fmt=rowfmt) rho(ix,:,kz,1)
     234!      end do
     235!    end do
     236!    close(500)
     237!  endif
    211238
    212239  lwindinterv=abs(memtime(2)-memtime(1))
  • src/init_domainfill.f90

    rb5127f9 r0a94e13  
    8686    endif
    8787  endif
     88
     89! Exit here if resuming a run from particle dump
     90!***********************************************
     91  if (gdomainfill.and.ipin.ne.0) return
    8892
    8993! Do not release particles twice (i.e., not at both in the leftmost and rightmost
     
    414418!***************************************************************************
    415419
    416   if (ipin.eq.1) then
     420  if ((ipin.eq.1).and.(.not.gdomainfill)) then
    417421    open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
    418422         form='unformatted')
  • src/init_domainfill_mpi.f90

    rb5127f9 r328fdf9  
    110110  endif
    111111
     112! Exit here if resuming a run from particle dump
     113!***********************************************
     114  if (gdomainfill.and.ipin.ne.0) return
     115
    112116! Do not release particles twice (i.e., not at both in the leftmost and rightmost
    113117! grid cell) for a global domain
     
    213217        colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
    214218        colmasstotal=colmasstotal+colmass(ix,jy)
    215 
    216219      end do
    217220    end do
     
    466469
    467470! eso TODO: only needed for root process
    468     if (ipin.eq.1) then
     471    if ((ipin.eq.1).and.(.not.gdomainfill)) then
    469472      open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
    470473           form='unformatted')
     
    474477    endif
    475478
    476     numpart = numpart/mp_partgroup_np
    477     if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
    478 
    479   else ! Allocate dummy arrays for receiving processes
    480     allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
    481          & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
    482          & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
    483          & xmass1_tmp(nullsize, nullsize))
     479    if (ipin.eq.0) then   
     480      numpart = numpart/mp_partgroup_np
     481      if (mod(numpart,mp_partgroup_np).ne.0) numpart=numpart+1
     482    end if
     483
     484  else ! Allocate dummy arrays for receiving processes
     485    if (ipin.eq.0) then   
     486      allocate(itra1_tmp(nullsize),npoint_tmp(nullsize),nclass_tmp(nullsize),&
     487           & idt_tmp(nullsize),itramem_tmp(nullsize),itrasplit_tmp(nullsize),&
     488           & xtra1_tmp(nullsize),ytra1_tmp(nullsize),ztra1_tmp(nullsize),&
     489           & xmass1_tmp(nullsize, nullsize))
     490    end if
    484491   
    485   end if ! end if(lroot) 
     492  end if ! end if(lroot)
     493
    486494
    487495
    488496! Distribute particles to other processes (numpart is 'per-process', not total)
    489   call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
    490 ! eso TODO: xmassperparticle: not necessary to send
    491   call MPI_Bcast(xmassperparticle, 1, mp_sp, id_root, mp_comm_used, mp_ierr)
    492   call mpif_send_part_properties(numpart)
     497! Only if not restarting from previous run
     498  if (ipin.eq.0) then
     499    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
     500    call mpif_send_part_properties(npart(1)/mp_partgroup_np)
    493501
    494502! Deallocate the temporary arrays used for all particles
    495   deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
     503    deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
    496504         & itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
     505  end if
    497506
    498507
  • src/makefile

    r7123c70 r2eefa58  
    118118OBJECTS_SERIAL = \
    119119        releaseparticles.o      partoutput.o \
     120        partoutput_average.o \
    120121        conccalc.o \
    121122        init_domainfill.o       concoutput.o  \
     
    127128        redist.o                \
    128129        concoutput_surf.o       concoutput_surf_nest.o  \
     130        concoutput_inversion_nest.o     \
     131        concoutput_inversion.o \
    129132        getfields.o \
    130133        readwind_ecmwf.o
     
    132135## For MPI version
    133136OBJECTS_MPI = releaseparticles_mpi.o partoutput_mpi.o \
    134         conccalc_mpi.o \
     137        partoutput_average_mpi.o conccalc_mpi.o \
    135138        init_domainfill_mpi.o concoutput_mpi.o  \
    136139        timemanager_mpi.o FLEXPART_MPI.o        \
     
    149152advance.o               initialize.o            \
    150153writeheader.o           writeheader_txt.o       \
    151 writeprecip.o \
     154partpos_average.o       writeprecip.o \
    152155writeheader_surf.o      assignland.o\
    153156part0.o                 gethourlyOH.o\
     
    199202drydepokernel_nest.o    zenithangle.o \
    200203ohreaction.o            getvdep_nests.o \
    201 initial_cond_calc.o     initial_cond_output.o \
     204initial_cond_calc.o     initial_cond_output.o initial_cond_output_inversion.o \
    202205dynamic_viscosity.o     get_settling.o  \
    203206initialize_cbl_vel.o    re_initialize_particle.o \
     
    271274conccalc_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o unc_mod.o
    272275concoutput.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
     276concoutput_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
    273277concoutput_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
    274278        unc_mod.o mean_mod.o
    275279concoutput_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
     280concoutput_inversion_nest.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o mean_mod.o
    276281concoutput_nest_mpi.o: com_mod.o mpi_mod.o outg_mod.o par_mod.o point_mod.o \
    277282        unc_mod.o mean_mod.o
     
    320325initial_cond_calc.o: com_mod.o outg_mod.o par_mod.o unc_mod.o
    321326initial_cond_output.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
     327initial_cond_output_inversion.o: com_mod.o outg_mod.o par_mod.o point_mod.o unc_mod.o
    322328initialize.o: com_mod.o hanna_mod.o interpol_mod.o par_mod.o random_mod.o
    323329initialize_cbl_vel.o: com_mod.o par_mod.o random_mod.o
     
    348354part0.o: par_mod.o
    349355partdep.o: par_mod.o
     356partpos_average.o: com_mod.o par_mod.o
    350357partoutput.o: com_mod.o par_mod.o
     358partoutput_average.o: com_mod.o par_mod.o
     359partoutput_average_mpi.o: com_mod.o par_mod.o mpi_mod.o
    351360partoutput_mpi.o: com_mod.o mpi_mod.o par_mod.o
    352361partoutput_short.o: com_mod.o par_mod.o
  • src/mpi_mod.f90

    r0ecc1fe r0c8c7f2  
    8888! Variables for MPI processes in the 'particle' group
    8989  integer, allocatable, dimension(:) :: mp_partgroup_rank
     90  integer, allocatable, dimension(:) :: npart_per_process
    9091  integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np
    9192
     
    125126! mp_time_barrier   Measure MPI barrier time
    126127! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
    127 !                   to standard output (this does *not* otherwise affect the simulation)
     128!                   to standard output (this does not otherwise affect the simulation)
    128129  logical, parameter :: mp_dbg_mode = .false.
    129130  logical, parameter :: mp_dev_mode = .false.
     
    190191!   mp_np       number of running processes, decided at run-time
    191192!***********************************************************************
    192     use par_mod, only: maxpart, numwfmem, dep_prec
    193     use com_mod, only: mpi_mode, verbosity
     193    use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec
     194    use com_mod, only: mpi_mode, verbosity, creceptor0
    194195
    195196    implicit none
     
    337338
    338339! Set maxpart per process
    339 ! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution
     340! ESO 08/2016: Increase maxpart per process, in case of unbalanced distribution
    340341    maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np))
    341342    if (mp_np == 1) maxpart_mpi = maxpart
     
    365366    end if
    366367
     368! Allocate array for number of particles per process   
     369    allocate(npart_per_process(0:mp_partgroup_np-1))
     370
     371! Write whether MPI_IN_PLACE is used or not
     372#ifdef USE_MPIINPLACE
     373    if (lroot) write(*,*) 'Using MPI_IN_PLACE operations'
     374#else
     375    if (lroot) allocate(creceptor0(maxreceptor,maxspec))
     376    if (lroot) write(*,*) 'Not using MPI_IN_PLACE operations'
     377#endif
    367378    goto 101
    368379
     
    559570! invalid particles at the end of the arrays
    560571
    561 601 do i=num_part, 1, -1
     572601 do i=numpart, 1, -1
    562573      if (itra1(i).eq.-999999999) then
    563574        numpart=numpart-1
     
    598609    integer :: i,jj,nn, num_part=1,m,imin, num_trans
    599610    logical :: first_iter
    600     integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
     611    integer,allocatable,dimension(:) :: idx_arr
    601612    real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
    602613
     
    607618! All processes exchange information on number of particles
    608619!****************************************************************************
    609     allocate(numparticles_mpi(0:mp_partgroup_np-1), &
    610          &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
    611 
    612     call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &
     620    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
     621
     622    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
    613623         & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
    614624
     
    616626! Sort from lowest to highest
    617627! Initial guess: correct order
    618     sorted(:) = numparticles_mpi(:)
     628    sorted(:) = npart_per_process(:)
    619629    do i=0, mp_partgroup_np-1
    620630      idx_arr(i) = i
    621631    end do
     632
     633! Do not rebalance particles for ipout=3   
     634    if (ipout.eq.3) return
    622635
    623636! For each successive element in index array, see if a lower value exists
     
    645658    m=mp_partgroup_np-1 ! index for last sorted process (most particles)
    646659    do i=0,mp_partgroup_np/2-1
    647       num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
     660      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
    648661      if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
    649         if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&
    650              & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then
     662        if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.&
     663             & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then
    651664! DBG
    652           ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &
    653           !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi
     665          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', &
     666          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process
    654667! DBG
    655668          call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
     
    659672    end do
    660673
    661     deallocate(numparticles_mpi, idx_arr, sorted)
     674    deallocate(idx_arr, sorted)
    662675
    663676  end subroutine mpif_calculate_part_redist
     
    19611974    if (readclouds) then
    19621975      j=j+1
    1963       call MPI_Irecv(ctwc(:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
     1976      call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
    19641977           &MPI_COMM_WORLD,reqs(j),mp_ierr)
    19651978      if (mp_ierr /= 0) goto 600
     
    23262339      if (readclouds) then
    23272340        j=j+1
    2328         call MPI_Irecv(ctwcn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
     2341        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
    23292342             &MPI_COMM_WORLD,reqs(j),mp_ierr)
    23302343        if (mp_ierr /= 0) goto 600
     
    24622475    end if
    24632476
     2477  ! Receptor concentrations   
     2478    if (lroot) then
     2479      call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
     2480           & mp_comm_used,mp_ierr)
     2481      if (mp_ierr /= 0) goto 600
     2482    else
     2483      call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
     2484           & mp_comm_used,mp_ierr)
     2485    end if
     2486
    24642487#else
    24652488
    24662489      call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
    24672490           & mp_comm_used, mp_ierr)
     2491      if (mp_ierr /= 0) goto 600
    24682492      if (lroot) gridunc = gridunc0
     2493
     2494      call MPI_Reduce(creceptor, creceptor0,rcpt_size,mp_sp,MPI_SUM,id_root, &
     2495           & mp_comm_used,mp_ierr)
     2496      if (mp_ierr /= 0) goto 600
     2497      if (lroot) creceptor = creceptor0
    24692498
    24702499#endif
     
    24822511    end if
    24832512
    2484 ! Receptor concentrations   
    2485     if (lroot) then
    2486       call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
    2487            & mp_comm_used,mp_ierr)
    2488       if (mp_ierr /= 0) goto 600
    2489     else
    2490       call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
    2491            & mp_comm_used,mp_ierr)
    2492     end if
    24932513
    24942514    if (mp_measure_time) call mpif_mtime('commtime',1)
     
    27002720      end if
    27012721
    2702     case ('readwind')
    2703       if (imode.eq.0) then
    2704         call cpu_time(mp_readwind_time_beg)
    2705         mp_readwind_wtime_beg = mpi_wtime()
    2706       else
    2707         call cpu_time(mp_readwind_time_end)
    2708         mp_readwind_wtime_end = mpi_wtime()
    2709 
    2710         mp_readwind_time_total = mp_readwind_time_total + &
    2711              &(mp_readwind_time_end - mp_readwind_time_beg)
    2712         mp_readwind_wtime_total = mp_readwind_wtime_total + &
    2713              &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
    2714       end if
     2722   case ('readwind')
     2723     if (imode.eq.0) then
     2724       call cpu_time(mp_readwind_time_beg)
     2725       mp_readwind_wtime_beg = mpi_wtime()
     2726     else
     2727       call cpu_time(mp_readwind_time_end)
     2728       mp_readwind_wtime_end = mpi_wtime()
     2729
     2730       mp_readwind_time_total = mp_readwind_time_total + &
     2731            &(mp_readwind_time_end - mp_readwind_time_beg)
     2732       mp_readwind_wtime_total = mp_readwind_wtime_total + &
     2733            &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
     2734     end if
    27152735
    27162736    case ('commtime')
     
    27882808          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',&
    27892809               & mp_getfields_time_total
    2790           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
    2791                & mp_readwind_wtime_total
    2792           write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
    2793                & mp_readwind_time_total
     2810!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
     2811!               & mp_readwind_wtime_total
     2812!          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
     2813!               & mp_readwind_time_total
    27942814          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',&
    27952815               & mp_io_wtime_total
  • src/netcdf_output_mod.f90

    r4ad96c5 r0a94e13  
    9393  character(len=255), parameter :: institution = 'NILU'
    9494
    95   integer            :: tpointer
     95  integer            :: tpointer=0
    9696  character(len=255) :: ncfname, ncfnamen
    9797
  • src/outg_mod.f90

    r4c64400 r2eefa58  
    3737  real,allocatable, dimension (:,:,:) :: areanorth
    3838  real,allocatable, dimension (:,:,:) :: densityoutgrid
     39  real,allocatable, dimension (:,:,:) :: densitydrygrid ! added RLT
     40  real,allocatable, dimension (:,:,:) :: factor_drygrid ! added RLT
    3941  real,allocatable, dimension (:,:,:) :: factor3d
    4042  real,allocatable, dimension (:,:,:) :: grid
  • src/outgrid_init.f90

    rd2a5a83 r2eefa58  
    268268    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
    269269  allocate(densityoutgrid(0:max(numxgrid,numxgridn)-1, &
     270       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
     271    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
     272! RLT
     273  allocate(densitydrygrid(0:max(numxgrid,numxgridn)-1, &
     274       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
     275    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
     276  allocate(factor_drygrid(0:max(numxgrid,numxgridn)-1, &
    270277       0:max(numygrid,numygridn)-1,numzgrid),stat=stat)
    271278    if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
  • src/par_mod.f90

    rdf96ea65 r2eefa58  
    3030!        Update 15 August 2013 IP                                              *
    3131!                                                                              *
    32 !        ESO 2016:                                                             *
    33 !          GFS specific parameters moved to gfs_mod.f90                        *
    34 !          ECMWF specific parameters moved to ecmwf_mod.f90                    *
    3532!                                                                              *
    3633!*******************************************************************************
     
    8077  real,parameter :: pi=3.14159265, r_earth=6.371e6, r_air=287.05, ga=9.81
    8178  real,parameter :: cpa=1004.6, kappa=0.286, pi180=pi/180., vonkarman=0.4
     79  ! additional constants RLT Aug-2017
     80  real,parameter :: rgas=8.31447
     81  real,parameter :: r_water=461.495
    8282
    8383  ! pi                      number "pi"
     
    8989  ! kappa                   exponent of formula for potential temperature
    9090  ! vonkarman               von Karman constant
     91  ! rgas                    universal gas constant [J/mol/K]
     92  ! r_water                 specific gas constant for water vapor [J/kg/K]
    9193
    9294  real,parameter :: karman=0.40, href=15., convke=2.0
     
    280282
    281283  integer,parameter :: unitpath=1, unitcommand=1, unitageclasses=1, unitgrid=1
    282   integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93
     284  integer,parameter :: unitavailab=1, unitreleases=88, unitpartout=93, unitpartout_average=105
    283285  integer,parameter :: unitpartin=93, unitflux=98, unitouttraj=96
    284286  integer,parameter :: unitvert=1, unitoro=1, unitpoin=1, unitreceptor=1
     
    290292  integer,parameter :: unitboundcond=89
    291293  integer,parameter :: unittmp=101
     294! RLT
     295  integer,parameter :: unitoutfactor=102
    292296
    293297!******************************************************
  • src/partoutput.f90

    rd2a5a83 r0a94e13  
    7171  !**************************************
    7272
    73   if (ipout.eq.1) then
     73  if (ipout.eq.1.or.ipout.eq.3) then
    7474    open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
    7575         atime,form='unformatted')
  • src/partoutput_mpi.f90

    rd2a5a83 r0a94e13  
    7878  !**************************************
    7979
    80   if (ipout.eq.1) then
     80  if (ipout.eq.1.or.ipout.eq.3) then
    8181    open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
    8282         atime,form='unformatted',status=file_stat,position='append')
  • src/readcommand.f90

    r20963b1 r2eefa58  
    5050  ! ipin                 1 continue simulation with dumped particle data, 0 no *
    5151  ! ipout                0 no particle dump, 1 every output time, 3 only at end*
     52  ! ipoutfac             increase particle dump interval by factor (default 1) *
    5253  ! itsplit [s]          time constant for particle splitting                  *
    5354  ! loutaver [s]         concentration output is an average over loutaver      *
     
    9798  iout, &
    9899  ipout, &
     100  ipoutfac, &
    99101  lsubgrid, &
    100102  lconvection, &
     
    112114  surf_only, &
    113115  cblflag, &
     116  linversionout, &
    114117  ohfields_path
    115118
     
    129132  iout=3
    130133  ipout=0
     134  ipoutfac=1
    131135  lsubgrid=1
    132136  lconvection=1
     
    144148  surf_only=0
    145149  cblflag=0 ! if using old-style COMMAND file, set to 1 here to use mc cbl routine
     150  linversionout=0
    146151  ohfields_path="../../flexin/"
    147152
     
    411416  !**********************************************************************
    412417
    413   if ((iout.lt.1).or.(iout.gt.6)) then
     418  if ((iout.lt.1).or.(iout.gt.5)) then
    414419    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    415420    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4 OR 5 FOR        #### '
     
    471476      write(*,*) '#### FOR DOMAIN FILLING RUNS OUTPUT FOR      ####'
    472477      write(*,*) '#### EACH RELEASE IS FORBIDDEN !             ####'
     478      stop
     479  endif
     480
     481  ! Inversion output format only for backward runs
     482  !*****************************************************************************
     483 
     484  if ((linversionout.eq.1).and.(ldirect.eq.1)) then
     485      write(*,*) '#### FLEXPART MODEL ERROR! FILE COMMAND:     ####'
     486      write(*,*) '#### INVERSION OUTPUT FORMAT ONLY FOR        ####'
     487      write(*,*) '#### BACKWARD RUNS                           ####'
    473488      stop
    474489  endif
     
    507522  !****************************************************************
    508523
    509   if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2)) then
     524  if ((ipout.ne.0).and.(ipout.ne.1).and.(ipout.ne.2).and.(ipout.ne.3)) then
    510525    write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
    511     write(*,*) ' #### IPOUT MUST BE 1, 2 OR 3!                #### '
     526    write(*,*) ' #### IPOUT MUST BE 0, 1, 2 OR 3!             #### '
    512527    stop
    513528  endif
  • src/readpaths.f90

    r62e65c7 r2eefa58  
    2020!**********************************************************************
    2121
    22 subroutine readpaths !(pathfile)
     22subroutine readpaths
    2323
    2424  !*****************************************************************************
  • src/readwind_gfs.f90

    rdb91eb7 r5d74ed9  
    8383
    8484  ! NCEP
    85   integer :: numpt,numpu,numpv,numpw,numprh
     85  integer :: numpt,numpu,numpv,numpw,numprh,numpclwch
    8686  real :: help, temp, ew
    8787  real :: elev
     
    134134  numpw=0
    135135  numprh=0
     136  numpclwch=0
    136137  ifield=0
    13713810   ifield=ifield+1
     
    557558      endif
    558559! SEC & IP 12/2018 read GFS clouds
    559       if(isec1(6).eq.153) then  !! CLWCR  Cloud liquid water content [kg/kg]
    560         clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     560      if((isec1(6).eq.153).and.(isec1(7).eq.100)) then  !! CLWCR  Cloud liquid water content [kg/kg]
     561         if((i.eq.0).and.(j.eq.0)) then
     562            do ii=1,nuvz
     563              if ((isec1(8)*100.0).eq.akz(ii)) numpclwch=ii
     564            end do
     565        endif
     566        help=zsec4(nxfield*(ny-j-1)+i+1)
     567        if(i.le.i180) then
     568          clwch(i179+i,j,numpclwch,n)=help
     569        else
     570          clwch(i-i181,j,numpclwch,n)=help
     571        endif
    561572        readclouds=.true.
    562573        sumclouds=.true.
     574!        readclouds=.false.
     575!       sumclouds=.false.
    563576      endif
    564577
  • src/releaseparticles.f90

    r75a4ded r7873bf7  
    114114      average_timecorrect=0.
    115115      do k=1,nspec
    116         if (zpoint1(i).gt.0.5) then      ! point source
     116        if(abs(xpoint2(i)-xpoint1(i)).lt.1.E-4.and.abs(ypoint2(i)-ypoint1(i)).lt.1.E-4) then
     117!        if (zpoint1(i).gt.0.5) then      ! point source
    117118          timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
    118119        else                             ! area source
  • src/timemanager.f90

    rc7d1052 r2eefa58  
    408408#endif
    409409            else
     410              if (linversionout.eq.1) then
     411                call concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     412                if (verbosity.eq.1) then
     413                  print*,'called concoutput_inversion'
     414                  call system_clock(count_clock)
     415                  write(*,*) 'system clock',count_clock - count_clock0
     416                endif
     417              else
    410418              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     419              endif
    411420              if (verbosity.eq.1) then
    412421                print*,'called concoutput_surf '
     
    422431                call concoutput_nest(itime,outnum)
    423432              else
     433                if(linversionout.eq.1) then
     434                  call concoutput_inversion_nest(itime,outnum)
     435                else
    424436                call concoutput_surf_nest(itime,outnum)
     437              endif
    425438              endif
    426439            else
     
    45146445      format(i13,' Seconds simulated: ',i13, ' Particles:    Uncertainty: ',3f7.3)
    45246546      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
    453         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
     466        if (ipout.ge.1) then
     467          if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
     468          if (ipout.eq.3) call partoutput_average(itime) ! dump particle positions
     469        endif
    454470        loutnext=loutnext+loutstep
    455471        loutstart=loutnext-loutaver/2
     
    609625!        write (*,*) 'advance: ',prob(1),xmass1(j,1),ztra1(j)
    610626
     627  ! Calculate average position for particle dump output
     628  !****************************************************
     629
     630        if (ipout.eq.3) call partpos_average(itime,j)
     631
     632
    611633  ! Calculate the gross fluxes across layer interfaces
    612634  !***************************************************
     
    730752  if (ipout.eq.2) call partoutput(itime)     ! dump particle positions
    731753
    732   if (linit_cond.ge.1) call initial_cond_output(itime)   ! dump initial cond. field
     754  if (linit_cond.ge.1) then
     755    if(linversionout.eq.1) then
     756      call initial_cond_output_inversion(itime)   ! dump initial cond. field
     757    else
     758      call initial_cond_output(itime)   ! dump initial cond. fielf
     759    endif
     760  endif
    733761
    734762  !close(104)
  • src/timemanager_mpi.f90

    r20963b1 r0c8c7f2  
    113113  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0
    114114! integer :: ksp
    115   integer :: ip
     115  integer :: ip,irec
    116116  integer :: loutnext,loutstart,loutend
    117117  integer :: ix,jy,ldeltat,itage,nage,idummy
     
    129129! Measure time spent in timemanager
    130130  if (mp_measure_time) call mpif_mtime('timemanager',0)
     131
    131132
    132133! First output for time 0
     
    532533                end if
    533534
    534               else  ! :TODO: check for zeroing in the netcdf module
     535              else
    535536                call concoutput_surf_nest(itime,outnum)
    536537              end if
     
    59359446      format(' Simulated ',f7.1,' hours (',i13,' s), ',i13, ' particles')
    594595        if (ipout.ge.1) then
     596          if (mp_measure_time) call mpif_mtime('iotime',0)
     597          irec=0
    595598          do ip=0, mp_partgroup_np-1
    596             if (ip.eq.mp_partid) call partoutput(itime) ! dump particle positions
     599            if (ip.eq.mp_partid) then
     600              if (mod(itime,ipoutfac*loutstep).eq.0) call partoutput(itime) ! dump particle positions
     601              if (ipout.eq.3) call partoutput_average(itime,irec) ! dump particle positions
     602            endif
     603            if (ipout.eq.3) irec=irec+npart_per_process(ip)
    597604            call mpif_mpi_barrier
    598605          end do
     606          if (mp_measure_time) call mpif_mtime('iotime',1)
    599607        end if
    600608
     
    757765        if (mp_measure_time) call mpif_mtime('advance',1)
    758766
     767  ! Calculate average position for particle dump output
     768  !****************************************************
     769
     770        if (ipout.eq.3) call partpos_average(itime,j)
     771
    759772
    760773! Calculate the gross fluxes across layer interfaces
     
    895908    do ip=0, mp_partgroup_np-1
    896909      if (ip.eq.mp_partid) then
    897         !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
     910        if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
    898911        call partoutput(itime)    ! dump particle positions
    899912      end if
  • src/verttransform_ecmwf.f90

    r79e0349 r2eefa58  
    7373  use com_mod
    7474  use cmapf_mod, only: cc2gll
    75 !  use mpi_mod
    7675
    7776  implicit none
     
    8281  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
    8382  real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
     83  ! RLT added pressure
     84  real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: prsh
    8485  real,dimension(0:nxmax-1,0:nymax-1) ::  tvold,pold,pint,tv
    8586  real,dimension(0:nymax-1) :: cosf
     
    219220!*************************
    220221
    221 
    222222  do jy=0,nymin1
    223223    do ix=0,nxmin1
     
    230230  wzlev(:,:,1)=0.
    231231  rhoh(:,:,1)=pold/(r_air*tvold)
     232  ! RLT add pressure
     233  prsh(:,:,1)=ps(:,:,1,n)
    232234
    233235
     
    237239  do kz=2,nuvz
    238240    pint=akz(kz)+bkz(kz)*ps(:,:,1,n)
     241    ! RLT add pressure
     242    prsh(:,:,kz)=pint
    239243    tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n))
    240244    rhoh(:,:,kz)=pint(:,:)/(r_air*tv)
     
    288292  pv(:,:,1,n)=pvh(:,:,1)
    289293  rho(:,:,1,n)=rhoh(:,:,1)
     294! RLT add pressure
     295  prs(:,:,1,n)=prsh(:,:,1)
    290296
    291297  uu(:,:,nz,n)=uuh(:,:,nuvz)
     
    301307  pv(:,:,nz,n)=pvh(:,:,nuvz)
    302308  rho(:,:,nz,n)=rhoh(:,:,nuvz)
    303 
     309! RLT
     310  prs(:,:,nz,n)=prsh(:,:,nuvz)
    304311
    305312  kmin=2
     
    321328          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
    322329          rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
     330! RLT
     331          prs(ix,jy,iz,n)=prs(ix,jy,nz,n)
    323332        else
    324333          innuvz: do kz=idx(ix,jy),nuvz
     
    354363          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
    355364          rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz
     365! RLT add pressure
     366          prs(ix,jy,iz,n)=(prsh(ix,jy,kz-1)*dz2+prsh(ix,jy,kz)*dz1)/dz
    356367        endif
    357368      enddo
     
    654665        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    655666
    656           do kz=nz,1,-1 !go Bottom up!
     667          do kz=nz,2,-1 !go Bottom up!
    657668            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
    658669              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
  • src/verttransform_gfs.f90

    rdb91eb7 r437c545  
    548548        if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
    549549
    550           do kz=nz,1,-1 !go Bottom up!
     550          do kz=nz,2,-1 !go Bottom up!
    551551            if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
    552552              cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
    553553              clouds(ix,jy,kz,n)=1                               ! is a cloud
    554554              if (lsp.ge.convp) then
    555                 clouds(ix,jy,kz,n)=3                            ! lsp in-cloud
     555                clouds(ix,jy,kz,n)=3                             ! lsp in-cloud
    556556              else
    557557                clouds(ix,jy,kz,n)=2                             ! convp in-cloud
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG