Changeset 3274 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- Sep 24, 2018 3:42:55 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (37 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r3272 r3274 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Modularization of all bulk cloud physics code components 31 ! 32 ! 3272 2018-09-24 10:16:32Z suehring 30 33 ! - split direct and diffusion shortwave radiation using RRTMG rather than using 31 34 ! calc_diffusion_radiation, in case of RRTMG … … 103 106 ! 104 107 ! 3117 2018-07-11 09:59:11Z maronga 105 ! Bugfix: water vapor was not transfered to RRTMG when cloud_physics= .F.108 ! Bugfix: water vapor was not transfered to RRTMG when bulk_cloud_model = .F. 106 109 ! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim) 107 110 ! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere … … 437 440 438 441 USE arrays_3d, & 439 ONLY: dzw, hyp, nc, pt, q, ql, zu, zw 442 ONLY: dzw, hyp, nc, pt, q, ql, zu, zw, exner, d_exner 443 444 USE basic_constants_and_equations_mod, & 445 ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, & 446 barometric_formula 440 447 441 448 USE calc_mean_profile_mod, & 442 449 ONLY: calc_mean_profile 443 450 444 USE cloud_parameters, &445 ONLY: cp, l_d_cp, l_v, r_d, rho_l446 447 USE constants, &448 ONLY: pi449 450 451 USE control_parameters, & 451 ONLY: cloud_droplets, c loud_physics, coupling_char, dz, g,&452 ONLY: cloud_droplets, coupling_char, dz, & 452 453 humidity, & 453 454 initializing_actions, io_blocks, io_group, & 454 455 latitude, longitude, large_scale_forcing, lsf_surf, & 455 message_string, microphysics_morrison,plant_canopy, pt_surface,&456 message_string, plant_canopy, pt_surface,& 456 457 rho_surface, surface_pressure, time_since_reference_point, & 457 458 urban_surface, land_surface, end_time, spinup_time, dt_spinup … … 475 476 USE kinds 476 477 477 USE microphysics_mod,&478 ONLY: na_init, nc_const, sigma_gc478 USE bulk_cloud_model_mod, & 479 ONLY: bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc 479 480 480 481 #if defined ( __netcdf ) … … 561 562 /) 562 563 564 REAL(wp), PARAMETER :: sigma_sb = 5.67037321E-8_wp !< Stefan-Boltzmann constant 565 563 566 INTEGER(iwp) :: albedo_type = 9999999, & !< Albedo surface type 564 567 dots_rad = 0 !< starting index for timeseries output … … 578 581 !< will be considered. However fewer SVFs are expected. 579 582 radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees 580 581 582 REAL(wp), PARAMETER :: sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant583 solar_constant = 1368.0_wp !< solar constant at top of atmosphere584 583 585 584 REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha … … 2274 2273 ! 2275 2274 !-- Initialize RRTMG 2276 IF ( lw_radiation ) CALL rrtmg_lw_ini ( c p )2277 IF ( sw_radiation ) CALL rrtmg_sw_ini ( c p )2275 IF ( lw_radiation ) CALL rrtmg_lw_ini ( c_p ) 2276 IF ( sw_radiation ) CALL rrtmg_sw_ini ( c_p ) 2278 2277 2279 2278 ! … … 2335 2334 2336 2335 INTEGER(iwp) :: l !< running index for surface orientation 2337 2338 REAL(wp) :: exn !< Exner functions at surface2339 REAL(wp) :: exn1 !< Exner functions at first grid level or at urban layer top2340 2336 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 2341 2337 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain … … 2355 2351 ! 2356 2352 !-- Calculate value of the Exner function at model surface 2357 exn = (surface_pressure / 1000.0_wp )**0.286_wp2358 2353 ! 2359 2354 !-- In case averaged radiation is used, calculate mean temperature and … … 2361 2356 IF ( average_radiation ) THEN 2362 2357 pt1 = 0.0_wp 2363 IF ( cloud_physics.OR. cloud_droplets ) ql1 = 0.0_wp2358 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp 2364 2359 2365 2360 pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) ) 2366 IF ( cloud_physics.OR. cloud_droplets ) ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )2361 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) ) 2367 2362 2368 2363 #if defined( __parallel ) 2369 2364 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 2370 2365 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2371 IF ( cloud_physics.OR. cloud_droplets ) &2366 IF ( bulk_cloud_model .OR. cloud_droplets ) & 2372 2367 CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2373 2368 #else 2374 2369 pt1 = pt1_l 2375 IF ( cloud_physics.OR. cloud_droplets ) ql1 = ql1_l2370 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l 2376 2371 #endif 2377 2378 exn1 = ( hyp(nzut) / 100000.0_wp )**0.286_wp 2379 IF ( cloud_physics .OR. cloud_droplets ) pt1 = pt1 + l_d_cp / exn1 * ql1 2372 2373 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nzut) * ql1 2380 2374 ! 2381 2375 !-- Finally, divide by number of grid points … … 2418 2412 k = nzut 2419 2413 2420 exn1 = ( hyp(k+1) / 100000.0_wp )**0.286_wp2421 2422 2414 surf%rad_sw_in = solar_constant * sky_trans * zenith(0) 2423 2415 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 2424 2416 2425 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exn 1)**42417 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(k+1))**4 2426 2418 2427 2419 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 2443 2435 j = surf%j(m) 2444 2436 k = surf%k(m) 2445 2446 exn1 = (hyp(k) / 100000.0_wp )**0.286_wp2447 2437 2448 2438 surf%rad_sw_in(m) = solar_constant * sky_trans * zenith(0) … … 2469 2459 ) & 2470 2460 * sigma_sb & 2471 * ( surf%pt_surface(m) * exn )**42461 * ( surf%pt_surface(m) * exner(nzb) )**4 2472 2462 2473 2463 surf%rad_lw_out_change_0(m) = & … … 2479 2469 surf%emissivity(ind_wat_win,m) & 2480 2470 ) * 3.0_wp * sigma_sb & 2481 * ( surf%pt_surface(m) * exn )** 32482 2483 2484 IF ( cloud_physics.OR. cloud_droplets ) THEN2485 pt1 = pt(k,j,i) + l _d_cp / exn1* ql(k,j,i)2486 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exn 1)**42471 * ( surf%pt_surface(m) * exner(nzb) )** 3 2472 2473 2474 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 2475 pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i) 2476 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exner(k))**4 2487 2477 ELSE 2488 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt(k,j,i) * exn 1)**42478 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt(k,j,i) * exner(k))**4 2489 2479 ENDIF 2490 2480 … … 2524 2514 INTEGER(iwp) :: l !< running index for surface orientation 2525 2515 2526 REAL(wp) :: exn !< Exner functions at surface2527 REAL(wp) :: exn1 !< Exner functions at first grid level2528 2516 REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top 2529 2517 REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain … … 2531 2519 REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain 2532 2520 2533 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 2534 2535 ! 2536 !-- Calculate value of the Exner function 2537 exn = (surface_pressure / 1000.0_wp )**0.286_wp 2521 TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine 2522 2538 2523 ! 2539 2524 !-- In case averaged radiation is used, calculate mean temperature and … … 2541 2526 IF ( average_radiation ) THEN 2542 2527 pt1 = 0.0_wp 2543 IF ( cloud_physics.OR. cloud_droplets ) ql1 = 0.0_wp2528 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp 2544 2529 2545 2530 pt1_l = SUM( pt(nzut,nys:nyn,nxl:nxr) ) 2546 IF ( cloud_physics.OR. cloud_droplets ) ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) )2531 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nzut,nys:nyn,nxl:nxr) ) 2547 2532 2548 2533 #if defined( __parallel ) 2549 2534 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 2550 2535 CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2551 IF ( cloud_physics.OR. cloud_droplets ) &2536 IF ( bulk_cloud_model .OR. cloud_droplets ) & 2552 2537 CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2553 2538 #else 2554 2539 pt1 = pt1_l 2555 IF ( cloud_physics.OR. cloud_droplets ) ql1 = ql1_l2540 IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l 2556 2541 #endif 2557 IF ( cloud_physics .OR. cloud_droplets ) pt1 = pt1 + l_d_cp / exn1* ql12542 IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nzut+1) * ql1 2558 2543 ! 2559 2544 !-- Finally, divide by number of grid points … … 2595 2580 IF ( average_radiation ) THEN 2596 2581 2597 ! set height above canopy2598 k = nzut2599 2600 2582 surf%rad_net = net_radiation 2601 ! MS: Wyh k + 1 ? 2602 exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp 2603 2604 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exn1)**4 2583 2584 surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(nzut+1))**4 2605 2585 2606 2586 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 2636 2616 surf%rad_net(m) = net_radiation 2637 2617 2638 exn1 = (hyp(k) / 100000.0_wp )**0.286_wp 2639 2640 IF ( cloud_physics .OR. cloud_droplets ) THEN 2641 pt1 = pt(k,j,i) + l_d_cp / exn1 * ql(k,j,i) 2642 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exn1)**4 2618 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 2619 pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i) 2620 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exner(k))**4 2643 2621 ELSE 2644 2622 surf%rad_lw_in(m) = 0.8_wp * sigma_sb * & 2645 ( pt(k,j,i) * exn 1)**42623 ( pt(k,j,i) * exner(k) )**4 2646 2624 ENDIF 2647 2625 … … 2656 2634 ) & 2657 2635 * sigma_sb & 2658 * ( surf%pt_surface(m) * exn )**42636 * ( surf%pt_surface(m) * exner(nzb) )**4 2659 2637 2660 2638 surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m) & … … 2954 2932 rrtm_tlev(0,nzb+1) = t_rad_urb 2955 2933 2956 IF ( cloud_physics) THEN2934 IF ( bulk_cloud_model ) THEN 2957 2935 2958 2936 CALL calc_mean_profile( ql, 54 ) … … 2962 2940 DO k = nzb+1, nzt+1 2963 2941 rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp & 2964 )**.286_wp + l _d_cp * ql_av(k)2942 )**.286_wp + lv_d_cp * ql_av(k) 2965 2943 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k)) 2966 2944 ENDDO … … 3013 2991 rrtm_icld = 0 3014 2992 3015 IF ( cloud_physics) THEN2993 IF ( bulk_cloud_model ) THEN 3016 2994 DO k = nzb+1, nzt+1 3017 2995 rrtm_cliqwp(0,k) = ql_av(k) * 1000._wp * & … … 3133 3111 !-- Prepare profiles of temperature and H2O volume mixing ratio 3134 3112 DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i) 3135 rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * & 3136 ( surface_pressure / 1000.0_wp )**0.286_wp 3113 rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * exner(nzb) 3137 3114 ENDDO 3138 3115 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) 3139 rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * & 3140 ( surface_pressure / 1000.0_wp )**0.286_wp 3116 rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * exner(nzb) 3141 3117 ENDDO 3142 3118 3143 3119 3144 IF ( cloud_physics) THEN3120 IF ( bulk_cloud_model ) THEN 3145 3121 DO k = nzb+1, nzt+1 3146 rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp&3147 )**0.286_wp + l_d_cp * ql(k,j,i)3122 rrtm_tlay(0,k) = pt(k,j,i) * exner(k) & 3123 + lv_d_cp * ql(k,j,i) 3148 3124 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i)) 3149 3125 ENDDO 3150 3126 ELSEIF ( cloud_droplets ) THEN 3151 3127 DO k = nzb+1, nzt+1 3152 rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp&3153 )**0.286_wp + l_d_cp * ql(k,j,i)3128 rrtm_tlay(0,k) = pt(k,j,i) * exner(k) & 3129 + lv_d_cp * ql(k,j,i) 3154 3130 rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i) 3155 3131 ENDDO 3156 3132 ELSE 3157 3133 DO k = nzb+1, nzt+1 3158 rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp & 3159 )**0.286_wp 3134 rrtm_tlay(0,k) = pt(k,j,i) * exner(k) 3160 3135 ENDDO 3161 3136 … … 3202 3177 rrtm_icld = 0 3203 3178 3204 IF ( cloud_physics.OR. cloud_droplets ) THEN3179 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 3205 3180 DO k = nzb+1, nzt+1 3206 3181 rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * & … … 3214 3189 ! 3215 3190 !-- Calculate cloud droplet effective radius 3216 IF ( cloud_physics) THEN3191 IF ( bulk_cloud_model ) THEN 3217 3192 ! 3218 3193 !-- Calculete effective droplet radius. In case of using … … 3285 3260 surf_lsm_h%frac(ind_wat_win,m) * & 3286 3261 surf_lsm_h%emissivity(ind_wat_win,m) 3287 rrtm_tsfc = surf_lsm_h%pt_surface(m) * & 3288 (surface_pressure / 1000.0_wp )**0.286_wp 3262 rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb) 3289 3263 ENDDO 3290 3264 DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i) … … 3295 3269 surf_usm_h%frac(ind_wat_win,m) * & 3296 3270 surf_usm_h%emissivity(ind_wat_win,m) 3297 rrtm_tsfc = surf_usm_h%pt_surface(m) * & 3298 (surface_pressure / 1000.0_wp )**0.286_wp 3271 rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb) 3299 3272 ENDDO 3300 3273 ! … … 3985 3958 ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) ) 3986 3959 3987 t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp3960 t_surface = pt_surface * exner(nzb) 3988 3961 DO k = nzb+1, nzt+1 3989 3962 rrtm_play(0,k) = hyp(k) * 0.01_wp 3990 rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / & 3991 t_surface )**(1.0_wp/0.286_wp) 3963 rrtm_plev(0,k) = barometric_formula(zw(k-1), & 3964 pt_surface * exner(nzb), & 3965 surface_pressure ) 3992 3966 ENDDO 3993 3967 … … 4406 4380 SUBROUTINE radiation_tendency_ij ( i, j, tend ) 4407 4381 4408 USE cloud_parameters, &4409 ONLY: pt_d_t4410 4411 4382 IMPLICIT NONE 4412 4383 … … 4421 4392 DO k = nzb+1, nzt+1 4422 4393 tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) & 4423 * pt_d_t(k) * d_seconds_hour4394 * d_exner(k) * d_seconds_hour 4424 4395 ENDDO 4425 4396 #endif … … 4437 4408 SUBROUTINE radiation_tendency ( tend ) 4438 4409 4439 USE cloud_parameters, &4440 ONLY: pt_d_t4441 4442 4410 USE indices, & 4443 4411 ONLY: nxl, nxr, nyn, nys … … 4457 4425 DO k = nzb+1, nzt+1 4458 4426 tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) & 4459 + rad_sw_hr(k,j,i) ) * pt_d_t(k)&4427 + rad_sw_hr(k,j,i) ) * d_exner(k) & 4460 4428 * d_seconds_hour 4461 4429 ENDDO … … 4525 4493 #if ! defined( __nopointer ) 4526 4494 IF ( plant_canopy ) THEN 4527 pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp&4528 / (c p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)4529 pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp&4495 pchf_prep(:) = r_d * exner(nzub:nzut) & 4496 / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T) 4497 pctf_prep(:) = r_d * exner(nzub:nzut) & 4530 4498 / (l_v * hyp(nzub:nzut) * dx*dy*dz(1)) 4531 4499 ENDIF … … 5032 5000 ! k = surf_lsm_h%k(m) 5033 5001 ! pt_surf_urb_l = pt_surf_urb_l + surf_lsm_h%pt_surface(m) & 5034 ! * ( hyp(k) / 100000.0_wp )**0.286_wp5002 ! * exner(k) 5035 5003 ! count_surfaces_l = count_surfaces_l + 1.0_wp 5036 5004 ! ENDDO … … 5038 5006 ! k = surf_usm_h%k(m) 5039 5007 ! pt_surf_urb_l = pt_surf_urb_l + surf_usm_h%pt_surface(m) & 5040 ! * ( hyp(k) / 100000.0_wp )**0.286_wp5008 ! * exner(k) 5041 5009 ! count_surfaces_l = count_surfaces_l + 1.0_wp 5042 5010 ! ENDDO … … 5045 5013 ! k = surf_lsm_v(l)%k(m) 5046 5014 ! pt_surf_urb_l = pt_surf_urb_l + surf_lsm_v(l)%pt_surface(m) & 5047 ! * ( hyp(k) / 100000.0_wp )**0.286_wp5015 ! * exner(k) 5048 5016 ! count_surfaces_l = count_surfaces_l + 1.0_wp 5049 5017 ! ENDDO 5050 5018 ! DO m = 1, surf_usm_v(l)%ns 5051 5019 ! k = surf_usm_v(l)%k(m) 5052 ! pt_surf_urb_l = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m) & 5053 ! * ( hyp(k) / 100000.0_wp )**0.286_wp 5020 ! pt_surf_urb_l = pt_surf_urb_l + surf_usm_v(l)%pt_surface(m) * exner(k) 5054 5021 ! count_surfaces_l = count_surfaces_l + 1.0_wp 5055 5022 ! ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
