Ignore:
Timestamp:
Apr 7, 2016 12:01:39 PM (5 years ago)
Author:
maronga
Message:

further modularization of radiation model and plant canopy model

File:
1 moved

Legend:

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

    r1823 r1826  
    1 !> @file plant_canopy_model.f90
     1!> @file plant_canopy_model_mod.f90
    22!--------------------------------------------------------------------------------!
    33! This file is part of PALM.
     
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Further modularization
    2222!
    2323! Former revisions:
     
    3535! Changes due to new module structure of the plant canopy model:
    3636!   module plant_canopy_model_mod now contains a subroutine for the
    37 !   initialization of the canopy model (init_plant_canopy),
     37!   initialization of the canopy model (pcm_init),
    3838!   limitation of the canopy drag (previously accounted for by calculation of
    3939!   a limiting canopy timestep for the determination of the maximum LES timestep
    4040!   in subroutine timestep) is now realized by the calculation of pre-tendencies
    41 !   and preliminary velocities in subroutine plant_canopy_model,
    42 !   some redundant MPI communication removed in subroutine init_plant_canopy
     41!   and preliminary velocities in subroutine pcm_tendency,
     42!   some redundant MPI communication removed in subroutine pcm_init
    4343!   (was previously in init_3d_model),
    4444!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
     
    7373! ------------
    7474!> 1) Initialization of the canopy model, e.g. construction of leaf area density
    75 !> profile (subroutine init_plant_canopy).
     75!> profile (subroutine pcm_init).
    7676!> 2) Calculation of sinks and sources of momentum, heat and scalar concentration
    77 !> due to canopy elements (subroutine plant_canopy_model).
     77!> due to canopy elements (subroutine pcm_tendency).
    7878!------------------------------------------------------------------------------!
    7979 MODULE plant_canopy_model_mod
     
    134134
    135135    PRIVATE
    136     PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,      &
    137            canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad,    &
    138            lad_s, lad_surface, lad_vertical_gradient,                          &
    139            lad_vertical_gradient_level, lad_vertical_gradient_level_ind,       &
    140            lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec,     &
    141            pch_index, plant_canopy, plant_canopy_model
    142 
    143 
    144     INTERFACE init_plant_canopy
    145        MODULE PROCEDURE init_plant_canopy
    146     END INTERFACE init_plant_canopy
    147 
    148     INTERFACE plant_canopy_model
    149        MODULE PROCEDURE plant_canopy_model
    150        MODULE PROCEDURE plant_canopy_model_ij
    151     END INTERFACE plant_canopy_model
    152 
    153 
     136 
     137!
     138!-- Public functions
     139    PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
     140
     141!
     142!-- Public variables and constants
     143    PUBLIC canopy_mode, cthf, dt_plant_canopy, plant_canopy
     144
     145
     146    INTERFACE pcm_check_parameters
     147       MODULE PROCEDURE pcm_check_parameters
     148    END INTERFACE pcm_check_parameters     
     149   
     150     INTERFACE pcm_header
     151       MODULE PROCEDURE pcm_header
     152    END INTERFACE pcm_header       
     153   
     154    INTERFACE pcm_init
     155       MODULE PROCEDURE pcm_init
     156    END INTERFACE pcm_init
     157
     158    INTERFACE pcm_parin
     159       MODULE PROCEDURE pcm_parin
     160    END INTERFACE pcm_parin 
     161   
     162    INTERFACE pcm_tendency
     163       MODULE PROCEDURE pcm_tendency
     164       MODULE PROCEDURE pcm_tendency_ij
     165    END INTERFACE pcm_tendency
    154166
    155167
    156168 CONTAINS
    157169
    158 
     170 
     171!------------------------------------------------------------------------------!
     172! Description:
     173! ------------
     174!> Check parameters routine for plant canopy model
     175!------------------------------------------------------------------------------!
     176    SUBROUTINE pcm_check_parameters
     177
     178       USE control_parameters,                                                 &
     179           ONLY: cloud_physics, message_string, microphysics_seifert
     180                 
     181   
     182       IMPLICIT NONE
     183
     184   
     185       IF ( canopy_drag_coeff == 0.0_wp )  THEN
     186          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
     187                           'coefficient & given value is canopy_drag_coeff = 0.0'
     188          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
     189       ENDIF
     190   
     191       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
     192              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
     193          message_string = 'using the beta function for the construction ' //  &
     194                           'of the leaf area density profile requires '    //  &
     195                           'both alpha_lad and beta_lad to be /= 9999999.9'
     196          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
     197       ENDIF
     198   
     199       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
     200          message_string = 'using the beta function for the construction ' //  &
     201                           'of the leaf area density profile requires '    //  &
     202                           'a non-zero lai_beta, but given value is '      //  &
     203                           'lai_beta = 0.0'
     204          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
     205       ENDIF
     206
     207       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
     208          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
     209                           'and lad_surface /= 0.0 is not possible, '       // &
     210                           'use either vertical gradients or the beta '     // &
     211                           'function for the construction of the leaf area '// &
     212                           'density profile'
     213          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
     214       ENDIF
     215
     216       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     217          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
     218                          ' seifert_beheng'
     219          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
     220       ENDIF
     221
     222 
     223    END SUBROUTINE pcm_check_parameters
     224 
     225
     226!------------------------------------------------------------------------------!
     227! Description:
     228! ------------
     229!> Header output for plant canopy model
     230!------------------------------------------------------------------------------!
     231    SUBROUTINE pcm_header ( io )
     232
     233       USE control_parameters,                                                 &
     234           ONLY: dz, passive_scalar
     235
     236
     237       IMPLICIT NONE
     238 
     239       CHARACTER (LEN=10) ::  coor_chr            !<
     240
     241       CHARACTER (LEN=86) ::  coordinates         !<
     242       CHARACTER (LEN=86) ::  gradients           !<
     243       CHARACTER (LEN=86) ::  leaf_area_density   !<
     244       CHARACTER (LEN=86) ::  slices              !<
     245 
     246       INTEGER(iwp) :: i                !<
     247       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
     248       INTEGER(iwp) :: k                !<       
     249   
     250       REAL(wp) ::  canopy_height       !< canopy height (in m)
     251       
     252       canopy_height = pch_index * dz
     253
     254       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
     255                          canopy_drag_coeff
     256       IF ( passive_scalar )  THEN
     257          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
     258                             leaf_surface_conc
     259       ENDIF
     260
     261!
     262!--    Heat flux at the top of vegetation
     263       WRITE ( io, 3 )  cthf
     264
     265!
     266!--    Leaf area density profile, calculated either from given vertical
     267!--    gradients or from beta probability density function.
     268       IF (  .NOT.  calc_beta_lad_profile )  THEN
     269
     270!--       Building output strings, starting with surface value
     271          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
     272          gradients = '------'
     273          slices = '     0'
     274          coordinates = '   0.0'
     275          i = 1
     276          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i)        &
     277                      /= -9999 )
     278
     279             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
     280             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
     281                                 TRIM( coor_chr )
     282 
     283             WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
     284             gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     285
     286             WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
     287             slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     288
     289             WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
     290             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     291
     292             i = i + 1
     293          ENDDO
     294
     295          WRITE ( io, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
     296                             TRIM( gradients ), TRIM( slices )
     297
     298       ELSE
     299       
     300          WRITE ( leaf_area_density, '(F7.4)' )  lad_surface
     301          coordinates = '   0.0'
     302         
     303          DO  k = 1, pch_index
     304
     305             WRITE (coor_chr,'(F7.2)')  lad(k)
     306             leaf_area_density = TRIM( leaf_area_density ) // ' ' //           &
     307                                 TRIM( coor_chr )
     308 
     309             WRITE (coor_chr,'(F7.1)')  zu(k)
     310             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     311
     312          ENDDO       
     313
     314          WRITE ( io, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
     315                          alpha_lad, beta_lad, lai_beta
     316
     317       ENDIF 
     318
     3191 FORMAT (//' Vegetation canopy (drag) model:'/                                &
     320              ' ------------------------------'//                              &
     321              ' Canopy mode: ', A /                                            &
     322              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
     323              ' Leaf drag coefficient: ',F6.2 /)
     3242 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
     325              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
     3263 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
     327          ' K m/s')
     3284 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
     329              ' Height:              ',A,'  m'/                                &
     330              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
     331              ' Gradient:            ',A,'  m**2/m**4'/                        &
     332              ' Gridpoint:           ',A)
     3335 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'&
     334          //  ' Height:              ',A,'  m'/                                &
     335              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
     336              ' Coefficient alpha: ',F6.2 /                                    &
     337              ' Coefficient beta: ',F6.2 /                                     &
     338              ' Leaf area index: ',F6.2,'  m**2/m**2' /)   
     339       
     340    END SUBROUTINE pcm_header
     341 
     342 
    159343!------------------------------------------------------------------------------!
    160344! Description:
     
    162346!> Initialization of the plant canopy model
    163347!------------------------------------------------------------------------------!
    164     SUBROUTINE init_plant_canopy
     348    SUBROUTINE pcm_init
    165349   
    166350
     
    187371       pre_lad = 0.0_wp
    188372
     373!
     374!--    Set flag that indicates that the lad-profile shall be calculated by using
     375!--    a beta probability density function
     376       IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )  THEN
     377          calc_beta_lad_profile = .TRUE.
     378       ENDIF
     379       
     380       
    189381!
    190382!--    Compute the profile of leaf area density used in the plant
     
    244436          DO k = nzb, pch_index
    245437             int_bpdf = int_bpdf +                                             &
    246                            ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * &
    247                            ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) *   &
    248                            ( ( zw(k+1)-zw(k) ) / canopy_height ) )
     438                      ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
     439                      ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(              &
     440                          beta_lad-1.0_wp ) )                                  &
     441                      * ( ( zw(k+1)-zw(k) ) / canopy_height ) )
    249442          ENDDO
    250443
     
    252445!--       Preliminary lad profile (defined on w-grid)
    253446          DO k = nzb, pch_index
    254              pre_lad(k) =  lai_beta *                                              &
    255                            (   ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *          &
    256                                ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) /   &
    257                                int_bpdf                                            &
    258                            ) / canopy_height
     447             pre_lad(k) =  lai_beta *                                          &
     448                        ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
     449                        * ( ( 1.0_wp - ( zw(k) / canopy_height ) )**(          &
     450                              beta_lad-1.0_wp ) ) / int_bpdf                   &
     451                        ) / canopy_height
    259452          ENDDO
    260453
     
    376569
    377570
    378     END SUBROUTINE init_plant_canopy
     571    END SUBROUTINE pcm_init
    379572
    380573
     
    387580!>
    388581!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
    389 !> (defined on scalar grid), as initialized in subroutine init_plant_canopy.
     582!> (defined on scalar grid), as initialized in subroutine pcm_init.
    390583!> The lad on the w-grid is vertically interpolated from the surrounding
    391584!> lad_s. The upper boundary of the canopy is defined on the w-grid at
     
    404597!> Call for all grid points
    405598!------------------------------------------------------------------------------!
    406     SUBROUTINE plant_canopy_model( component )
     599    SUBROUTINE pcm_tendency( component )
    407600
    408601
     
    672865
    673866             WRITE( message_string, * ) 'wrong component: ', component
    674              CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
     867             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
    675868
    676869       END SELECT
    677870
    678     END SUBROUTINE plant_canopy_model
     871    END SUBROUTINE pcm_tendency
     872
     873
     874!------------------------------------------------------------------------------!
     875! Description:
     876! ------------
     877!> Parin for &canopy_par for plant canopy model
     878!------------------------------------------------------------------------------!
     879    SUBROUTINE pcm_parin
     880
     881
     882       IMPLICIT NONE
     883
     884       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
     885       
     886        NAMELIST /canopy_par/     alpha_lad, beta_lad, canopy_drag_coeff,      &
     887                                  canopy_mode, cthf,                           &
     888                                  lad_surface,                                 &
     889                                  lad_vertical_gradient,                       &
     890                                  lad_vertical_gradient_level,                 &
     891                                  lai_beta,                                    &
     892                                  leaf_scalar_exch_coeff,                      &
     893                                  leaf_surface_conc, pch_index
     894       
     895       line = ' '
     896       
     897!
     898!--    Try to find radiation model package
     899       REWIND ( 11 )
     900       line = ' '
     901       DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
     902          READ ( 11, '(A)', END=10 )  line
     903       ENDDO
     904       BACKSPACE ( 11 )
     905
     906!
     907!--    Read user-defined namelist
     908       READ ( 11, canopy_par )
     909
     910!
     911!--    Set flag that indicates that the radiation model is switched on
     912       plant_canopy = .TRUE.
     913
     914 10    CONTINUE
     915       
     916
     917    END SUBROUTINE pcm_parin   
     918   
    679919
    680920
     
    686926!>
    687927!> The canopy is located where the leaf area density lad_s(k,j,i) > 0.0
    688 !> (defined on scalar grid), as initialized in subroutine init_plant_canopy.
     928!> (defined on scalar grid), as initialized in subroutine pcm_init.
    689929!> The lad on the w-grid is vertically interpolated from the surrounding
    690930!> lad_s. The upper boundary of the canopy is defined on the w-grid at
     
    703943!> Call for grid point i,j
    704944!------------------------------------------------------------------------------!
    705     SUBROUTINE plant_canopy_model_ij( i, j, component )
     945    SUBROUTINE pcm_tendency_ij( i, j, component )
    706946
    707947
     
    9441184
    9451185          WRITE( message_string, * ) 'wrong component: ', component
    946           CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 )
     1186          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 )
    9471187
    9481188       END SELECT
    9491189
    950     END SUBROUTINE plant_canopy_model_ij
     1190    END SUBROUTINE pcm_tendency_ij
    9511191
    9521192 END MODULE plant_canopy_model_mod
Note: See TracChangeset for help on using the changeset viewer.