Changeset 5f9d14a in flexpart.git for src


Ignore:
Timestamp:
Apr 8, 2015, 2:23:27 PM (9 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
1585284
Parents:
cd85138
Message:

Updated wet depo scheme

Location:
src
Files:
53 edited

Legend:

Unmodified
Added
Removed

  • src/FLEXPART_MPI.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    424424
    425425! NIK 16.02.2005
    426   write(*,*) '**********************************************'
    427   write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
    428   write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
    429   write(*,*) '**********************************************'
    430 
    431426  if (lroot) then
     427! eso TODO: do MPI_Reduce (sum) of total occurences across processes
     428    call MPI_Reduce(MPI_IN_PLACE, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
     429         & mp_comm_used, mp_ierr)
     430    call MPI_Reduce(MPI_IN_PLACE, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
     431         & mp_comm_used, mp_ierr)
     432  else
     433    call MPI_Reduce(tot_blc_count, tot_blc_count, 1, mp_pp, MPI_SUM, id_root, &
     434         & mp_comm_used, mp_ierr)
     435    call MPI_Reduce(tot_inc_count, tot_inc_count, 1, mp_pp, MPI_SUM, id_root, &
     436         & mp_comm_used, mp_ierr)
     437  end if
     438
     439  if (lroot) then
     440    write(*,*) '**********************************************'
     441    write(*,*) 'Total number of occurences of below-cloud scavenging', tot_blc_count
     442    write(*,*) 'Total number of occurences of in-cloud    scavenging', tot_inc_count
     443    write(*,*) '**********************************************'
     444
    432445    write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
    433446         &XPART MODEL RUN!'

  • src/advance.f90

    r8a65cb0 r5f9d14a  
    229229  if (ngrid.le.0) then
    230230    do k=1,2
    231 ! eso: compatibility with 3-field version
    232       mind=memind(k)
     231      mind=memind(k) ! eso: compatibility with 3-field version
    233232      do j=jy,jyp
    234233        do i=ix,ixp

  • src/boundcond_domainfill_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/cbl.f90

    • Property mode changed from 100755 to 100644

  • src/com_mod.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    269269  !                         disc or on tape
    270270
    271   integer :: memtime(numwfmem),memind(numwfmem) ! eso: or memind(3) and change
    272                                                 ! interpol_rain
     271  integer :: memtime(numwfmem),memind(3) ! eso: or memind(numwfmem)
    273272
    274273  ! memtime [s]             validation times of wind fields in memory

  • src/conccalc_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    6262  real :: xl,yl,wx,wy,w
    6363  real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
    64 
     64  integer :: mind2
     65  ! mind2        eso: pointer to 2nd windfield in memory
    6566
    6667  ! For forward simulations, make a loop over the number of species;
    … …  
    7071
    7172  if (mp_measure_time) call mpif_mtime('conccalc',0)
     73
     74  mind2=memind(2)
    7275
    7376  do i=1,numpart
    … …  
    131134  !*****************************************************************************
    132135      do ind=indz,indzp
    133         rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
    134              +p2*rho(ixp,jy ,ind,2) &
    135              +p3*rho(ix ,jyp,ind,2) &
    136              +p4*rho(ixp,jyp,ind,2)
     136        ! rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
     137        !      +p2*rho(ixp,jy ,ind,2) &
     138        !      +p3*rho(ix ,jyp,ind,2) &
     139        !      +p4*rho(ixp,jyp,ind,2)
     140        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
     141             +p2*rho(ixp,jy ,ind,mind2) &
     142             +p3*rho(ix ,jyp,ind,mind2) &
     143             +p4*rho(ixp,jyp,ind,mind2)
     144
    137145      end do
    138146      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz

  • src/concoutput.f90

    rf13406c r5f9d14a  
    9999  character :: adate*8,atime*6
    100100  character(len=3) :: anspec
     101  integer :: mind
     102! mind        eso:added to ensure identical results between 2&3-fields versions
    101103
    102104
    … …  
    140142  !*******************************************************************
    141143
     144  mind=memind(2)
    142145  do kz=1,numzgrid
    143146    if (kz.eq.1) then
    … …  
    162165        iix=max(min(nint(xl),nxmin1),0)
    163166        jjy=max(min(nint(yl),nymin1),0)
    164         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    165              rho(iix,jjy,kzz-1,2)*dz2)/dz
     167        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     168        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     169        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     170             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    166171      end do
    167172    end do
    168173  end do
    169174
    170     do i=1,numreceptor
    171       xl=xreceptor(i)
    172       yl=yreceptor(i)
    173       iix=max(min(nint(xl),nxmin1),0)
    174       jjy=max(min(nint(yl),nymin1),0)
    175       densityoutrecept(i)=rho(iix,jjy,1,2)
    176     end do
     175  do i=1,numreceptor
     176    xl=xreceptor(i)
     177    yl=yreceptor(i)
     178    iix=max(min(nint(xl),nxmin1),0)
     179    jjy=max(min(nint(yl),nymin1),0)
     180    !densityoutrecept(i)=rho(iix,jjy,1,2)
     181    densityoutrecept(i)=rho(iix,jjy,1,mind)
     182  end do
    177183
    178184

  • src/concoutput_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    104104  character :: adate*8,atime*6
    105105  character(len=3) :: anspec
     106  integer :: mind
     107! mind        eso:added to ensure identical results between 2&3-fields versions
    106108
    107109! Measure execution time
    … …  
    148150!*******************************************************************
    149151
     152  mind=memind(2)
    150153  do kz=1,numzgrid
    151154    if (kz.eq.1) then
    … …  
    170173        iix=max(min(nint(xl),nxmin1),0)
    171174        jjy=max(min(nint(yl),nymin1),0)
    172         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    173              rho(iix,jjy,kzz-1,2)*dz2)/dz
     175        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     176        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     177        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     178             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    174179      end do
    175180    end do
    … …  
    181186    iix=max(min(nint(xl),nxmin1),0)
    182187    jjy=max(min(nint(yl),nymin1),0)
    183     densityoutrecept(i)=rho(iix,jjy,1,2)
     188    !densityoutrecept(i)=rho(iix,jjy,1,2)
     189    densityoutrecept(i)=rho(iix,jjy,1,mind)
    184190  end do
    185191

  • src/concoutput_nest.f90

    re200b7a r5f9d14a  
    9595  character :: adate*8,atime*6
    9696  character(len=3) :: anspec
     97  integer :: mind
     98! mind        eso:added to ensure identical results between 2&3-fields versions
    9799
    98100
    … …  
    134136  !*******************************************************************
    135137
     138  mind=memind(2)
    136139  do kz=1,numzgrid
    137140    if (kz.eq.1) then
    … …  
    156159        iix=max(min(nint(xl),nxmin1),0)
    157160        jjy=max(min(nint(yl),nymin1),0)
    158         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    159              rho(iix,jjy,kzz-1,2)*dz2)/dz
     161        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     162        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     163        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     164             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    160165      end do
    161166    end do
    162167  end do
    163168
    164     do i=1,numreceptor
    165       xl=xreceptor(i)
    166       yl=yreceptor(i)
    167       iix=max(min(nint(xl),nxmin1),0)
    168       jjy=max(min(nint(yl),nymin1),0)
    169       densityoutrecept(i)=rho(iix,jjy,1,2)
    170     end do
     169  do i=1,numreceptor
     170    xl=xreceptor(i)
     171    yl=yreceptor(i)
     172    iix=max(min(nint(xl),nxmin1),0)
     173    jjy=max(min(nint(yl),nymin1),0)
     174    !densityoutrecept(i)=rho(iix,jjy,1,2)
     175    densityoutrecept(i)=rho(iix,jjy,1,mind)
     176  end do
    171177
    172178

  • src/concoutput_nest_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    9999  character :: adate*8,atime*6
    100100  character(len=3) :: anspec
     101  integer :: mind
     102! mind        eso:added to ensure identical results between 2&3-fields versions
    101103
    102104! Measure execution time
    … …  
    149151  !*******************************************************************
    150152
     153  mind=memind(2)
    151154  do kz=1,numzgrid
    152155    if (kz.eq.1) then
    … …  
    171174        iix=max(min(nint(xl),nxmin1),0)
    172175        jjy=max(min(nint(yl),nymin1),0)
    173         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    174              rho(iix,jjy,kzz-1,2)*dz2)/dz
     176        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     177        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     178        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     179             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    175180      end do
    176181    end do
    … …  
    182187    iix=max(min(nint(xl),nxmin1),0)
    183188    jjy=max(min(nint(yl),nymin1),0)
    184     densityoutrecept(i)=rho(iix,jjy,1,2)
     189    !densityoutrecept(i)=rho(iix,jjy,1,2)
     190    densityoutrecept(i)=rho(iix,jjy,1,mind)
    185191  end do
    186192

  • src/concoutput_surf_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    105105  character :: adate*8,atime*6
    106106  character(len=3) :: anspec
     107  integer :: mind
     108! mind        eso:added to get consistent results between 2&3-fields versions
    107109
    108110! Measure execution time
    109   if (mp_measure_time) then
    110     call cpu_time(mp_root_time_beg)
    111     mp_root_wtime_beg = mpi_wtime()
    112   end if
    113 
    114   if (verbosity.eq.1) then
    115      print*,'inside concoutput_surf '
    116      CALL SYSTEM_CLOCK(count_clock)
    117      WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
    118   endif
     111  if (mp_measure_time) call mpif_mtime('rootonly',0)
     112
    119113
    120114  ! Determine current calendar date, needed for the file name
    … …  
    165159  !*******************************************************************
    166160
     161  mind=memind(2)
    167162  do kz=1,numzgrid
    168163    if (kz.eq.1) then
    … …  
    187182        iix=max(min(nint(xl),nxmin1),0)
    188183        jjy=max(min(nint(yl),nymin1),0)
    189         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    190              rho(iix,jjy,kzz-1,2)*dz2)/dz
     184        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     185        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     186        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     187             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    191188      end do
    192189    end do
    193190  end do
    194191
    195     do i=1,numreceptor
    196       xl=xreceptor(i)
    197       yl=yreceptor(i)
    198       iix=max(min(nint(xl),nxmin1),0)
    199       jjy=max(min(nint(yl),nymin1),0)
    200       densityoutrecept(i)=rho(iix,jjy,1,2)
    201     end do
     192  do i=1,numreceptor
     193    xl=xreceptor(i)
     194    yl=yreceptor(i)
     195    iix=max(min(nint(xl),nxmin1),0)
     196    jjy=max(min(nint(yl),nymin1),0)
     197    !densityoutrecept(i)=rho(iix,jjy,1,2)
     198    densityoutrecept(i)=rho(iix,jjy,1,mind)
     199  end do
    202200
    203201
    … …  
    655653  end do
    656654
    657   if (mp_measure_time) then
    658     call cpu_time(mp_root_time_end)
    659     mp_root_wtime_end = mpi_wtime()
    660     mp_root_time_total = mp_root_time_total + (mp_root_time_end - mp_root_time_beg)
    661     mp_root_wtime_total = mp_root_wtime_total + (mp_root_wtime_end - mp_root_wtime_beg)
    662   end if
    663 
     655  if (mp_measure_time) call mpif_mtime('rootonly',1)
     656 
    664657end subroutine concoutput_surf

  • src/concoutput_surf_nest_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    9898  character :: adate*8,atime*6
    9999  character(len=3) :: anspec
     100  integer :: mind
     101! mind        eso:added to get consistent results between 2&3-fields versions
    100102
    101103! Measure execution time
    … …  
    148150  !*******************************************************************
    149151
     152  mind=memind(2)
    150153  do kz=1,numzgrid
    151154    if (kz.eq.1) then
    … …  
    170173        iix=max(min(nint(xl),nxmin1),0)
    171174        jjy=max(min(nint(yl),nymin1),0)
    172         densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
    173              rho(iix,jjy,kzz-1,2)*dz2)/dz
     175        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
     176        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
     177        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
     178             rho(iix,jjy,kzz-1,mind)*dz2)/dz
    174179      end do
    175180    end do
    176181  end do
    177182
    178     do i=1,numreceptor
    179       xl=xreceptor(i)
    180       yl=yreceptor(i)
    181       iix=max(min(nint(xl),nxmin1),0)
    182       jjy=max(min(nint(yl),nymin1),0)
    183       densityoutrecept(i)=rho(iix,jjy,1,2)
    184     end do
     183  do i=1,numreceptor
     184    xl=xreceptor(i)
     185    yl=yreceptor(i)
     186    iix=max(min(nint(xl),nxmin1),0)
     187    jjy=max(min(nint(yl),nymin1),0)
     188    !densityoutrecept(i)=rho(iix,jjy,1,2)
     189    densityoutrecept(i)=rho(iix,jjy,1,mind)
     190  end do
    185191
    186192

  • src/erf.f90

    • Property mode changed from 100755 to 100644

  • src/getfields_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    149149! 3 fields in memory
    150150! Note that the reader process never needs to keep 3 fields in memory,
    151 ! so it is possible to save some memory here
     151! (2 is enough) so it is possible to save some memory
    152152!*********************************************************************
    153153      if (mind3.eq.32.or.mind3.eq.19) then
    … …  
    245245        end if
    246246        memtime(2)=wftime(indj+1)
    247 ! :DEV: not used
     247! DEV: not used
    248248        if (.not.lmp_sync) memind(3)=3 ! indicate position of next read
    249249        nstop = 1

  • src/gethourlyOH.f90

    • Property mode changed from 100755 to 100644

  • src/getvdep.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_emos.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_gfs.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_gfs_emos.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_nests.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_nests_emos.f90

    • Property mode changed from 100755 to 100644

  • src/gridcheck_orig_ecmwf.f90

    • Property mode changed from 100755 to 100644

  • src/init_domainfill_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/initial_cond_calc.f90

    re200b7a r5f9d14a  
    4444  real :: ddx,ddy
    4545  real :: rhoprof(2),rhoi,xl,yl,wx,wy,w
     46  integer :: mind2
     47  ! mind2        eso: pointer to 2nd windfield in memory
    4648
    4749
    … …  
    8890  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
    8991  !*****************************************************************************
     92    mind2=memind(2)
     93
    9094    do ind=indz,indzp
    91       rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
    92            +p2*rho(ixp,jy ,ind,2) &
    93            +p3*rho(ix ,jyp,ind,2) &
    94            +p4*rho(ixp,jyp,ind,2)
     95      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
     96           +p2*rho(ixp,jy ,ind,mind2) &
     97           +p3*rho(ix ,jyp,ind,mind2) &
     98           +p4*rho(ixp,jyp,ind,mind2)
    9599    end do
    96100    rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz

  • src/initialize_cbl_vel.f90

    • Property mode changed from 100755 to 100644

  • src/makefile

    rcd85138 r5f9d14a  
    6262
    6363FFLAGS   = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) -mtune=native $(FUSER) # -march=native
    64 DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -fdump-core -Warray-bounds -ffpe-trap=invalid,overflow,denormal -Wall -fcheck=all $(FUSER)  # ,underflow,zero
     64DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Warray-bounds  -Wall -fcheck=all $(FUSER)  # -ffpe-trap=invalid,overflow,denormal,underflow,zero -fdump-core
    6565
    6666LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -L$(LIBPATH2) $(LIBS)

  • src/mpi_mod.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    4444! mp_pid                  Process ID of each MPI process                     *
    4545! mp_seed                 Parameter for random number seed                   *
     46! read_grp_min            Minimum number of processes at which one will be   *
     47!                         used as reader                                     *
    4648! numpart_mpi,            Number of particles per node                       *
    4749! maxpart_mpi                                                                *
    … …  
    5052!                         loops over particles. Will be all processes        *
    5153!                         unless a dedicated process runs getfields/readwind *
    52 !                                                                            *
     54! lmp_sync                If .false., use asynchronous MPI                   *
    5355!                                                                            *
    5456!                                                                            *
    … …  
    6971! Set aside a process for reading windfields if using at least these many processes
    7072!==================================================
    71   integer, parameter, private :: read_grp_min=99
     73  integer, parameter, private :: read_grp_min=4
    7274!==================================================
    7375
    … …  
    8688
    8789! MPI tags/requests for send/receive operation
    88   integer :: tm1!=100
    89   integer, parameter :: nvar_async=27
     90  integer :: tm1
     91  integer, parameter :: nvar_async=27 !29 :DBG:
    9092  !integer, dimension(:), allocatable :: tags
    9193  integer, dimension(:), allocatable :: reqs
    … …  
    106108!===============================================================================
    107109
    108 ! mp_dbg_mode       MPI related output for debugging etc.
     110! mp_dbg_mode       Used for debugging MPI.
    109111! mp_dev_mode       various checks related to debugging the parallel code
    110112! mp_dbg_out        write some arrays to data file for debugging
    … …  
    197199           & 'numwfmem should be set to 2 for syncronous reading'
    198200      write(*,FMT='(80("#"))')
    199     end if
     201! Force "syncronized" version if all processes will call getfields
     202    else if (.not.lmp_sync.and.mp_np.lt.read_grp_min) then
     203      if (lroot) then
     204        write(*,FMT='(80("#"))')
     205        write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
     206             & 'all procs call getfields. Setting lmp_sync=.true.'
     207        write(*,FMT='(80("#"))')
     208      end if
     209      lmp_sync=.true. ! :DBG: eso fix this...
     210    end if
     211! TODO: Add warnings for unimplemented flexpart features
    200212
    201213! Set ID of process that calls getfield/readwind.
    … …  
    211223    call MPI_Comm_group (MPI_COMM_WORLD, world_group_id, mp_ierr)
    212224
    213 ! Create a MPI group/communiactor that will calculate trajectories.
     225! Create a MPI group/communicator that will calculate trajectories.
    214226! Skip this step if program is run with only a few processes
    215227!************************************************************************
    … …  
    228240      lmp_use_reader = .true.
    229241
    230 ! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping 'readwind' process
     242! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping reader process
    231243      j=-1
    232244      do i=0, mp_np-2 !loop over all (subgroup) IDs
    … …  
    289301    end if
    290302    if (mp_dev_mode) write(*,*) 'PID, mp_seed: ',mp_pid, mp_seed
     303    if (mp_dbg_mode) then
     304! :DBG: For debugging, set all seed to 0 and maxrand to e.g. 20
     305      mp_seed=0
     306      if (lroot) write(*,*) 'MPI: setting seed=0'
     307    end if
    291308
    292309! Allocate request array for asynchronous MPI
    293310    if (.not.lmp_sync) then
    294311      allocate(reqs(0:nvar_async*mp_np-1))
    295       reqs=MPI_REQUEST_NULL
     312      reqs(:)=MPI_REQUEST_NULL
    296313    else
    297314      allocate(reqs(0:1))
    298       reqs=MPI_REQUEST_NULL
     315      reqs(:)=MPI_REQUEST_NULL
    299316    end if
    300317
    … …  
    389406    integer :: add_part=0
    390407
    391     call MPI_BCAST(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
     408    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
    392409
    393410! MPI subgroup does the particle-loop
    … …  
    815832!*******************************************************************************
    816833! DESCRIPTION
    817 !   Distribute 'getfield' variables from reader process to all processes.
     834!   Distribute 'getfield' variables from reader process
    818835!
    819836!   Called from timemanager
    … …  
    887904!**********************************************************************
    888905
     906! The non-reader processes need to know if clouds were read.
     907! TODO: only at first step or always?
     908    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
     909    if (mp_ierr /= 0) goto 600
     910
    889911! Static fields/variables sent only at startup
    890912    if (first_call) then
    … …  
    904926      if (mp_ierr /= 0) goto 600
    905927      call MPI_Bcast(height,nzmax,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
    906     if (mp_ierr /= 0) goto 600
     928      if (mp_ierr /= 0) goto 600
    907929
    908930      first_call=.false.
    … …  
    934956    if (mp_ierr /= 0) goto 600
    935957    call MPI_Bcast(clouds(:,:,:,li:ui),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
    936     if (mp_ierr /= 0) goto 600
     958    if (mp_ierr /= 0) goto 600
     959
     960! cloud water/ice:
     961    if (readclouds) then
     962      call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
     963      if (mp_ierr /= 0) goto 600
     964      call MPI_Bcast(ciwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
     965      if (mp_ierr /= 0) goto 600
     966    end if
    937967
    938968! 2D fields
    … …  
    969999    if (mp_ierr /= 0) goto 600
    9701000
    971 
    9721001    if (mp_measure_time) call mpif_mtime('commtime',1)
    9731002
    … …  
    9951024!   step
    9961025!
     1026! TODO
     1027!   Transfer cloud water/ice if and when available for nested
    9971028!
    9981029!***********************************************************************
    … …  
    11401171! DESCRIPTION
    11411172!   Distribute 'getfield' variables from reader process to all processes.
    1142 !
    1143 !   Called from timemanager by root process only
     1173!   Called from timemanager by reader process only
    11441174!
    11451175! NOTE
    … …  
    11541184!   memstat -- input, for resolving pointer to windfield index being read
    11551185!   mind    -- index where to place new fields
    1156 !   !!! Under development, don't use yet !!!
    1157 !
     1186!
     1187! TODO
     1188!   Transfer cloud water/ice
    11581189!
    11591190!*******************************************************************************
    … …  
    11731204    integer :: d1s1 = maxwf
    11741205
    1175     !integer :: d3s1,d3s2,d2s1,d2s2
    1176 
    11771206!*******************************************************************************
    1178 
    1179 ! :TODO: don't need these
    1180     ! d3s1=d3_size1
    1181     ! d3s2=d3_size2
    1182     ! d2s1=d2_size1
    1183     ! d2s2=d2_size2
    11841207
    11851208! At the time the send is posted, the reader process is one step further
    … …  
    12201243
    12211244    do dest=0,mp_np-1 ! mp_np-2 will also work if last proc reserved for reading
    1222                       ! :TODO: use mp_partgroup_np here
     1245                      ! TODO: use mp_partgroup_np here
    12231246      if (dest.eq.id_read) cycle
    12241247      i=dest*nvar_async
    … …  
    13051328      if (mp_ierr /= 0) goto 600
    13061329
     1330! Send cloud water if it exists. Increment counter always (as on receiving end)
     1331      if (readclouds) then
     1332        i=i+1
     1333        call MPI_Isend(clwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
     1334             &MPI_COMM_WORLD,reqs(i),mp_ierr)
     1335        if (mp_ierr /= 0) goto 600
     1336        i=i+1
     1337
     1338        call MPI_Isend(ciwc(:,:,:,mind),d3s1,mp_pp,dest,tm1,&
     1339             &MPI_COMM_WORLD,reqs(i),mp_ierr)
     1340        if (mp_ierr /= 0) goto 600
     1341      ! else
     1342      !   i=i+2
     1343      end if
     1344
    13071345    end do
    13081346
    … …  
    13111349    goto 601
    13121350
    1313 600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
     1351600 write(*,*) "#### mpi_mod::mpif_gf_send_vars_async> mp_ierr \= 0", mp_ierr
    13141352    stop
    13151353
    … …  
    13201358!*******************************************************************************
    13211359! DESCRIPTION
    1322 !   Receive 'getfield' variables from reader process to all processes.
    1323 !
    1324 !   Called from timemanager by all processes except root
     1360!   Receive 'getfield' variables from reader process.
     1361!   Called from timemanager by all processes except reader
    13251362!
    13261363! NOTE
    … …  
    13301367!   memstat -- input, used to resolve windfield index being received
    13311368!
     1369! TODO
     1370!   Transfer cloud water/ice
    13321371!
    13331372!*******************************************************************************
    … …  
    13881427! Get MPI tags/requests for communications
    13891428    j=mp_pid*nvar_async
    1390     call MPI_Irecv(uu(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1391     if (mp_ierr /= 0) goto 600
    1392     j=j+1
    1393     call MPI_Irecv(vv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1394     if (mp_ierr /= 0) goto 600
    1395     j=j+1
    1396     call MPI_Irecv(uupol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1397     if (mp_ierr /= 0) goto 600
    1398     j=j+1
    1399     call MPI_Irecv(vvpol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1400     if (mp_ierr /= 0) goto 600
    1401     j=j+1
    1402     call MPI_Irecv(ww(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1429    call MPI_Irecv(uu(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1430         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1431    if (mp_ierr /= 0) goto 600
     1432    j=j+1
     1433    call MPI_Irecv(vv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1434         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1435    if (mp_ierr /= 0) goto 600
     1436    j=j+1
     1437    call MPI_Irecv(uupol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1438         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1439    if (mp_ierr /= 0) goto 600
     1440    j=j+1
     1441    call MPI_Irecv(vvpol(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1442         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1443    if (mp_ierr /= 0) goto 600
     1444    j=j+1
     1445    call MPI_Irecv(ww(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1446         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14031447    if (mp_ierr /= 0) goto 600
    14041448    j=j+1
    1405     call MPI_Irecv(tt(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1406     if (mp_ierr /= 0) goto 600
    1407     j=j+1
    1408     call MPI_Irecv(rho(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1409     if (mp_ierr /= 0) goto 600
    1410     j=j+1
    1411     call MPI_Irecv(drhodz(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1412     if (mp_ierr /= 0) goto 600
    1413     j=j+1
    1414     call MPI_Irecv(tth(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1415     if (mp_ierr /= 0) goto 600
    1416     j=j+1
    1417     call MPI_Irecv(qvh(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1418     if (mp_ierr /= 0) goto 600
    1419     j=j+1
    1420 
    1421     call MPI_Irecv(qv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1449    call MPI_Irecv(tt(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1450         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1451    if (mp_ierr /= 0) goto 600
     1452    j=j+1
     1453    call MPI_Irecv(rho(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1454         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1455    if (mp_ierr /= 0) goto 600
     1456    j=j+1
     1457    call MPI_Irecv(drhodz(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1458         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1459    if (mp_ierr /= 0) goto 600
     1460    j=j+1
     1461    call MPI_Irecv(tth(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,&
     1462         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1463    if (mp_ierr /= 0) goto 600
     1464    j=j+1
     1465    call MPI_Irecv(qvh(:,:,:,mind),d3s2,mp_pp,id_read,MPI_ANY_TAG,&
     1466         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1467    if (mp_ierr /= 0) goto 600
     1468    j=j+1
     1469
     1470    call MPI_Irecv(qv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1471         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14221472    if (mp_ierr /= 0) goto 600
    1423 j=j+1
    1424     call MPI_Irecv(pv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1425     if (mp_ierr /= 0) goto 600
    1426     j=j+1
    1427     call MPI_Irecv(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)   
    1428     if (mp_ierr /= 0) goto 600
    1429     j=j+1
    1430 
    1431     call MPI_Irecv(cloudsh(:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1432     if (mp_ierr /= 0) goto 600
    1433     j=j+1
    1434     call MPI_Irecv(vdep(:,:,:,mind),d2s2,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1435     if (mp_ierr /= 0) goto 600
    1436     j=j+1
    1437 
    1438     call MPI_Irecv(ps(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1473    j=j+1
     1474    call MPI_Irecv(pv(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1475         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1476    if (mp_ierr /= 0) goto 600
     1477    j=j+1
     1478    call MPI_Irecv(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,&
     1479         &MPI_COMM_WORLD,reqs(j),mp_ierr)   
     1480    if (mp_ierr /= 0) goto 600
     1481    j=j+1
     1482
     1483    call MPI_Irecv(cloudsh(:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1484         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1485    if (mp_ierr /= 0) goto 600
     1486    j=j+1
     1487    call MPI_Irecv(vdep(:,:,:,mind),d2s2,mp_pp,id_read,MPI_ANY_TAG,&
     1488         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1489    if (mp_ierr /= 0) goto 600
     1490    j=j+1
     1491    call MPI_Irecv(ps(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1492         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14391493    if (mp_ierr /= 0) goto 600
    14401494    j=j+1
    1441     call MPI_Irecv(sd(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1495    call MPI_Irecv(sd(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1496         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14421497    if (mp_ierr /= 0) goto 600
    14431498    j=j+1
    1444     call MPI_Irecv(tcc(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1499    call MPI_Irecv(tcc(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1500         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14451501    if (mp_ierr /= 0) goto 600
    14461502    j=j+1
    1447     call MPI_Irecv(tt2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1448     if (mp_ierr /= 0) goto 600
    1449     j=j+1
    1450     call MPI_Irecv(td2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1451     if (mp_ierr /= 0) goto 600
    1452     j=j+1
    1453     call MPI_Irecv(lsprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1454     if (mp_ierr /= 0) goto 600
    1455     j=j+1
    1456     call MPI_Irecv(convprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1457     if (mp_ierr /= 0) goto 600
    1458 
    1459     call MPI_Irecv(ustar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1460     if (mp_ierr /= 0) goto 600
    1461     j=j+1
    1462     call MPI_Irecv(wstar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1463     if (mp_ierr /= 0) goto 600
    1464     j=j+1
    1465     call MPI_Irecv(hmix(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1466     if (mp_ierr /= 0) goto 600
    1467     j=j+1
    1468     call MPI_Irecv(tropopause(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
     1503    call MPI_Irecv(tt2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1504         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1505    if (mp_ierr /= 0) goto 600
     1506    j=j+1
     1507    call MPI_Irecv(td2(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1508         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1509    if (mp_ierr /= 0) goto 600
     1510    j=j+1
     1511    call MPI_Irecv(lsprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1512         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1513    if (mp_ierr /= 0) goto 600
     1514    j=j+1
     1515    call MPI_Irecv(convprec(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1516         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1517    if (mp_ierr /= 0) goto 600
     1518
     1519    call MPI_Irecv(ustar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1520         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1521    if (mp_ierr /= 0) goto 600
     1522    j=j+1
     1523    call MPI_Irecv(wstar(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1524         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1525    if (mp_ierr /= 0) goto 600
     1526    j=j+1
     1527    call MPI_Irecv(hmix(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1528         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1529    if (mp_ierr /= 0) goto 600
     1530    j=j+1
     1531    call MPI_Irecv(tropopause(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1532         &MPI_COMM_WORLD,reqs(j),mp_ierr)
    14691533    if (mp_ierr /= 0) goto 600
    14701534    j=j+1
    1471     call MPI_Irecv(oli(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,MPI_COMM_WORLD,reqs(j),mp_ierr)
    1472     if (mp_ierr /= 0) goto 600
     1535    call MPI_Irecv(oli(:,:,:,mind),d2s1,mp_pp,id_read,MPI_ANY_TAG,&
     1536         &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1537    if (mp_ierr /= 0) goto 600
     1538
     1539
     1540! Post request for clwc. These data are possibly not sent, request must then be cancelled
     1541! For now assume that data at all steps either have or do not have water
     1542    if (readclouds) then
     1543      j=j+1
     1544      call MPI_Irecv(clwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1545           &MPI_COMM_WORLD,reqs(j),mp_ierr)   
     1546      if (mp_ierr /= 0) goto 600
     1547      j=j+1
     1548      call MPI_Irecv(ciwc(:,:,:,mind),d3s1,mp_pp,id_read,MPI_ANY_TAG,&
     1549           &MPI_COMM_WORLD,reqs(j),mp_ierr)   
     1550      if (mp_ierr /= 0) goto 600
     1551    end if
    14731552
    14741553
    … …  
    14771556    goto 601
    14781557
    1479 600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
     1558600 write(*,*) "#### mpi_mod::mpif_gf_recv_vars_async> MPI ERROR ", mp_ierr
    14801559    stop
    14811560
    … …  
    14941573!
    14951574!*******************************************************************************
     1575    use com_mod, only: readclouds
     1576
    14961577    implicit none
    14971578
    1498 !    if (mp_measure_time) call mpif_mtime('commtime',0)
     1579
     1580    integer :: n_req
     1581    integer :: i
     1582
     1583!***********************************************************************
     1584
     1585    n_req=nvar_async*mp_np
     1586
     1587    if (mp_measure_time) call mpif_mtime('commtime',0)
    14991588
    15001589!    call MPI_Wait(rm1,MPI_STATUS_IGNORE,mp_ierr)
    1501     call MPI_Waitall(nvar_async*mp_np,reqs,MPI_STATUSES_IGNORE,mp_ierr)
    1502     if (mp_ierr /= 0) goto 600
    1503 
    1504 !    if (mp_measure_time) call mpif_mtime('commtime',1)
     1590
     1591! TODO: cancel recv request if readclouds=.false.
     1592!    if (readclouds) then
     1593    call MPI_Waitall(n_req,reqs,MPI_STATUSES_IGNORE,mp_ierr)
     1594!    endif
     1595    ! else
     1596    !   do i = 0, nvar_async*mp_np-1
     1597    !     if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
     1598    !       call MPI_Cancel(reqs(i),mp_ierr)
     1599    !       cycle
     1600    !     end if
     1601    !     call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
     1602    !   end do
     1603    ! end if
     1604
     1605    if (mp_ierr /= 0) goto 600
     1606
     1607    if (mp_measure_time) call mpif_mtime('commtime',1)
    15051608
    15061609    goto 601
    15071610
    1508 600 write(*,*) "#### mpi_mod::mpif_gf_request> ERROR, mp_ierr \= 0 ", mp_ierr
     1611600 write(*,*) "#### mpi_mod::mpif_gf_request> MPI ERROR ", mp_ierr
    15091612    stop
    15101613
    … …  
    16061709
    16071710!**********************************************************************
     1711
    16081712    grid_size3d=numxgridn*numygridn*numzgrid*maxspec* &
    16091713         & maxpointspec_act*nclassunc*maxageclass
    … …  
    16661770    character(LEN=*), intent(in) :: ident
    16671771    integer, intent(in) :: imode
     1772
     1773!***********************************************************************
    16681774
    16691775    select case(ident)
    … …  
    18001906
    18011907    integer :: ip,j,r
     1908
     1909!***********************************************************************
    18021910
    18031911    if (mp_measure_time) then
    … …  
    18571965! In the implementation with 3 fields, the processes may have posted
    18581966! MPI_Irecv requests that should be cancelled here
    1859 !! :TODO:
     1967!! TODO:
    18601968    ! if (.not.lmp_sync) then
    18611969    !   r=mp_pid*nvar_async
    … …  
    18681976    call MPI_FINALIZE(mp_ierr)
    18691977    if (mp_ierr /= 0) then
    1870       write(*,*) '#### mpif_finalize::MPI_FINALIZE> ERROR ####'
     1978      write(*,*) '#### mpif_finalize::MPI_FINALIZE> MPI ERROR ', mp_ierr, ' ####'
    18711979      stop
    18721980    end if
    … …  
    18871995    integer, save :: free_lun=100
    18881996    logical :: exists, iopen
     1997
     1998!***********************************************************************
    18891999
    18902000    loop1: do
    … …  
    19122022    character(LEN=40) :: fn_1, fn_2
    19132023
     2024!***********************************************************************
     2025
    19142026    write(c_ts, FMT='(I8.8,BZ)') tstep
    19152027    fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)

  • src/netcdf_output_mod.f90

    • Property mode changed from 100755 to 100644

  • src/par_mod.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    216216  ! ---------
    217217  integer,parameter :: maxwf=50000, maxtable=1000, numclass=13, ni=11
    218   !integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
    219   integer,parameter :: numwfmem=3 ! MPI with 3 fields
     218  integer,parameter :: numwfmem=2 ! Serial version/MPI with 2 fields
     219  !integer,parameter :: numwfmem=3 ! MPI with 3 fields
    220220
    221221  ! maxwf                   maximum number of wind fields to be used for simulation

  • src/partdep.f90

    • Property mode changed from 100755 to 100644

  • src/partoutput_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/partoutput_short_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/pathnames

    r0aded10 r5f9d14a  
    22../output/
    33/
    4 /flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
     4/xnilu_wrk/flex_wrk/WIND_FIELDS/AVAILABLE_ECMWF_OPER_fields_global
    55============================================
    66

  • src/photo_O1D.f90

    • Property mode changed from 100755 to 100644

  • src/random_mod.f90

    • Property mode changed from 100755 to 100644

  • src/re_initialize_particle.f90

    • Property mode changed from 100755 to 100644

  • src/readageclasses.f90

    • Property mode changed from 100755 to 100644

  • src/readavailable.f90

    • Property mode changed from 100755 to 100644

  • src/readdepo.f90

    • Property mode changed from 100755 to 100644

  • src/readpartpositions_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/readpaths.f90

    • Property mode changed from 100755 to 100644

  • src/readreleases.f90

    • Property mode changed from 100755 to 100644

  • src/readspecies.f90

    r8a65cb0 r5f9d14a  
    216216! Check scavenging parameters given in SPECIES file
    217217  if (dquer(pos_spec).gt.0) then !is particle
    218     write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
    219          &(Rain collection efficiency)  ', weta(pos_spec)
    220     write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
    221          &(Snow collection efficiency)  ', wetb(pos_spec)
     218    if (lroot) then
     219      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter A &
     220           &(Rain collection efficiency)  ', weta(pos_spec)
     221      write(*,'(a,f5.2)') '  Particle below-cloud scavenging parameter B &
     222           &(Snow collection efficiency)  ', wetb(pos_spec)
     223    end if
    222224    if (weta(pos_spec).gt.1.0 .or. wetb(pos_spec).gt.1.0) then
    223       write(*,*) '*******************************************'
    224       write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
    225            &is out of likely range'
    226       write(*,*) '          Likely range is 0.0-1.0'
    227       write(*,*) '*******************************************'
     225      if (lroot) then
     226        write(*,*) '*******************************************'
     227        write(*,*) ' WARNING: Particle below-cloud scavenging parameter A or B &
     228             &is out of likely range'
     229        write(*,*) '          Likely range is 0.0-1.0'
     230        write(*,*) '*******************************************'
     231      end if
    228232    endif
    229     write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
    230          &(CCN efficiency)  ', weta_in(pos_spec)
    231     write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
    232          &(IN efficiency)  ', wetb_in(pos_spec)
    233     if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7 )then
    234       write(*,*) '*******************************************'
    235       write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
    236       write(*,*) '          The default value is 0.9 for CCN '
    237       write(*,*) '*******************************************'
     233    if (lroot) then
     234      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Ai &
     235           &(CCN efficiency)  ', weta_in(pos_spec)
     236      write(*,'(a,f5.2)') '  Particle in-cloud scavenging parameter Bi &
     237           &(IN efficiency)  ', wetb_in(pos_spec)
     238    end if
     239    if (weta_in(pos_spec).gt.1.0 .or. weta_in(pos_spec).lt.0.7) then
     240      if (lroot) then
     241        write(*,*) '*******************************************'
     242        write(*,*) ' WARNING: Particle in-cloud scavenging parameter A is out of likely range'
     243        write(*,*) '          The default value is 0.9 for CCN '
     244        write(*,*) '*******************************************'
     245      end if
    238246    endif
    239247    if (wetb_in(pos_spec).gt.0.2 .or. wetb_in(pos_spec).lt.0.01) then
    240       write(*,*) '*******************************************'
    241       write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
    242       write(*,*) '          The default value is 0.1 for IN '
    243       write(*,*) '*******************************************'
     248      if (lroot) then
     249        write(*,*) '*******************************************'
     250        write(*,*) ' WARNING: Particle in-cloud scavenging parameter B is out of likely range'
     251        write(*,*) '          The default value is 0.1 for IN '
     252        write(*,*) '*******************************************'
     253      end if
    244254    endif
    245255
    246256  else !is gas
    247     write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
    248     write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
    249     write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
     257    if (lroot) then
     258      write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter A  ', weta(pos_spec)
     259      write(*,'(a,f5.2)') '  Gas below-cloud scavenging parameter B  ', wetb(pos_spec)
     260      write(*,*) ' Gas in-cloud scavenging uses default values as in Hertel et al 1995'
     261    end if
     262    if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
     263      if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
     264        if (lroot) then
     265          write(*,*) '*******************************************'
     266          write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
     267          write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
     268          write(*,*) '*******************************************'
     269        end if
     270      endif
     271    end if
    250272    if (wetb(pos_spec).gt.0.) then !if wet deposition is turned on
    251273      if (wetb(pos_spec).gt.0.8 .or. wetb(pos_spec).lt.0.6) then
    252         write(*,*) '*******************************************'
    253         write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
    254         write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
    255         write(*,*) '*******************************************'
    256       endif
    257     end if
    258     if (weta(pos_spec).gt.0.) then !if wet deposition is turned on
    259       if (weta(pos_spec).gt.1E-04 .or. weta(pos_spec).lt.1E-09) then
    260         write(*,*) '*******************************************'
    261         write(*,*) ' WARNING: Gas below-cloud scavenging parameter A is out of likely range'
    262         write(*,*) '          Likely range is 1E-04 to 1E-08 (see Hertel et al 1995)'
    263         write(*,*) '*******************************************'
     274        if (lroot) then
     275          write(*,*) '*******************************************'
     276          write(*,*) ' WARNING: Gas below-cloud scavenging parameter B is out of likely range'
     277          write(*,*) '          Likely range is 0.6 to 0.8 (see Hertel et al 1995)'
     278          write(*,*) '*******************************************'
     279        end if
    264280      endif
    265281    end if

  • src/readwind.f90

    r8a65cb0 r5f9d14a  
    7777!HSO  end
    7878
    79   real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    80   real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    81   real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
     79  real(sp) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     80  real(sp) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
     81  real(sp) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8282  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax
    8383
    … …  
    9191
    9292  integer :: isec1(56),isec2(22+nxmax+nymax)
    93   real(kind=4) :: zsec4(jpunp)
    94   real(kind=4) :: xaux,yaux,xaux0,yaux0
    95   real(kind=8) :: xauxin,yauxin
    96   real,parameter :: eps=1.e-4
    97   real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
    98   real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    99 
    100   logical :: hflswitch,strswitch, readcloud
     93  real(sp) :: zsec4(jpunp)
     94  real(sp) :: xaux,yaux
     95  real(dp) :: xauxin,yauxin
     96  real(sp),parameter :: eps=1.e-4
     97  real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
     98  real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
     99
     100  logical :: hflswitch,strswitch!,readcloud
    101101
    102102!HSO  grib api error messages
    … …  
    107107  strswitch=.false.
    108108!hg
    109   readcloud=.false.
     109!  readcloud=.false.
    110110!hg end
    111111  levdiff2=nlev_ec-nwz+1
    … …  
    362362        !write(*,*) 'found water!'
    363363      endif
    364 
    365364      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
    366365        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    … …  
    379378
    38037950 call grib_close_file(ifile)
    381 
    382380
    383381!error message if no fields found with correct first longitude in it

  • src/readwind_emos.f90

    • Property mode changed from 100755 to 100644

  • src/readwind_gfs_emos.f90

    • Property mode changed from 100755 to 100644

  • src/readwind_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    4141!  eso 2014:
    4242!    This version for reading windfields in advance by a separate
    43 !    MPI process.
     43!    MPI process.
     44!    TODO: can be merged with readwind.f90 if using numwfmem in
     45!          shift_field
    4446!**********************************************************************
    4547!                                                                     *
    … …  
    9799  integer :: isec1(56),isec2(22+nxmax+nymax)
    98100  real(kind=4) :: zsec4(jpunp)
    99   real(kind=4) :: xaux,yaux,xaux0,yaux0
     101  real(kind=4) :: xaux,yaux
    100102  real(kind=8) :: xauxin,yauxin
    101103  real,parameter :: eps=1.e-4
    … …  
    103105  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
    104106
    105   logical :: hflswitch,strswitch, readcloud
     107  logical :: hflswitch,strswitch !,readcloud
    106108
    107109!HSO  grib api error messages
    … …  
    115117  strswitch=.false.
    116118!hg
    117   readcloud=.false.
     119!  readcloud=.false.
    118120!hg end
    119121  levdiff2=nlev_ec-nwz+1
    … …  
    234236    elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
    235237      isec1(6)=176         ! indicatorOfParameter
    236     elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS
     238!    elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
     239    elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
    237240      isec1(6)=180         ! indicatorOfParameter
    238     elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS
     241!    elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
     242    elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
    239243      isec1(6)=181         ! indicatorOfParameter
    240244    elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO

  • src/readwind_nests_emos.f90

    • Property mode changed from 100755 to 100644

  • src/redist_mpi.f90

    • Property mode changed from 100755 to 100644

  • src/releaseparticles_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    6969  real(kind=dp) :: juldate,julmonday,jul,jullocal,juldiff
    7070  real,parameter :: eps=nxmax/3.e5,eps2=1.e-6
    71  
     71  integer :: mind2
     72! mind2        eso: pointer to 2nd windfield in memory
    7273
    7374  integer :: idummy = -7
    … …  
    291292              do kz=1,nz
    292293                if (ngrid.gt.0) then
    293                   r=p1*rhon(ix ,jy ,kz,2,ngrid) &
    294                        +p2*rhon(ixp,jy ,kz,2,ngrid) &
    295                        +p3*rhon(ix ,jyp,kz,2,ngrid) &
    296                        +p4*rhon(ixp,jyp,kz,2,ngrid)
    297                   t=p1*ttn(ix ,jy ,kz,2,ngrid) &
    298                        +p2*ttn(ixp,jy ,kz,2,ngrid) &
    299                        +p3*ttn(ix ,jyp,kz,2,ngrid) &
    300                        +p4*ttn(ixp,jyp,kz,2,ngrid)
     294                  r=p1*rhon(ix ,jy ,kz,mind2,ngrid) &
     295                       +p2*rhon(ixp,jy ,kz,mind2,ngrid) &
     296                       +p3*rhon(ix ,jyp,kz,mind2,ngrid) &
     297                       +p4*rhon(ixp,jyp,kz,mind2,ngrid)
     298                  t=p1*ttn(ix ,jy ,kz,mind2,ngrid) &
     299                       +p2*ttn(ixp,jy ,kz,mind2,ngrid) &
     300                       +p3*ttn(ix ,jyp,kz,mind2,ngrid) &
     301                       +p4*ttn(ixp,jyp,kz,mind2,ngrid)
    301302                else
    302                   r=p1*rho(ix ,jy ,kz,2) &
    303                        +p2*rho(ixp,jy ,kz,2) &
    304                        +p3*rho(ix ,jyp,kz,2) &
    305                        +p4*rho(ixp,jyp,kz,2)
    306                   t=p1*tt(ix ,jy ,kz,2) &
    307                        +p2*tt(ixp,jy ,kz,2) &
    308                        +p3*tt(ix ,jyp,kz,2) &
    309                        +p4*tt(ixp,jyp,kz,2)
     303                  r=p1*rho(ix ,jy ,kz,mind2) &
     304                       +p2*rho(ixp,jy ,kz,mind2) &
     305                       +p3*rho(ix ,jyp,kz,mind2) &
     306                       +p4*rho(ixp,jyp,kz,mind2)
     307                  t=p1*tt(ix ,jy ,kz,mind2) &
     308                       +p2*tt(ixp,jy ,kz,mind2) &
     309                       +p3*tt(ix ,jyp,kz,mind2) &
     310                       +p4*tt(ixp,jyp,kz,mind2)
    310311                endif
    311312                press=r*r_air*t/100.
    … …  
    377378              if (ngrid.gt.0) then
    378379                do n=1,2
    379                   rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,2,ngrid) &
    380                        +p2*rhon(ixp,jy ,indz+n-1,2,ngrid) &
    381                        +p3*rhon(ix ,jyp,indz+n-1,2,ngrid) &
    382                        +p4*rhon(ixp,jyp,indz+n-1,2,ngrid)
     380                  rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,mind2,ngrid) &
     381                       +p2*rhon(ixp,jy ,indz+n-1,mind2,ngrid) &
     382                       +p3*rhon(ix ,jyp,indz+n-1,mind2,ngrid) &
     383                       +p4*rhon(ixp,jyp,indz+n-1,mind2,ngrid)
    383384                end do
    384385              else
    385386                do n=1,2
    386                   rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,2) &
    387                        +p2*rho(ixp,jy ,indz+n-1,2) &
    388                        +p3*rho(ix ,jyp,indz+n-1,2) &
    389                        +p4*rho(ixp,jyp,indz+n-1,2)
     387                  rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,mind2) &
     388                       +p2*rho(ixp,jy ,indz+n-1,mind2) &
     389                       +p3*rho(ix ,jyp,indz+n-1,mind2) &
     390                       +p4*rho(ixp,jyp,indz+n-1,mind2)
    390391                end do
    391392              endif

  • src/timemanager_mpi.f90

    • Property mode changed from 100755 to 100644
    r8a65cb0 r5f9d14a  
    180180    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) then
    181181! readwind process skips this step
    182       if (.not.(lmpreader.and.lmp_use_reader)) then !.or..not.lmp_use_reader)
     182      if (.not.(lmpreader.and.lmp_use_reader)) then
    183183        call ohreaction(itime,lsynctime,loutnext)
    184184      endif
    … …  
    223223      call mpif_gf_send_vars(memstat)
    224224      call mpif_gf_send_vars_nest(memstat)
    225 ! This version is also used whenever 2 fields are needed (in this case,
    226 ! async send/recv is impossible)
     225! Version 2  (lmp_sync=.false.) is also used whenever 2 new fields are read (in which
     226! case async send/recv is impossible.
    227227    else if (.not.lmp_sync.and.lmp_use_reader.and.memstat.ge.32) then
    228228      call mpif_gf_send_vars(memstat)
    … …  
    231231
    232232! Version 2 (lmp_sync=.false.) is for holding three fields in memory. Uses a
    233 ! read-ahead process where sending/receiving of the 3rd fields are done in
     233! read-ahead process where sending/receiving of the 3rd fields is done in
    234234! the background in parallel with performing computations with fields 1&2
    235235!********************************************************************************
    … …  
    241241      end if
    242242
    243 
    244243! COMPLETION CHECK:
    245244! Issued at start of each new field period.
    246245      if (memstat.ne.0.and.memstat.lt.32.and.lmp_use_reader) then
    247 ! :DEV: z0(7) changes with time, so should be dimension (numclass,2)
    248 ! to allow transfer of the future value in the background
     246! TODO: z0(7) changes with time, so should be dimension (numclass,2) to
     247! allow transfer of the future value in the background
     248        call MPI_Bcast(z0,numclass,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
    249249        call mpif_gf_request
    250         call MPI_Bcast(z0,numclass,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
    251250      end if
    252251
    253 
    254252! RECVEIVING PROCESS(ES):
     253      ! eso TODO: at this point we do not know if clwc/ciwc will be available
     254      ! at next time step. Issue receive request anyway, cancel at mpif_gf_request
    255255      if (memstat.gt.0.and.lmp_use_reader.and..not.lmpreader) then
    256256        call mpif_gf_recv_vars_async(memstat)
    … …  
    755755    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
    756756
     757! TODO: delete for release version?
    757758!!-------------------------------------------------------------------------------
    758759! These lines below to test the well-mixed condition, modified by  mc, not to
    … …  
    778779
    779780! eso :TODO: this not implemented yet (transfer particles to PID 0 or rewrite)
    780 ! the tools to do this is in mpi_mod.f90
     781! the tools to do this are already in mpi_mod.f90
    781782  if (lroot) then
    782783    do j=1,numpart
    … …  
    787788
    788789  if (ipout.eq.2) then
    789 ! MPI: process 0 creates the file, the other processes append to it
     790! MPI process 0 creates the file, the other processes append to it
    790791    do ip=0, mp_partgroup_np-1
    791792      if (ip.eq.mp_partid) then
    792         if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
     793        !if (mp_dbg_mode) write(*,*) 'call partoutput(itime), proc, mp_partid',ip,mp_partid
    793794        call partoutput(itime)    ! dump particle positions
    794795      end if

  • src/wetdepo.f90

    r8a65cb0 r5f9d14a  
    7979  real :: wetdeposit(maxspec),restmass
    8080  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    81   save lfr,cfr
    82 
    83   real :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
    84   real :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
     81  !save lfr,cfr
     82
     83  real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
     84  real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
    8585
    8686!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
    87   real :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
    88   real :: bcls(6) = (/ 22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
     87  real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
     88  real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
    8989  real :: frac_act, liq_frac, dquer_m
    9090
    … …  
    159159  !********************************************************************
    160160    interp_time=nint(itime-0.5*ltsample)
    161 
     161   
    162162    if (ngrid.eq.0) then
    163163      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
    … …  
    171171
    172172!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
    173      if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
     173    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
    174174
    175175  ! get the level were the actual particle is in
    … …  
    231231
    232232
    233 !ZHG calculated for 1) both 2) lsp 3) convp
     233  !ZHG calculated for 1) both 2) lsp 3) convp
    234234    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
    235235    grfraction(2)=max(0.05,cc*(lfr(i)))
    … …  
    252252      wetscav=0.   
    253253
    254      if (ngrid.gt.0) then
    255        act_temp=ttn(ix,jy,hz,n,ngrid)
    256      else
    257        act_temp=tt(ix,jy,hz,n)
    258      endif
     254      if (ngrid.gt.0) then
     255        act_temp=ttn(ix,jy,hz,n,ngrid)
     256      else
     257        act_temp=tt(ix,jy,hz,n)
     258      endif
    259259     
    260260
    … …  
    264264      if (clouds_v.ge.4) then !below cloud
    265265
    266          if (weta(ks).gt.0. .or. wetb(ks).gt.0) then !if positive below-cloud parameters given in SPECIES file (either A or B)
    267             blc_count=blc_count+1
     266        if (weta(ks).gt.0. .or. wetb(ks).gt.0) then !if positive below-cloud parameters given in SPECIES file (either A or B)
     267          blc_count=blc_count+1
    268268
    269269!GAS
    270            if (dquer(ks) .le. 0.) then  !is gas
    271             !Gas scavenging coefficient based on Hertel et al 1995 using the below-cloud scavenging parameters A (=weta) and B (=wetb) from SPECIES file
     270          if (dquer(ks) .le. 0.) then  !is gas
     271  ! Gas scavenging coefficient based on Hertel et al 1995 using the below-cloud scavenging parameters A (=weta) and B (=wetb) from SPECIES file
    272272            wetscav=weta(ks)*prec(1)**wetb(ks)
    273273
    274 !AEROSOL
    275            else !is particle
    276              !NIK 17.02.2015
    277              ! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
    278              dquer_m=dquer(ks)/1000000. !conversion from um to m
    279 
    280              !ZHG snow or rain removal is applied based on the temperature.
    281              if (act_temp .ge. 273)  then !Rain
    282  
    283              !Particle RAIN scavenging coefficient based on Laakso et al 2003, the below-cloud scavenging (rain efficienty) parameter A (=weta) from SPECIES file
    284                  wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
    285                 (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
    286 
    287              elseif (act_temp .lt. 273)  then !snow
    288  
    289              !Particle SNOW scavenging coefficient based on Kyro et al 2009, the below-cloud scavenging (Snow efficiency) parameter B (=wetb) from SPECIES file
    290                 wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
    291                 (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
    292 
    293              endif         
     274  !AEROSOL
     275          else !is particle
     276  !NIK 17.02.2015
     277  ! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
     278            dquer_m=dquer(ks)/1000000. !conversion from um to m
     279           
     280  !ZHG snow or rain removal is applied based on the temperature.
     281            if (act_temp .ge. 273)  then !Rain
     282
     283  !Particle RAIN scavenging coefficient based on Laakso et al 2003, the below-cloud scavenging (rain efficienty) parameter A (=weta) from SPECIES file
     284              wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
     285                   (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
     286
     287            elseif (act_temp .lt. 273)  then !snow
     288
     289  !Particle SNOW scavenging coefficient based on Kyro et al 2009, the below-cloud scavenging (Snow efficiency) parameter B (=wetb) from SPECIES file
     290              wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
     291                   (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
     292
     293            endif
    294294
    295295!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart
    296296
    297            endif !gas or particle
    298          endif ! positive below-cloud scavenging parameters given in Species file
     297          endif !gas or particle
     298        endif ! positive below-cloud scavenging parameters given in Species file
    299299      endif !end below-cloud
    300300
    … …  
    305305      if (clouds_v.lt.4) then !in-cloud
    306306
    307           inc_count=inc_count+1
    308 
    309           !ZHG liquid water parameterization (CLWC+CIWC)
    310           if (readclouds) then !get cloud water clwc & ciwc units Kg/Kg
    311              cl=clwc(ix,jy,hz,n)+ciwc(ix,jy,hz,n)
    312           else !parameterize cloudwater
    313              cl=2E-7*prec(1)**0.36
     307        inc_count=inc_count+1
     308
     309  !ZHG liquid water parameterization (CLWC+CIWC)
     310        if (readclouds) then !get cloud water clwc & ciwc units Kg/Kg
     311          cl=clwc(ix,jy,hz,n)+ciwc(ix,jy,hz,n)
     312        else !parameterize cloudwater
     313          cl=2E-7*prec(1)**0.36
     314        endif
     315
     316  ! AEROSOL
     317        if (dquer(ks).gt. 0.) then ! is particle
     318          if (act_temp .le. 253) then
     319            liq_frac=0
     320          else if (act_temp .ge. 273) then
     321            liq_frac=1
     322          else
     323            liq_frac =((act_temp-273)/(273-253))**2
    314324          endif
    315325
    316 ! AEROSOL
    317           if (dquer(ks).gt. 0.) then ! is particle
    318             if (act_temp .le. 253) then
    319               liq_frac=0
    320             else if (act_temp .ge. 273) then
    321               liq_frac=1
    322             else
    323               liq_frac =((act_temp-273)/(273-253))**2
    324             endif
    325 
    326             !ZHG  calculate the activated fraction based on the In-cloud scavenging parameters Ai (=weta_in) and Bi (=wetb_in) from SPECIES file
    327             ! frac_act is the combined IN and CCN efficiency
    328             ! The default values are 0.9 for CCN and 0.1 IN
    329             ! This parameterization is based on Verheggen et al. (2007) & Cozich et al. (2006)
    330             frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
     326! ZHG  calculate the activated fraction based on the In-cloud scavenging parameters Ai (=weta_in) and Bi (=wetb_in) from SPECIES file
     327! frac_act is the combined IN and CCN efficiency
     328! The default values are 0.9 for CCN and 0.1 IN
     329! This parameterization is based on Verheggen et al. (2007) & Cozich et al. (2006)
     330          frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
    331331 
    332             !ZHG Use the activated fraction and the liqid water to calculate the washout
    333             S_i= frac_act/cl
    334 
    335 ! GAS
    336           else ! is gas
    337             cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
    338             S_i=1/cle
    339           endif
    340 
    341           ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
    342           wetscav=S_i*(prec(1)/3.6E6)/clouds_h
     332  !ZHG Use the activated fraction and the liqid water to calculate the washout
     333          S_i= frac_act/cl
     334
     335  ! GAS
     336        else ! is gas
     337          cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
     338          S_i=1/cle
     339        endif
     340
     341  ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
     342        wetscav=S_i*(prec(1)/3.6E6)/clouds_h
    343343
    344344!          write(*,*) 'in-cloud, act_temp=',act_temp,',prec=',prec(1),',wetscav=',wetscav,',jpart=',jpart,',clouds_h=,', &
    … …  
    346346
    347347      endif !incloud
    348 
     348     
    349349  !**************************************************
    350350  ! CALCULATE DEPOSITION
    … …  
    355355      if (wetscav.gt.0.) then
    356356        wetdeposit(ks)=xmass1(jpart,ks)* &
    357         (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
     357             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
    358358      else ! if no scavenging
    359359        wetdeposit(ks)=0.
    … …  
    389389    if (ldirect.eq.1) then
    390390      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
    391       real(ytra1(jpart)),nage,kp)
     391           real(ytra1(jpart)),nage,kp)
    392392      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
    393         wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
     393           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
    394394    endif
    395395
    … …  
    397397  end do ! all particles
    398398
    399 ! count the total number of below-cloud and in-cloud occurences:
    400    tot_blc_count=tot_blc_count+blc_count
    401    tot_inc_count=tot_inc_count+inc_count
     399  ! count the total number of below-cloud and in-cloud occurences:
     400  tot_blc_count=tot_blc_count+blc_count
     401  tot_inc_count=tot_inc_count+inc_count
    402402
    403403end subroutine wetdepo
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG