Changeset 41d8574 in flexpart.git


Ignore:
Timestamp:
Jan 15, 2016, 12:35:58 PM (8 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:
fdc0f03
Parents:
d8107c2
Message:

bugfix: MPI version gave wrong wet deposition when using ECMWF cloud water fields. Cloud water in ECMWF fields now uses parameter qc, or reads clwc+ciwc. Added minmass variable as limit for terminating particles.

Location:
src
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • src/com_mod.f90

    rf75967d r41d8574  
    351351  real :: qv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
    352352!ZHG adding cloud water
    353   real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !liquid   [kg/kg]
    354   real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem) !ice      [kg/kg]
    355   real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)  !combined [m3/m3]
     353  real :: clwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !liquid   [kg/kg]
     354  real :: ciwc(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0 !ice      [kg/kg]
     355  real :: clw(0:nxmax-1,0:nymax-1,nzmax,numwfmem)=0.0  !combined [m3/m3]
    356356
    357357  real :: pv(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
     
    360360  real :: tth(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
    361361  real :: qvh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
    362   real :: clwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
    363   real :: ciwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
     362  real :: clwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)=0.0
     363  real :: ciwch(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)=0.0
    364364
    365365  real :: pplev(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
     
    371371   real :: icloud_stats(0:nxmax-1,0:nymax-1,5,numwfmem)
    372372   real :: clh(0:nxmax-1,0:nymax-1,nuvzmax,numwfmem)
     373   real :: clw4(0:nxmax-1,0:nymax-1,numwfmem) ! eso: =icloud_stats(:,:,4,:)
    373374
    374375
  • src/mpi_mod.f90

    rd6a0977 r41d8574  
    119119  logical, parameter :: mp_dbg_out = .false.
    120120  logical, parameter :: mp_time_barrier=.true.
    121   logical, parameter :: mp_measure_time=.false.
     121  logical, parameter :: mp_measure_time=.true.
    122122  logical, parameter :: mp_exact_numpart=.true.
    123123
     
    907907!**********************************************************************
    908908
    909 ! The non-reader processes need to know if clouds were read.
     909! The non-reader processes need to know if cloud water were read.
    910910! TODO: only at first step or always?
    911911    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
     
    963963! cloud water/ice:
    964964    if (readclouds) then
    965       call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2_size1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
    966       if (mp_ierr /= 0) goto 600
    967 
     965      ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
     966      ! if (mp_ierr /= 0) goto 600
     967      call MPI_Bcast(clw4(:,:,li:ui),d2s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
     968      if (mp_ierr /= 0) goto 600
    968969      ! call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_pp,id_read,MPI_COMM_WORLD,mp_ierr)
    969970      ! if (mp_ierr /= 0) goto 600
     
    13361337      if (readclouds) then
    13371338        i=i+1
    1338         call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,dest,tm1,&
     1339        ! call MPI_Isend(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,dest,tm1,&
     1340        !      &MPI_COMM_WORLD,reqs(i),mp_ierr)
     1341        call MPI_Isend(clw4(:,:,mind),d2s1,mp_pp,dest,tm1,&
    13391342             &MPI_COMM_WORLD,reqs(i),mp_ierr)
     1343
    13401344        if (mp_ierr /= 0) goto 600
    13411345
     
    15431547      j=j+1
    15441548
    1545       call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
     1549      ! call MPI_Irecv(icloud_stats(:,:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
     1550      !      &MPI_COMM_WORLD,reqs(j),mp_ierr)
     1551      call MPI_Irecv(clw4(:,:,mind),d2s1*5,mp_pp,id_read,MPI_ANY_TAG,&
    15461552           &MPI_COMM_WORLD,reqs(j),mp_ierr)
    15471553      if (mp_ierr /= 0) goto 600
  • src/par_mod.f90

    rd6a0977 r41d8574  
    211211
    212212  integer,parameter :: maxpart=40000000
    213 !  integer,parameter :: maxpart=60000000
    214 !  integer,parameter :: maxpart=120000000
    215213  integer,parameter :: maxspec=6
     214  integer,parameter :: minmass=0.0 !0.0001
    216215
    217216  ! maxpart                 Maximum number of particles
    218217  ! maxspec                 Maximum number of chemical species per release
    219 
     218  ! minmass                 Terminate particles carrying less mass
    220219
    221220  ! maxpoint is also set dynamically during runtime
  • src/readwind.f90

    rd6a0977 r41d8574  
    190190    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    191191      isec1(6)=133         ! indicatorOfParameter
    192 !ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
     192! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc
    193193    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
    194194      isec1(6)=246         ! indicatorOfParameter
    195 ! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
    196 !    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
    197 !      isec1(6)=247         ! indicatorOfParameter
    198 !ZHG end
     195    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     196      isec1(6)=247         ! indicatorOfParameter
     197! ESO qc(=clwc+ciwc):
     198    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
     199      isec1(6)=201031         ! indicatorOfParameter
    199200    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    200201      isec1(6)=134         ! indicatorOfParameter
     
    360361      if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
    361362!ZHG READING CLOUD FIELDS ASWELL
     363! ESO TODO: add check for if one of clwc/ciwc missing (error),
     364! also if all 3 cw fields present, use qc and disregard the others
    362365      if(isec1(6).eq.246) then  !! CLWC  Cloud liquid water content [kg/kg]
    363366        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    364         readclouds = .true.
     367        readclouds=.true.
     368        sumclouds=.false.
     369      endif
     370      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
     371        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     372      endif
     373!ZHG end
     374!ESO read qc (=clwc+ciwc)
     375      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
     376        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     377        readclouds=.true.
     378        sumclouds=.true.
    365379!if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
    366380      endif
    367 !      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
    368 !        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    369         !write(*,*) 'found ice!'
    370 !      endif
    371 !ZHG end
    372381
    373382    end do
     
    427436!ZHG
    428437    call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,2,n)
    429     call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
     438    if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,2,n)
    430439!ZHG end
    431440
  • src/readwind_mpi.f90

    rd6a0977 r41d8574  
    8383!HSO  end
    8484
    85   real(kind=4), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
    86   real(kind=4), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
    87   real(kind=4), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
     85  real(sp), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
     86  real(sp), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
     87  real(sp), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
    8888  integer, intent(in) :: indj,n
    8989  integer :: i,j,k,levdiff2,ifield,iumax,iwmax
     
    116116  hflswitch=.false.
    117117  strswitch=.false.
    118 !hg
     118!ZHG test the grib fields that have lcwc without using them
    119119!  readcloud=.false.
    120120!hg end
     
    200200    elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
    201201      isec1(6)=133         ! indicatorOfParameter
    202 !ZHG ! READ CLOUD FIELD : assume these to be added together to one variable
     202! ESO Cloud water is in a) fields CLWC and CIWC, *or* b) field qc
    203203    elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
    204204      isec1(6)=246         ! indicatorOfParameter
    205 ! ICE AND WATER IS ADDED TOGETHER IN NEW WINDFIELDS
    206 !    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
    207 !      isec1(6)=247         ! indicatorOfParameter
    208 !ZHG end
     205    elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
     206      isec1(6)=247         ! indicatorOfParameter
     207! ESO qc(=clwc+ciwc):
     208    elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
     209      isec1(6)=201031         ! indicatorOfParameter
    209210    elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
    210211      isec1(6)=134         ! indicatorOfParameter
     
    373374        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    374375        readclouds = .true.
    375 !if (clwch(i,j,nlev_ec-k+2,n) .gt. 0)        write(*,*) 'readwind: found water!', clwch(i,j,nlev_ec-k+2,n)
    376       endif
    377 !      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
    378 !        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
    379         !write(*,*) 'found ice!'
    380 !      endif
     376        sumclouds=.false.
     377      endif
     378      if(isec1(6).eq.247) then  !! CIWC  Cloud ice water content
     379        ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     380      endif
    381381!ZHG end
     382!ESO read qc (=clwc+ciwc)
     383      if(isec1(6).eq.201031) then  !! QC  Cloud liquid water content [kg/kg]
     384        clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
     385        readclouds=.true.
     386        sumclouds=.true.
     387      endif
    382388
    383389    end do
  • src/timemanager.f90

    rf75967d r41d8574  
    607607          end do
    608608
    609           if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
     609          if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass
    610610            itra1(j)=-999999999
    611611            if (verbosity.gt.0) then
     
    644644  !***************************************************************
    645645   
    646           total_nan_intl=0
    647           i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
    648           sum_nan_count(i_nan)=nan_count
    649           if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
    650           do ii_nan=1, (maxtl/lsynctime)
    651               total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
    652           end do
     646    total_nan_intl=0
     647    i_nan=i_nan+1 ! added by mc to count nan during a time of maxtl (i.e. maximum tl fixed here to 20 minutes, see com_mod)
     648    sum_nan_count(i_nan)=nan_count
     649    if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
     650    do ii_nan=1, (maxtl/lsynctime)
     651      total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
     652    end do
    653653  ! Output to keep track of the numerical instabilities in CBL simulation and if
    654654  ! they are compromising the final result (or not)
    655           if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
     655    if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
    656656         
    657 !!------------------------------------------------------------------------------
    658 ! these lines below to test the well-mixed condition, modified by  mc, not to be
    659 ! included in final release:
    660     ! if (itime.eq.0) then
    661     ! open(551,file='/home/mc/test_cbl/out/WELLMIXEDTEST_CBL_lonlat_9_33_100_3hours_3htp_cd.DAT')
    662     ! open(552,file='/home/mc/test_cbl/out/avg_ol_h_wst_lonlat_9_33_100_3hours_3htp_cd.DAT')
    663     ! end if
    664     ! write(552,'(5F16.7)')itime*1./3600.,avg_wst/well_mixed_norm,avg_ol/well_mixed_norm,avg_h/well_mixed_norm
    665     ! do j=1,25
    666     !      !write(551,*))itime*1.,h_well/50.*j,well_mixed_vector(j)/well_mixed_norm*50.
    667     !     avg_air_dens(j)=avg_air_dens(j)/well_mixed_vector(j)
    668     !     write(551,'(5F16.7)')itime*1.,h_well/25.*j,well_mixed_vector(j)/well_mixed_norm*25.,avg_air_dens(j),0.04*j
    669     ! end do
    670 !!------------------------------------------------------------------------------
    671 
    672657  end do
    673658
  • src/timemanager_mpi.f90

    rf75967d r41d8574  
    103103  implicit none
    104104
    105   logical :: fc_mp=.true., fc_sync=.true.
    106105  logical :: reqv_state=.false. ! .true. if waiting for a MPI_Irecv to complete
    107106  integer :: j,ks,kp,l,n,itime=0,nstop,nstop1,memstat=0,mind
     
    352351      do ks=1,nspec
    353352        do kp=1,maxpointspec_act
    354           if (decay(ks).gt.0.) then
     353          if (decay(ks).gt.0.) then ! TODO move this statement up 2 levels
    355354            do nage=1,nageclass
    356355              do l=1,nclassunc
     
    485484            if (lroot) then
    486485              if (lnetcdfout.eq.1) then
    487 
    488486                call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,&
    489487                     &drygridtotalunc)
     
    772770          end do
    773771
    774           if (xmassfract.lt.0.00005 .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then   ! terminate all particles carrying less mass
     772          if (xmassfract.lt.minmass) then ! .and. sum(real(npart(npoint(j)))*xmass1(j,:)).lt.1.0) then   ! terminate all particles carrying less mass
    775773          !            print*,'terminated particle ',j,' for small mass (', sum(real(npart(npoint(j)))* &
    776774          !         xmass1(j,:)), ' of ', sum(xmass(npoint(j),:)),')'
  • src/verttransform.f90

    rf75967d r41d8574  
    251251!hg adding the cloud water
    252252  clwc(:,:,1,n)=clwch(:,:,1,n)
    253   ciwc(:,:,1,n)=ciwch(:,:,1,n)   
     253  if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)   
    254254!hg
    255255  pv(:,:,1,n)=pvh(:,:,1)
     
    262262!hg adding the cloud water
    263263  clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
    264   ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
     264  if (.not.sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
    265265!hg
    266266  pv(:,:,nz,n)=pvh(:,:,nuvz)
     
    280280!hg adding the cloud water
    281281          clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
    282           ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
     282          if (.not.sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
    283283!hg
    284284          pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
     
    310310!hg adding the cloud water
    311311          clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz
    312           ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
     312          if (.not.sumclouds) &
     313               &ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz
    313314!hg
    314315          pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
     
    578579!***********************************************************************************
    579580    write(*,*) 'using cloud water from ECMWF'
    580     clw(:,:,:,n)=0
    581     icloud_stats(:,:,:,n)=0
     581    clw(:,:,:,n)=0.0
     582    icloud_stats(:,:,:,n)=0.0
    582583    clouds(:,:,:,n)=0
     584! If water/ice are read separately into clwc and ciwc, store sum in clwc
     585    if (.not.sumclouds) then
     586      clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
     587    end if
    583588    do jy=0,nymin1
    584589      do ix=0,nxmin1
     
    635640!   create a cloud and rainout/washout field, clouds occur where rh>80%
    636641!   total cloudheight is stored at level 0
    637 ! if (.not.readclouds) write(*,*) 'using cloud water from Parameterization'
    638642    write(*,*) 'using cloud water from Parameterization'
    639643    do jy=0,nymin1
     
    676680    end do
    677681  endif !readclouds
     682
     683! eso: copy the relevant data to clw4 to reduce amount of communicated data for MPI
     684  clw4(:,:,n) = icloud_stats(:,:,4,n)
    678685
    679686     !********* TEST ***************
  • src/wetdepo.f90

    rf75967d r41d8574  
    317317          !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
    318318          if (readclouds) then                  !icloud_stats(ix,jy,4,n) has units kg/m2
    319             cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
     319!            cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
     320! ESO: stop using icloud_stats
     321            cl =clw4(ix,jy,n)*(grfraction(1)/cc)
    320322          else                                  !parameterize cloudwater m2/m3
    321323            !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG