Ignore:
Timestamp:
Sep 24, 2018 3:42:55 PM (3 years ago)
Author:
knoop
Message:

Modularization of all bulk cloud physics code components

File:
1 edited

Legend:

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

    r3248 r3274  
    2828! -----------------
    2929! $Id$
     30! Modularization of all bulk cloud physics code components
     31!
     32! 3248 2018-09-14 09:42:06Z sward
    3033! Minor formating changes
    3134!
     
    331334#if ! defined( __nopointer )
    332335    USE arrays_3d,                                                             &
    333         ONLY:  dzu, hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
     336        ONLY:  dzu, hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend, exner
    334337#endif
    335338
    336     USE cloud_parameters,                                                      &
    337         ONLY:  cp, r_d
    338 
    339     USE constants,                                                             &
    340         ONLY:  pi
    341    
     339    USE basic_constants_and_equations_mod,                                     &
     340        ONLY:  c_p, g, kappa, pi, r_d
     341
    342342    USE control_parameters,                                                    &
    343343        ONLY:  coupling_start_time, topography, dt_3d, humidity,               &
     
    346346               timestep_scheme, tsc, coupling_char, io_blocks, io_group,       &
    347347               message_string, time_since_reference_point, surface_pressure,   &
    348                g, pt_surface, large_scale_forcing, lsf_surf, spinup,           &
     348               pt_surface, large_scale_forcing, lsf_surf, spinup,              &
    349349               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
    350                average_count_3d, varnamelength, urban_surface, kappa,          &
     350               average_count_3d, varnamelength, urban_surface,                 &
    351351               plant_canopy
    352352
     
    376376        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,        &
    377377               radiation, rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,        &
    378                sigma_sb, solar_constant, sun_direction, sun_dir_lat,           &
    379                sun_dir_lon,                                                    &
     378               sigma_sb, sun_direction, sun_dir_lat, sun_dir_lon,              &
    380379               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
    381380               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
     
    36033602        REAL(wp)     ::  ground_floor_level_l         !< local height of ground floor level
    36043603        REAL(wp)     ::  z_agl                        !< height above ground
    3605         REAL(wp), DIMENSION(nzb:nzt)   ::  exn        !< value of the Exner function in layers
    36063604
    36073605!
     
    46904688        IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.        &
    46914689             TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    4692        
    4693 !--         Calculate initial surface temperature from pt of adjacent gridbox
    4694 #if ! defined( __nopointer )
    4695             exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp          !< Exner function
    4696 #endif
    46974690
    46984691!
     
    47054698               k = surf_usm_h%k(m)
    47064699
    4707                t_surf_h(m) = pt(k,j,i) * exn(k)
    4708                t_surf_window_h(m) = pt(k,j,i) * exn(k)
    4709                t_surf_green_h(m) = pt(k,j,i) * exn(k)
    4710                surf_usm_h%pt_surface(m) = pt(k,j,i) * exn(k)
     4700               t_surf_h(m) = pt(k,j,i) * exner(k)
     4701               t_surf_window_h(m) = pt(k,j,i) * exner(k)
     4702               t_surf_green_h(m) = pt(k,j,i) * exner(k)
     4703               surf_usm_h%pt_surface(m) = pt(k,j,i) * exner(k)
    47114704            ENDDO
    47124705!
     
    47184711                  k = surf_usm_v(l)%k(m)
    47194712
    4720                   t_surf_v(l)%t(m) = pt(k,j,i) * exn(k)
    4721                   t_surf_window_v(l)%t(m) = pt(k,j,i) * exn(k)
    4722                   t_surf_green_v(l)%t(m) = pt(k,j,i) * exn(k)
    4723                   surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exn(k)
     4713                  t_surf_v(l)%t(m) = pt(k,j,i) * exner(k)
     4714                  t_surf_window_v(l)%t(m) = pt(k,j,i) * exner(k)
     4715                  t_surf_green_v(l)%t(m) = pt(k,j,i) * exner(k)
     4716                  surf_usm_v(l)%pt_surface(m) = pt(k,j,i) * exner(k)
    47244717               ENDDO
    47254718            ENDDO
     
    73387331        REAL(wp)                              :: lambda_surface_window !< current value of lambda_surface (heat conductivity between air and window)
    73397332        REAL(wp)                              :: lambda_surface_green  !< current value of lambda_surface (heat conductivity between air and greeb wall)
    7340         REAL(wp), DIMENSION(nzb:nzt)          :: exn                !< value of the Exner function in layers
    73417333       
    73427334        REAL(wp)                              :: dtime              !< simulated time of day (in UTC)
     
    73457337
    73467338
    7347 #if ! defined( __nopointer )
    7348         exn(nzb:nzt) = (hyp(nzb:nzt) / 100000.0_wp )**0.286_wp          !< Exner function
    7349 #endif
    73507339!       
    73517340!--     First, treat horizontal surface elements
     
    73717360#if ! defined( __nopointer )
    73727361!
    7373 !--        calculate rho * cp coefficient at surface layer
    7374            rho_cp  = cp * hyp(k) / ( r_d * surf_usm_h%pt1(m) * exn(k) )
     7362!--        calculate rho * c_p coefficient at surface layer
     7363           rho_cp  = c_p * hyp(k) / ( r_d * surf_usm_h%pt1(m) * exner(k) )
    73757364#endif
    73767365!
     
    73877376           surf_usm_h%r_a_green(m)  = surf_usm_h%r_a(m)
    73887377
    7389 !            r_a = ( surf_usm_h%pt1(m) - t_surf_h(m) / exn(k) ) /                              &
     7378!            r_a = ( surf_usm_h%pt1(m) - t_surf_h(m) / exner(k) ) /                              &
    73907379!                  ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-20_wp )
    7391 !            r_a_window = ( surf_usm_h%pt1(m) - t_surf_window_h(m) / exn(k) ) /                &
     7380!            r_a_window = ( surf_usm_h%pt1(m) - t_surf_window_h(m) / exner(k) ) /                &
    73927381!                  ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-20_wp )
    7393 !            r_a_green = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exn(k) ) /                  &
     7382!            r_a_green = ( surf_usm_h%pt1(m) - t_surf_green_h(m) / exner(k) ) /                  &
    73947383!                  ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-20_wp )
    73957384               
     
    74467435           coef_2 = 4.0_wp * surf_usm_h%emissivity(ind_veg_wall,m) *           &
    74477436                             sigma_sb * t_surf_h(m) ** 3                       &
    7448                            + lambda_surface + f_shf / exn(k)
     7437                           + lambda_surface + f_shf / exner(k)
    74497438           coef_window_2 = 4.0_wp * surf_usm_h%emissivity(ind_wat_win,m) *     &
    74507439                             sigma_sb * t_surf_window_h(m) ** 3                &
    7451                            + lambda_surface_window + f_shf_window / exn(k)
     7440                           + lambda_surface_window + f_shf_window / exner(k)
    74527441           coef_green_2 = 4.0_wp * surf_usm_h%emissivity(ind_pav_green,m) *    &
    74537442                             sigma_sb * t_surf_green_h(m) ** 3                 &
    7454                            + lambda_surface_green + f_shf_green / exn(k)
     7443                           + lambda_surface_green + f_shf_green / exner(k)
    74557444
    74567445!--        implicit solution when the surface layer has no heat capacity,
     
    74807469                               + surf_usm_h%frac(ind_wat_win,m) * t_surf_window_h_p(m)   &
    74817470                               + surf_usm_h%frac(ind_pav_green,m) * t_surf_green_h_p(m) )  &
    7482                                / exn(k)
     7471                               / exner(k)
    74837472                               
    74847473           IF ( humidity )  surf_usm_h%vpt_surface(m) =                        &
     
    75407529!
    75417530!--        ground/wall/roof surface heat flux
    7542            surf_usm_h%wshf_eb(m)   = - f_shf  * ( surf_usm_h%pt1(m) - t_surf_h_p(m) / exn(k) ) *               &
     7531           surf_usm_h%wshf_eb(m)   = - f_shf  * ( surf_usm_h%pt1(m) - t_surf_h_p(m) / exner(k) ) *               &
    75437532                                       surf_usm_h%frac(ind_veg_wall,m)         &
    7544                                      - f_shf_window  * ( surf_usm_h%pt1(m) - t_surf_window_h_p(m) / exn(k) ) * &
     7533                                     - f_shf_window  * ( surf_usm_h%pt1(m) - t_surf_window_h_p(m) / exner(k) ) * &
    75457534                                       surf_usm_h%frac(ind_wat_win,m)          &
    7546                                      - f_shf_green  * ( surf_usm_h%pt1(m) - t_surf_green_h_p(m) / exn(k) ) *   &
     7535                                     - f_shf_green  * ( surf_usm_h%pt1(m) - t_surf_green_h_p(m) / exner(k) ) *   &
    75477536                                       surf_usm_h%frac(ind_pav_green,m)
    75487537!           
    75497538!--        store kinematic surface heat fluxes for utilization in other processes
    75507539!--        diffusion_s, surface_layer_fluxes,...
    7551            surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / cp
     7540           surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p
    75527541
    75537542       ENDDO
     
    75737562#if ! defined( __nopointer )         
    75747563!
    7575 !--          calculate rho * cp coefficient at surface layer
    7576              rho_cp  = cp * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exn(k) )
     7564!--          calculate rho * c_p coefficient at surface layer
     7565             rho_cp  = c_p * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exner(k) )
    75777566#endif
    75787567
     
    76477636              coef_2 = 4.0_wp * surf_usm_v(l)%emissivity(ind_veg_wall,m) *     &
    76487637                                sigma_sb * t_surf_v(l)%t(m) ** 3               &
    7649                               + lambda_surface + f_shf / exn(k)
     7638                              + lambda_surface + f_shf / exner(k)
    76507639              coef_window_2 = 4.0_wp * surf_usm_v(l)%emissivity(ind_wat_win,m) *&
    76517640                                sigma_sb * t_surf_window_v(l)%t(m) ** 3        &
    7652                               + lambda_surface_window + f_shf / exn(k)
     7641                              + lambda_surface_window + f_shf / exner(k)
    76537642              coef_green_2 = 4.0_wp * surf_usm_v(l)%emissivity(ind_pav_green,m) *&
    76547643                                sigma_sb * t_surf_green_v(l)%t(m) ** 3         &
    7655                               + lambda_surface_green + f_shf / exn(k)
     7644                              + lambda_surface_green + f_shf / exner(k)
    76567645
    76577646!--           implicit solution when the surface layer has no heat capacity,
     
    76847673                                      + surf_usm_v(l)%frac(ind_wat_win,m) * t_surf_window_v_p(l)%t(m)  &
    76857674                                      + surf_usm_v(l)%frac(ind_pav_green,m) * t_surf_green_v_p(l)%t(m) ) &
    7686                                       / exn(k)
     7675                                      / exner(k)
    76877676                                     
    76887677              IF ( humidity )  surf_usm_v(l)%vpt_surface(m) =                  &
     
    77517740              surf_usm_v(l)%wshf_eb(m)   =                                     &
    77527741                 - f_shf  * ( surf_usm_v(l)%pt1(m) -                           &
    7753                  t_surf_v_p(l)%t(m) / exn(k) ) * surf_usm_v(l)%frac(ind_veg_wall,m)       &
     7742                 t_surf_v_p(l)%t(m) / exner(k) ) * surf_usm_v(l)%frac(ind_veg_wall,m)       &
    77547743                 - f_shf_window  * ( surf_usm_v(l)%pt1(m) -                    &
    7755                  t_surf_window_v_p(l)%t(m) / exn(k) ) * surf_usm_v(l)%frac(ind_wat_win,m)&
     7744                 t_surf_window_v_p(l)%t(m) / exner(k) ) * surf_usm_v(l)%frac(ind_wat_win,m)&
    77567745                 - f_shf_green  * ( surf_usm_v(l)%pt1(m) -                     &
    7757                  t_surf_green_v_p(l)%t(m) / exn(k) ) * surf_usm_v(l)%frac(ind_pav_green,m)
     7746                 t_surf_green_v_p(l)%t(m) / exner(k) ) * surf_usm_v(l)%frac(ind_pav_green,m)
    77587747
    77597748!           
    77607749!--           store kinematic surface heat fluxes for utilization in other processes
    77617750!--           diffusion_s, surface_layer_fluxes,...
    7762               surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / cp
     7751              surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / c_p
    77637752
    77647753           ENDDO
     
    77907779                         (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
    77917780                 IF ( aheat(k,j,i) > 0.0_wp )  THEN
    7792                     pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exn(k)*rho_cp*dzu(k))
     7781                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dzu(k))
    77937782                 ENDIF
    77947783              ENDIF
Note: See TracChangeset for help on using the changeset viewer.