Changeset 8a65cb0 in flexpart.git for src/timemanager.f90


Ignore:
Timestamp:
Mar 2, 2015, 3:11:55 PM (9 years ago)
Author:
Espen Sollum ATMOS <espen@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
1d207bb
Parents:
60403cd
Message:

Added code, makefile for dev branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/timemanager.f90

    r1a08571 r8a65cb0  
    4343  !     call convection BEFORE new fields are read in BWD mode                 *
    4444  !  Changes Caroline Forster, Feb 2005                                        *
    45   !new interface between flexpart and convection scheme                        *
    46   !Emanuel's latest subroutine convect43c.f is used                            *
     45  !   new interface between flexpart and convection scheme                     *
     46  !   Emanuel's latest subroutine convect43c.f is used                         *
     47  !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
     48  !   added netcdf output code                                                 *
     49  !  Changes Espen Sollum 2014                                                 *
     50  !   For compatibility with MPI version,                                      *
     51  !   variables uap,ucp,uzp,us,vs,ws,cbt now in module com_mod                 *
    4752  !*****************************************************************************
    4853  !                                                                            *
    … …  
    9297  use par_mod
    9398  use com_mod
     99  use netcdf_output_mod, only: concoutput_netcdf,concoutput_nest_netcdf,&
     100       &concoutput_surf_netcdf,concoutput_surf_nest_netcdf
    94101
    95102  implicit none
    … …  
    99106  integer :: loutnext,loutstart,loutend
    100107  integer :: ix,jy,ldeltat,itage,nage
    101   real :: outnum,weight,prob(maxspec)
    102   real :: uap(maxpart),ucp(maxpart),uzp(maxpart),decfact
    103   real :: us(maxpart),vs(maxpart),ws(maxpart)
    104   integer(kind=2) :: cbt(maxpart)
     108  integer :: i_nan=0,ii_nan,total_nan_intl=0  !added by mc to check instability in CBL scheme
     109  real :: outnum,weight,prob(maxspec),decfact
     110  ! real :: uap(maxpart),ucp(maxpart),uzp(maxpart)
     111  ! real :: us(maxpart),vs(maxpart),ws(maxpart)
     112  ! integer(kind=2) :: cbt(maxpart)
    105113  real :: drydeposit(maxspec),gridtotalunc,wetgridtotalunc
    106114  real :: drygridtotalunc,xold,yold,zold,xmassfract
    … …  
    129137
    130138
    131   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    132   itime=0 ! initialise to avoid random numbers on output IP 2015-03-02
    133   write(*,46) float(itime)/3600,itime,numpart
    134139  if (verbosity.gt.0) then
    135140    write (*,*) 'timemanager> starting simulation'
    … …  
    141146
    142147  do itime=0,ideltas,lsynctime
    143     if (verbosity.gt.0) then
    144            write (*,*) 'timemanager>  itime=', itime
    145     endif
    146 
    147148
    148149  ! Computation of wet deposition, OH reaction and mass transfer
    … …  
    198199    if (nstop1.gt.1) stop 'NO METEO FIELDS AVAILABLE'
    199200
     201  ! Get hourly OH fields if not available
     202  !****************************************************
     203    if (OHREA) then
     204      if (verbosity.gt.0) then
     205             write (*,*) 'timemanager> call gethourlyOH'
     206      endif
     207      call gethourlyOH(itime)
     208          if (verbosity.gt.1) then
     209            CALL SYSTEM_CLOCK(count_clock)
     210            WRITE(*,*) 'timemanager> SYSTEM CLOCK',(count_clock - count_clock0)/real(count_rate)
     211          endif
     212    endif
     213       
    200214  ! Release particles
    201215  !******************
    202216
    203     !if (verbosity.gt.0) then
    204     !       write (*,*) 'timemanager>  Release particles'
    205     !endif
     217    if (verbosity.gt.0) then
     218           write (*,*) 'timemanager>  Release particles'
     219    endif
    206220
    207221    if (mdomainfill.ge.1) then
    … …  
    219233    else
    220234      if (verbosity.gt.0) then
    221         print*,'timemanager> call releaseparticles' 
     235        print*,'call releaseparticles' 
    222236      endif
    223237      call releaseparticles(itime)
    … …  
    354368        if ((iout.le.3.).or.(iout.eq.5)) then
    355369          if (surf_only.ne.1) then
    356           call concoutput(itime,outnum,gridtotalunc, &
    357                wetgridtotalunc,drygridtotalunc)
     370            if (lnetcdfout.eq.1) then
     371              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     372            else
     373              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     374            endif
    358375          else 
    359   if (verbosity.eq.1) then
    360      print*,'call concoutput_surf '
    361      CALL SYSTEM_CLOCK(count_clock)
    362      WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
    363   endif
    364           call concoutput_surf(itime,outnum,gridtotalunc, &
    365                wetgridtotalunc,drygridtotalunc)
    366   if (verbosity.eq.1) then
    367      print*,'called concoutput_surf '
    368      CALL SYSTEM_CLOCK(count_clock)
    369      WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
    370   endif
     376            if (verbosity.eq.1) then
     377             print*,'call concoutput_surf '
     378             call system_clock(count_clock)
     379             write(*,*) 'system clock',count_clock - count_clock0   
     380            endif
     381            if (lnetcdfout.eq.1) then
     382              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     383            else
     384              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
     385              if (verbosity.eq.1) then
     386                print*,'called concoutput_surf '
     387                call system_clock(count_clock)
     388                write(*,*) 'system clock',count_clock - count_clock0   
     389              endif
     390            endif
    371391          endif
    372392
    373           if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
    374           if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
     393          if (nested_output .eq. 1) then
     394            if (lnetcdfout.eq.0) then
     395              if (surf_only.ne.1) then
     396                call concoutput_nest(itime,outnum)
     397              else
     398                call concoutput_surf_nest(itime,outnum)
     399              endif
     400            else
     401              if (surf_only.ne.1) then
     402                call concoutput_nest_netcdf(itime,outnum)
     403              else
     404                call concoutput_surf_nest_netcdf(itime,outnum)
     405              endif
     406            endif
     407          endif
    375408          outnum=0.
    376409        endif
    377410        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    378411        if (iflux.eq.1) call fluxoutput(itime)
    379         !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
    380         write(*,46) float(itime)/3600,itime,numpart
     412        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
     413        !write(*,46) float(itime)/3600,itime,numpart
    38141445      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
    38241546      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
    … …  
    448481  ! Loop over all particles
    449482  !************************
    450 
     483  ! Various variables for testing reason of CBL scheme, by mc
     484    well_mixed_vector=0. !erase vector to test well mixed condition: modified by mc
     485    well_mixed_norm=0.   !erase normalization to test well mixed condition: modified by mc
     486    avg_ol=0.
     487    avg_wst=0.
     488    avg_h=0.
     489    avg_air_dens=0.  !erase vector to obtain air density at particle positions: modified by mc
     490  !-----------------------------------------------------------------------------
    451491    do j=1,numpart
    452492
    … …  
    542582          if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
    543583            itra1(j)=-999999999
     584            if (verbosity.gt.0) then
     585              print*,'terminated particle ',j,' for small mass'
     586            endif
    544587          endif
    545588
    … …  
    558601
    559602          if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
    560             if (linit_cond.ge.1) &
    561                  call initial_cond_calc(itime+lsynctime,j)
     603            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
    562604            itra1(j)=-999999999
     605            if (verbosity.gt.0) then
     606              print*,'terminated particle ',j,' for age'
     607            endif
    563608          endif
    564609        endif
    … …  
    566611      endif
    567612
    568     end do
     613    end do !loop over particles
     614   
     615  ! Counter of "unstable" particle velocity during a time scale of
     616  ! maximumtl=20 minutes (defined in com_mod)
     617  !***************************************************************
     618   
     619          total_nan_intl=0
     620          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)
     621          sum_nan_count(i_nan)=nan_count
     622          if (i_nan > maxtl/lsynctime) i_nan=1 !lsynctime must be <= maxtl
     623          do ii_nan=1, (maxtl/lsynctime)
     624              total_nan_intl=total_nan_intl+sum_nan_count(ii_nan)
     625          end do
     626  ! Output to keep track of the numerical instabilities in CBL simulation and if
     627  ! they are compromising the final result (or not)
     628          if (cblflag.eq.1) print *,j,itime,'nan_synctime',nan_count,'nan_tl',total_nan_intl 
     629         
     630!!------------------------------------------------------------------------------
     631! these lines below to test the well-mixed condition, modified by  mc, not to be
     632! included in final release:
     633    ! if (itime.eq.0) then
     634    ! open(551,file='/home/mc/test_cbl/out/WELLMIXEDTEST_CBL_lonlat_9_33_100_3hours_3htp_cd.DAT')
     635    ! open(552,file='/home/mc/test_cbl/out/avg_ol_h_wst_lonlat_9_33_100_3hours_3htp_cd.DAT')
     636    ! end if
     637    ! write(552,'(5F16.7)')itime*1./3600.,avg_wst/well_mixed_norm,avg_ol/well_mixed_norm,avg_h/well_mixed_norm
     638    ! do j=1,25
     639    !      !write(551,*))itime*1.,h_well/50.*j,well_mixed_vector(j)/well_mixed_norm*50.
     640    !     avg_air_dens(j)=avg_air_dens(j)/well_mixed_vector(j)
     641    !     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
     642    ! end do
     643!!------------------------------------------------------------------------------
    569644
    570645  end do
    … …  
    590665      deallocate(flux)
    591666  endif
    592   if (OHREA.eqv..TRUE.) then
    593       deallocate(OH_field,OH_field_height)
     667  if (OHREA) then
     668      deallocate(OH_field,OH_hourly,lonOH,latOH,altOH)
    594669  endif
    595670  if (ldirect.gt.0) then
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG