Changeset 4230 for palm/trunk


Ignore:
Timestamp:
Sep 11, 2019 1:58:14 PM (23 months ago)
Author:
suehring
Message:

Several bugfixes: profile initialization of chemical species in restart runs; Runge-Kutta tendency array not initialized in chemistry model in restart runs; fixed determination of time indices for chemical emissions (introduced with commit -4227); Update chemistry profiles in offline nesting; initialize canopy resistances for greened building walls (even if green fraction is zero)

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/chem_emissions_mod.f90

    r4227 r4230  
    2727! -----------------
    2828! $Id$
     29! Bugfix, consider that time_since_reference_point can be also negative when
     30! time indices are determined.
     31!
     32! 4227 2019-09-10 18:04:34Z gronemeier
    2933! implement new palm_date_time_mod
    3034!
     
    10911095
    10921096          CALL get_date_time( 0.0_wp, second_of_day=time_utc_init )
    1093           CALL get_date_time( time_since_reference_point, hour=hour_of_day )
    1094 
    1095           days_since_reference_point = INT( FLOOR( ( time_utc_init + time_since_reference_point ) &
     1097          CALL get_date_time( MAX( 0.0_wp, time_since_reference_point ),       &
     1098                              hour=hour_of_day )
     1099
     1100          days_since_reference_point = INT( FLOOR( (                            &
     1101                    time_utc_init + MAX( 0.0_wp, time_since_reference_point ) ) &
    10961102                                                   / seconds_per_day ) )
    10971103
     
    11241130!-- Update time indices
    11251131
    1126              CALL get_date_time( time_since_reference_point, &
     1132             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),    &
    11271133                                 day_of_year=day_of_year, hour=hour_of_day )
    11281134             index_hh = ( day_of_year - 1_iwp ) * hour_of_day
     
    11561162!
    11571163!-- Update time indices
    1158              CALL get_date_time( time_since_reference_point, &
    1159                                  month=month_of_year,        &
    1160                                  day=day_of_month,           &
    1161                                  hour=hour_of_day,           &
     1164             CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ), &
     1165                                 month=month_of_year,                       &
     1166                                 day=day_of_month,                          &
     1167                                 hour=hour_of_day,                          &
    11621168                                 day_of_week=day_of_week     )
    11631169             index_mm = month_of_year
     
    12301236!
    12311237!--       Get time-factor for specific hour
    1232           CALL get_date_time( time_since_reference_point, hour=hour_of_day )
     1238          CALL get_date_time( MAX( time_since_reference_point, 0.0_wp ),              &
     1239                              hour=hour_of_day )
    12331240
    12341241          index_hh = hour_of_day
  • palm/trunk/SOURCE/chemistry_model_mod.f90

    r4227 r4230  
    2727! -----------------
    2828! $Id$
     29! Bugfix, initialize mean profiles also in restart runs. Also initialize
     30! array used for Runge-Kutta tendecies in restart runs. 
     31!
     32! 4227 2019-09-10 18:04:34Z gronemeier
    2933! implement new palm_date_time_mod
    3034!
     
    19481952          ENDDO
    19491953       ENDIF
    1950 !
    1951 !--    Initiale old and new time levels.
    1952        DO  lsp = 1, nvar
    1953           chem_species(lsp)%tconc_m = 0.0_wp                     
    1954           chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
    1955        ENDDO
    19561954
    19571955    ENDIF
     1956!
     1957!-- Initiale old and new time levels. Note, this has to be done also in restart runs
     1958    DO  lsp = 1, nvar
     1959       chem_species(lsp)%tconc_m = 0.0_wp                     
     1960       chem_species(lsp)%conc_p  = chem_species(lsp)%conc     
     1961    ENDDO
    19581962
    19591963    DO  lsp = 1, nphot
     
    19641968!--          WRITE(6,'(a,i4,3x,a)')  'Photolysis: ',lsp,TRIM( phot_names(lsp) )
    19651969!--       ENDIF
    1966           phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
     1970       phot_frequen(lsp)%freq(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  =>  freq_1(:,:,:,lsp)
    19671971    ENDDO
    19681972
     
    19821986!------------------------------------------------------------------------------!
    19831987 SUBROUTINE chem_init_profiles             
    1984 !
    1985 !-- SUBROUTINE is called from chem_init in case of TRIM( initializing_actions ) /= 'read_restart_data'
    1986 !< We still need to see what has to be done in case of restart run
    19871988
    19881989    USE chem_modules
     
    20002001!-- "cs_heights" are prescribed, their values will!override the constant profile given by
    20012002!-- "cs_surface".
    2002     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    2003        lsp_usr = 1
    2004        DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
    2005           DO  lsp = 1, nspec                                !
    2006 !
    2007 !--          create initial profile (conc_pr_init) for each chemical species
    2008              IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
    2009                 IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
    2010 !
    2011 !--               set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
    2012                    DO lpr_lev = 0, nzt+1
    2013                       chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
    2014                    ENDDO
    2015                 ELSE
    2016                    IF ( cs_heights(1,1) /= 0.0_wp )  THEN
    2017                       message_string = 'The surface value of cs_heights must be 0.0'
    2018                       CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
     2003!     IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     2004    lsp_usr = 1
     2005    DO  WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )   !'novalue' is the default
     2006       DO  lsp = 1, nspec                                !
     2007!
     2008!--       create initial profile (conc_pr_init) for each chemical species
     2009          IF ( TRIM( chem_species(lsp)%name ) == TRIM( cs_name(lsp_usr) ) )  THEN   !
     2010             IF ( cs_profile(lsp_usr,1) == 9999999.9_wp )  THEN
     2011!
     2012!--            set a vertically constant profile based on the surface conc (cs_surface(lsp_usr)) of each species
     2013                DO lpr_lev = 0, nzt+1
     2014                   chem_species(lsp)%conc_pr_init(lpr_lev) = cs_surface(lsp_usr)
     2015                ENDDO
     2016             ELSE
     2017                IF ( cs_heights(1,1) /= 0.0_wp )  THEN
     2018                   message_string = 'The surface value of cs_heights must be 0.0'
     2019                   CALL message( 'chem_check_parameters', 'CM0434', 1, 2, 0, 6, 0 )
     2020                ENDIF
     2021
     2022                use_prescribed_profile_data = .TRUE.
     2023
     2024                npr_lev = 1
     2025!                chem_species(lsp)%conc_pr_init(0) = 0.0_wp
     2026                DO  lpr_lev = 1, nz+1
     2027                   IF ( npr_lev < 100 )  THEN
     2028                      DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
     2029                         npr_lev = npr_lev + 1
     2030                         IF ( npr_lev == 100 )  THEN
     2031                            message_string = 'number of chem spcs exceeding the limit'
     2032                            CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
     2033                            EXIT
     2034                         ENDIF
     2035                      ENDDO
    20192036                   ENDIF
    2020 
    2021                    use_prescribed_profile_data = .TRUE.
    2022 
    2023                    npr_lev = 1
    2024 !                   chem_species(lsp)%conc_pr_init(0) = 0.0_wp
    2025                    DO  lpr_lev = 1, nz+1
    2026                       IF ( npr_lev < 100 )  THEN
    2027                          DO  WHILE ( cs_heights(lsp_usr, npr_lev+1) <= zu(lpr_lev) )
    2028                             npr_lev = npr_lev + 1
    2029                             IF ( npr_lev == 100 )  THEN
    2030                                message_string = 'number of chem spcs exceeding the limit'
    2031                                CALL message( 'chem_check_parameters', 'CM0435', 1, 2, 0, 6, 0 )               
    2032                                EXIT
    2033                             ENDIF
    2034                          ENDDO
    2035                       ENDIF
    2036                       IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
    2037                          chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
    2038                               ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
    2039                               ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
    2040                               ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
    2041                       ELSE
    2042                          chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
    2043                       ENDIF
    2044                    ENDDO
    2045                 ENDIF
    2046 !
    2047 !--          If a profile is prescribed explicity using cs_profiles and cs_heights, then 
    2048 !--          chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
    2049 !--          on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     2037                   IF ( npr_lev < 100  .AND.  cs_heights(lsp_usr,npr_lev+1) /= 9999999.9_wp )  THEN
     2038                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev) +       &
     2039                           ( zu(lpr_lev) - cs_heights(lsp_usr, npr_lev) ) /                          &
     2040                           ( cs_heights(lsp_usr, (npr_lev + 1)) - cs_heights(lsp_usr, npr_lev ) ) *  &
     2041                           ( cs_profile(lsp_usr, (npr_lev + 1)) - cs_profile(lsp_usr, npr_lev ) )
     2042                   ELSE
     2043                      chem_species(lsp)%conc_pr_init(lpr_lev) = cs_profile(lsp_usr, npr_lev)
     2044                   ENDIF
     2045                ENDDO
    20502046             ENDIF
    2051           ENDDO
    2052           lsp_usr = lsp_usr + 1
     2047!
     2048!--       If a profile is prescribed explicity using cs_profiles and cs_heights, then 
     2049!--       chem_species(lsp)%conc_pr_init is populated with the specific "lsp" based
     2050!--       on the cs_profiles(lsp_usr,:)  and cs_heights(lsp_usr,:).
     2051          ENDIF
     2052
    20532053       ENDDO
    2054     ENDIF
     2054
     2055       lsp_usr = lsp_usr + 1
     2056    ENDDO
     2057!     ENDIF
    20552058
    20562059 END SUBROUTINE chem_init_profiles
  • palm/trunk/SOURCE/nesting_offl_mod.f90

    r4227 r4230  
    2020! Current revisions:
    2121! ------------------
    22 ! implement new palm_date_time_mod
     22!
    2323!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Update mean chemistry profiles. These are also used for rayleigh damping.
     28!
     29! 4227 2019-09-10 18:04:34Z gronemeier
     30! implement new palm_date_time_mod
     31!
    2732! - Data input moved into nesting_offl_mod
    2833! - check rephrased
     
    403408!--    Check if dynamic driver data input is required.
    404409       IF ( nest_offl%time(nest_offl%tind_p) <=                                &
    405             MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.                   &
     410            MAX( time_since_reference_point, 0.0_wp) + time_utc_init  .OR.     &
    406411            .NOT.  nest_offl%init )  THEN
    407412          CONTINUE
     
    922927       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref   !< reference profile for potential temperature
    923928       REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_ref_l !< reference profile for potential temperature on subdomain
    924        REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio on subdomain
     929       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref    !< reference profile for mixing ratio
    925930       REAL(wp), DIMENSION(nzb:nzt+1) ::  q_ref_l  !< reference profile for mixing ratio on subdomain
    926        REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component on subdomain
     931       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref    !< reference profile for u-component
    927932       REAL(wp), DIMENSION(nzb:nzt+1) ::  u_ref_l  !< reference profile for u-component on subdomain
    928        REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component on subdomain
     933       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref    !< reference profile for v-component
    929934       REAL(wp), DIMENSION(nzb:nzt+1) ::  v_ref_l  !< reference profile for v-component on subdomain
    930        
     935
     936       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem   !< reference profile for chemical species
     937       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ref_chem_l !< reference profile for chemical species on subdomain
    931938
    932939       IF ( debug_output_timestep )  CALL debug_message( 'nesting_offl_bc', 'start' )
     
    934941       CALL  cpu_log( log_point(58), 'offline nesting', 'start' )     
    935942!
    936 !--    Set mean profiles, derived from boundary data, to zero
     943!--    Initialize mean profiles, derived from boundary data, to zero
    937944       pt_ref   = 0.0_wp
    938945       q_ref    = 0.0_wp
     
    944951       u_ref_l  = 0.0_wp
    945952       v_ref_l  = 0.0_wp
     953!
     954!--    If required, allocate temporary arrays to compute chemistry mean profiles
     955       IF ( air_chemistry )  THEN
     956          ALLOCATE( ref_chem(nzb:nzt+1,1:UBOUND( chem_species, 1 ) )   )
     957          ALLOCATE( ref_chem_l(nzb:nzt+1,1:UBOUND( chem_species, 1 ) ) )
     958          ref_chem   = 0.0_wp
     959          ref_chem_l = 0.0_wp
     960       ENDIF
    946961!
    947962!--    Set boundary conditions of u-, v-, w-component, as well as q, and pt.
     
    10241039                                                  fac_dt                   )
    10251040                      ENDDO
     1041                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
     1042                                         + chem_species(n)%conc(nzb+1:nzt,j,-1)
    10261043                   ENDDO
    10271044                ENDIF
     
    11001117                                                 fac_dt                       )
    11011118                      ENDDO
     1119                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
     1120                                       + chem_species(n)%conc(nzb+1:nzt,j,nxr+1)
    11021121                   ENDDO
    11031122                ENDIF
     
    11791198                                                 fac_dt                    )
    11801199                      ENDDO
     1200                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
     1201                                       + chem_species(n)%conc(nzb+1:nzt,-1,i)
    11811202                   ENDDO
    11821203                ENDIF
     
    12571278                                                 fac_dt                       )
    12581279                      ENDDO
     1280                      ref_chem_l(nzb+1:nzt,n) = ref_chem_l(nzb+1:nzt,n)        &
     1281                                       + chem_species(n)%conc(nzb+1:nzt,nyn+1,i)
    12591282                   ENDDO
    12601283                ENDIF
     
    13371360                   DO  j = nys, nyn
    13381361                      chem_species(n)%conc(nzt+1,j,i) = interpolate_in_time(   &
    1339                                               nest_offl%chem_top(0,j,i,n),   &
    1340                                               nest_offl%chem_top(1,j,i,n),   &
     1362                                              nest_offl%chem_top(0,j,i,n),     &
     1363                                              nest_offl%chem_top(1,j,i,n),     &
    13411364                                              fac_dt                       )
     1365                      ref_chem_l(nzt+1,n) = ref_chem_l(nzt+1,n) +              &
     1366                                            chem_species(n)%conc(nzt+1,j,i)
    13421367                   ENDDO
    13431368                ENDDO
     
    14041429                              comm2d, ierr )
    14051430       ENDIF
     1431       IF ( air_chemistry )  THEN
     1432          CALL MPI_ALLREDUCE( ref_chem_l, ref_chem,                            &
     1433                              ( nzt+1-nzb+1 ) * SIZE( ref_chem(nzb,:) ),       &
     1434                              MPI_REAL, MPI_SUM, comm2d, ierr )
     1435       ENDIF
    14061436#else
    14071437       u_ref  = u_ref_l
    14081438       v_ref  = v_ref_l
    1409        IF ( humidity )       q_ref  = q_ref_l
    1410        IF ( .NOT. neutral )  pt_ref = pt_ref_l
     1439       IF ( humidity )       q_ref    = q_ref_l
     1440       IF ( .NOT. neutral )  pt_ref   = pt_ref_l
     1441       IF ( air_chemistry )  ref_chem = ref_chem_l
    14111442#endif
    14121443!
     
    14271458                                                          ( ny + 1 + nx + 1 ), &
    14281459                                              KIND = wp )
     1460       IF ( air_chemistry )                                                    &
     1461          ref_chem(nzb:nzt,:) = ref_chem(nzb:nzt,:) / REAL( 2.0_wp *           &
     1462                                                          ( ny + 1 + nx + 1 ), &
     1463                                                            KIND = wp )
    14291464!
    14301465!--    Derived from top boundary.   
     
    14371472          pt_ref(nzt+1) = pt_ref(nzt+1) / REAL( ( ny + 1 ) * ( nx + 1 ),       &
    14381473                                                KIND = wp )
    1439 !
    1440 !--    Write onto init profiles, which are used for damping
     1474       IF ( air_chemistry )                                                    &
     1475          ref_chem(nzt+1,:) = ref_chem(nzt+1,:) /                              &
     1476                              REAL( ( ny + 1 ) * ( nx + 1 ),KIND = wp )
     1477!
     1478!--    Write onto init profiles, which are used for damping. Also set lower
     1479!--    boundary condition for scalars (not required for u and v as these are
     1480!--    zero at k=nzb.
    14411481       u_init = u_ref
    14421482       v_init = v_ref
    1443        IF ( humidity      )  q_init  = q_ref
    1444        IF ( .NOT. neutral )  pt_init = pt_ref
    1445 !
    1446 !--    Set bottom boundary condition
    1447        IF ( humidity      )  q_init(nzb)  = q_init(nzb+1)
    1448        IF ( .NOT. neutral )  pt_init(nzb) = pt_init(nzb+1)
     1483       IF ( humidity      )  THEN
     1484          q_init      = q_ref
     1485          q_init(nzb) = q_init(nzb+1)
     1486       ENDIF
     1487       IF ( .NOT. neutral )  THEN
     1488          pt_init      = pt_ref
     1489          pt_init(nzb) = pt_init(nzb+1)
     1490       ENDIF
     1491
     1492       IF ( air_chemistry )  THEN
     1493          DO  n = 1, UBOUND( chem_species, 1 )
     1494             IF ( nest_offl%chem_from_file_t(n) )  THEN
     1495                chem_species(n)%conc_pr_init(:) = ref_chem(:,n)
     1496                chem_species(n)%conc_pr_init(nzb) =                            &
     1497                                            chem_species(n)%conc_pr_init(nzb+1)
     1498             ENDIF
     1499          ENDDO
     1500       ENDIF
     1501
     1502       DEALLOCATE( ref_chem   )
     1503       DEALLOCATE( ref_chem_l )     
    14491504!
    14501505!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r4227 r4230  
    2828! -----------------
    2929! $Id$
     30! Bugfix, initialize canopy resistance. Even if no green fraction is set,
     31! r_canopy must be initialized for output purposes.
     32!
     33! 4227 2019-09-10 18:04:34Z gronemeier
    3034! implement new palm_date_time_mod
    3135!
     
    35553559!--    Map values onto horizontal elemements
    35563560       DO  m = 1, surf_usm_h%ns
    3557              surf_usm_h%r_canopy_min(m)     = 200.0_wp !< min_canopy_resistance
    3558              surf_usm_h%g_d(m)              = 0.0_wp   !< canopy_resistance_coefficient
     3561          surf_usm_h%r_canopy(m)     = 200.0_wp !< canopy_resistance
     3562          surf_usm_h%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
     3563          surf_usm_h%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
    35593564       ENDDO
    35603565!
     
    35633568       DO  l = 0, 3
    35643569          DO  m = 1, surf_usm_v(l)%ns
    3565                 surf_usm_v(l)%r_canopy_min(m)     = 200.0_wp !< min_canopy_resistance
    3566                 surf_usm_v(l)%g_d(m)              = 0.0_wp   !< canopy_resistance_coefficient
     3570             surf_usm_v(l)%r_canopy(m)     = 200.0_wp !< canopy_resistance
     3571             surf_usm_v(l)%r_canopy_min(m) = 200.0_wp !< min_canopy_resistance
     3572             surf_usm_v(l)%g_d(m)          = 0.0_wp   !< canopy_resistance_coefficient
    35673573          ENDDO
    35683574       ENDDO
     
    75857591!--           vegetated surfaces and bare soils.
    75867592              m_liq_max = m_max_depth * ( surf_usm_h%lai(m) )
     7593
    75877594              surf_usm_h%c_liq(m) = MIN( 1.0_wp, ( m_liq_usm_h%var_usm_1d(m) / m_liq_max )**0.67 )
    75887595!
     
    79877994 
    79887995
    7989               IF ( surf_usm_v(l)%frac(1,m) > 0.0_wp ) THEN
     7996              IF ( surf_usm_v(l)%frac(ind_pav_green,m) > 0.0_wp ) THEN
    79907997!
    79917998!--             Adapted from LSM:
Note: See TracChangeset for help on using the changeset viewer.