Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (4 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r2547 r2696  
    11!> @file surface_layer_fluxes_mod.f90
    22!------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    1818!
    1919!------------------------------------------------------------------------------!
     20!
    2021! Current revisions:
    2122! ------------------
     
    2526! -----------------
    2627! $Id$
     28! - Implementation of chemistry module (FK)
     29! - Added calculation of pt1 and qv1 for all surface types. Added calculation of
     30!   pt_surface for default-type surfaces (BM)
     31! - Add flag to disable computation of qsws in case of urban surface (MS)
     32!
     33! 2547 2017-10-16 12:41:56Z schwenkel
    2734! extended by cloud_droplets option
    2835!
     
    188195!> @todo (re)move large_scale_forcing actions
    189196!> @todo check/optimize OpenMP directives
     197!> @todo simplify if conditions (which flux need to be computed in which case)
    190198!------------------------------------------------------------------------------!
    191199 MODULE surface_layer_fluxes_mod
     
    195203               drho_air_zw, rho_air_zw
    196204
     205#if defined( __chem )
     206    USE chem_modules,                                                          &
     207        ONLY:  constant_csflux, nvar
     208#endif
     209
    197210    USE cloud_parameters,                                                      &
    198211        ONLY:  l_d_cp, pt_d_t
     
    204217
    205218    USE control_parameters,                                                    &
    206         ONLY:  cloud_droplets, cloud_physics, constant_heatflux,               &
    207                constant_scalarflux, constant_waterflux, coupling_mode, g,      &
    208                humidity, ibc_e_b, ibc_pt_b, initializing_actions, kappa,       &
     219        ONLY:  air_chemistry, cloud_droplets, cloud_physics,                   &
     220               constant_heatflux, constant_scalarflux,                         &     
     221               constant_waterflux, coupling_mode, g, humidity, ibc_e_b,        &
     222               ibc_pt_b, initializing_actions, kappa,                          &
    209223               intermediate_timestep_count, intermediate_timestep_count_max,   &
    210224               land_surface, large_scale_forcing, lsf_surf,                    &
     
    291305       downward      = .FALSE.
    292306!
    293 !--    In case cloud physics is used, it is required to derive potential
    294 !--    temperature and specific humidity at first grid level from the fields pt
    295 !--    and q
    296        IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    297 !
    298 !--       First call for horizontal default-type surfaces (l=0 - upward facing,
    299 !--       l=1 - downward facing)
    300           DO  l = 0, 1
    301              IF ( surf_def_h(l)%ns >= 1  )  THEN
    302                 surf => surf_def_h(l)
    303                 CALL calc_pt_q
    304              ENDIF
    305           ENDDO
    306 !
    307 !--       Call for natural-type surfaces
    308           IF ( surf_lsm_h%ns >= 1  )  THEN
    309              surf => surf_lsm_h
     307!--    Derive potential temperature and specific humidity at first grid level
     308!--    from the fields pt and q
     309!
     310!--    First call for horizontal default-type surfaces (l=0 - upward facing,
     311!--    l=1 - downward facing)
     312       DO  l = 0, 1
     313          IF ( surf_def_h(l)%ns >= 1  )  THEN
     314             surf => surf_def_h(l)
    310315             CALL calc_pt_q
    311           ENDIF
    312 !
    313 !--       Call for urban-type surfaces
    314           IF ( surf_usm_h%ns >= 1  )  THEN
    315              surf => surf_usm_h
     316             IF ( .NOT. neutral )  CALL calc_pt_surface
     317          ENDIF
     318       ENDDO
     319!
     320!--    Call for natural-type horizontal surfaces
     321       IF ( surf_lsm_h%ns >= 1  )  THEN
     322          surf => surf_lsm_h
     323          CALL calc_pt_q
     324       ENDIF
     325
     326!
     327!--    Call for urban-type horizontal surfaces
     328       IF ( surf_usm_h%ns >= 1  )  THEN
     329          surf => surf_usm_h
     330          CALL calc_pt_q
     331       ENDIF
     332
     333!
     334!--    Call for natural-type vertical surfaces
     335       DO  l = 0, 3
     336          IF ( surf_lsm_v(l)%ns >= 1  )  THEN
     337             surf => surf_lsm_v(l)
    316338             CALL calc_pt_q
    317339          ENDIF
    318        ENDIF
     340
     341!
     342!--    Call for urban-type vertical surfaces
     343          IF ( surf_usm_v(l)%ns >= 1  )  THEN
     344             surf => surf_usm_v(l)
     345             CALL calc_pt_q
     346          ENDIF
     347       ENDDO
    319348
    320349!
     
    923952!--       However, this do not apply for downward-facing walls. To mask this,
    924953!--       use ibit, which checks for upward/downward-facing surfaces.
    925 
    926954          surf%uvw_abs(m) = SQRT(                                              &
    927955                              ( 0.5_wp * (   u(k,j,i)   + u(k,j,i+1)           &
     
    15381566             surf%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    15391567             surf%qv1(m) = q(k,j,i)
     1568          ELSE
     1569             surf%pt1(m) = pt(k,j,i)
     1570             IF ( humidity )  THEN
     1571                surf%qv1(m) = q(k,j,i)
     1572             ELSE
     1573                surf%qv1(m) = 0.0_wp
     1574             ENDIF
    15401575          ENDIF
    15411576
     
    15431578
    15441579    END SUBROUTINE calc_pt_q
     1580
     1581
     1582!
     1583!-- Calculate potential temperature and specific humidity at first grid level
     1584!-- ( only for upward-facing surfs )
     1585    SUBROUTINE calc_pt_surface
     1586
     1587       IMPLICIT NONE
     1588
     1589       INTEGER(iwp) ::  koff    !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls)
     1590       INTEGER(iwp) ::  m       !< loop variable over all horizontal surf elements
     1591       
     1592       koff = surf%koff
     1593       !$OMP PARALLEL DO PRIVATE( i, j, k )
     1594       DO  m = 1, surf%ns
     1595
     1596          i   = surf%i(m)           
     1597          j   = surf%j(m)
     1598          k   = surf%k(m)
     1599
     1600          surf%pt_surface(m) = pt(k+koff,j,i)
     1601
     1602       ENDDO
     1603
     1604    END SUBROUTINE calc_pt_surface
    15451605
    15461606!
     
    15521612
    15531613       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
    1554 
     1614       INTEGER(iwp)  ::  lsp     !< running index for chemical species
    15551615!
    15561616!--    Compute theta* at horizontal surfaces
     
    15921652          ENDIF
    15931653
    1594           IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    1595              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
    1596              DO  m = 1, surf%ns   
    1597 
    1598                 i   = surf%i(m)           
    1599                 j   = surf%j(m)
    1600                 k   = surf%k(m)
    1601 
    1602                 z_mo = surf%z_mo(m)
    1603 
    1604                 surf%ts(m) = kappa * ( surf%pt1(m) - pt(k-1,j,i) )             &
    1605                                      / ( LOG( z_mo / surf%z0h(m) )             &
    1606                                          - psi_h( z_mo / surf%ol(m) )          &
    1607                                          + psi_h( surf%z0h(m) / surf%ol(m) ) )
    1608 
    1609              ENDDO
    1610           ELSE
    1611              !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
    1612              DO  m = 1, surf%ns   
    1613 
    1614                 i   = surf%i(m)           
    1615                 j   = surf%j(m)
    1616                 k   = surf%k(m)
    1617 
    1618                 z_mo = surf%z_mo(m)
    1619 
    1620                 surf%ts(m) = kappa * ( pt(k,j,i) - pt(k-1,j,i) )               &
    1621                                      / ( LOG( z_mo / surf%z0h(m) )             &
    1622                                          - psi_h( z_mo / surf%ol(m) )          &
    1623                                          + psi_h( surf%z0h(m) / surf%ol(m) ) )
    1624              ENDDO
    1625           ENDIF
     1654          !$OMP PARALLEL DO PRIVATE( z_mo )
     1655          DO  m = 1, surf%ns   
     1656
     1657             z_mo = surf%z_mo(m)
     1658
     1659             surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) )      &
     1660                                  / ( LOG( z_mo / surf%z0h(m) )             &
     1661                                      - psi_h( z_mo / surf%ol(m) )          &
     1662                                      + psi_h( surf%z0h(m) / surf%ol(m) ) )
     1663
     1664          ENDDO
     1665
    16261666       ENDIF
    16271667!
     
    17691809
    17701810!
     1811!--    If required compute cs* (chemical species)
     1812#if defined( __chem )
     1813       IF ( air_chemistry  )  THEN 
     1814!
     1815!--       At horizontal surfaces                             
     1816          DO  lsp = 1, nvar
     1817             IF ( constant_csflux(lsp)  .AND.  .NOT.  surf_vertical )  THEN
     1818!--             For a given chemical species' flux in the surface layer
     1819                !$OMP PARALLEL DO PRIVATE( i, j )
     1820                DO  m = 1, surf%ns 
     1821                   i   = surf%i(m)           
     1822                   j   = surf%j(m)
     1823                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
     1824                ENDDO
     1825             ENDIF
     1826          ENDDO
     1827!
     1828!--       At vertical surfaces
     1829          IF ( surf_vertical )  THEN
     1830             DO  lsp = 1, nvar
     1831                !$OMP PARALLEL DO PRIVATE( i, j )
     1832                DO  m = 1, surf%ns 
     1833                   i   = surf%i(m)           
     1834                   j   = surf%j(m)
     1835                   surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )
     1836                ENDDO
     1837             ENDDO
     1838          ENDIF
     1839       ENDIF
     1840#endif
     1841
     1842!
    17711843!--    If required compute qc* and nc*
    17721844       IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
     
    18301902       IMPLICIT NONE
    18311903
    1832        INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
     1904       INTEGER(iwp)  ::  m       !< loop variable over all horizontal surf elements
     1905       INTEGER(iwp)  ::  lsp     !< running index for chemical species
    18331906
    18341907       REAL(wp)                            ::  dum     !< dummy to precalculate logarithm
     
    18831956
    18841957                surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )
    1885 !
    1886 !--        To Do: Is the sign correct???
    18871958                surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)
    18881959
     
    19292000
    19302001                surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )
    1931  
    19322002                surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)
    19332003             ENDDO
     
    19372007          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
    19382008               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
    1939                .NOT.  land_surface )  .AND.  .NOT.  urban_surface  .AND.       &
     2009               .NOT.  land_surface )  .AND.  .NOT. urban_surface  .AND.        &
    19402010               .NOT. downward )  THEN
    19412011             !$OMP PARALLEL DO PRIVATE( i, j, k )
     
    19512021          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
    19522022               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
    1953                simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface )  .AND.    &
    1954                .NOT. downward )  THEN
     2023               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface  )  .AND.   &
     2024               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
    19552025             !$OMP PARALLEL DO PRIVATE( i, j, k )
    19562026             DO  m = 1, surf%ns
     
    19742044             ENDDO
    19752045          ENDIF   
     2046!
     2047!--       Compute the vertical chemical species' flux
     2048#if defined( __chem )
     2049          DO  lsp = 1, nvar
     2050             IF (  .NOT.  constant_csflux(lsp)  .AND.  air_chemistry  .AND.    &
     2051                   .NOT.  downward )  THEN
     2052                !$OMP PARALLEL DO PRIVATE( i, j )
     2053                DO  m = 1, surf%ns   
     2054                   i     = surf%i(m)           
     2055                   j     = surf%j(m)
     2056                   surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)
     2057                ENDDO
     2058             ENDIF   
     2059          ENDDO
     2060#endif
     2061
    19762062!
    19772063!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
Note: See TracChangeset for help on using the changeset viewer.