 Timestamp:
 Mar 29, 2019 12:05:13 PM (6 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/radiation_model_mod.f90
r3826 r3843 23 23 ! Current revisions: 24 24 !  25 ! 25 ! Code review. 26 26 ! 27 27 ! Former revisions: … … 761 761 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_xy_av !< average of incoming shortwave radiation at surface 762 762 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface 763 764 REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp !< emissivity of the clearsky atmosphere 763 765 ! 764 766 ! Land surface albedos for solar zenith angle of 60degree after Briegleb (1992) … … 934 936 INTEGER(iwp), PARAMETER :: iy = 3 !< position of jindex in surfl and surf 935 937 INTEGER(iwp), PARAMETER :: ix = 4 !< position of iindex in surfl and surf 938 INTEGER(iwp), PARAMETER :: im = 5 !< position of surface mindex in surfl and surf 939 INTEGER(iwp), PARAMETER :: nidx_surf = 5 !< number of indices in surfl and surf 936 940 937 941 INTEGER(iwp), PARAMETER :: nsurf_type = 10 !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces  1 … … 974 978 975 979 ! indices and sizes of urban and land surface models 976 INTEGER(iwp), DIMENSION(: ), ALLOCATABLE,TARGET :: surfl_l !< coordinates of ith local surface in local grid  surfl[:,k] = [d, z, y, x]977 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surf _l !< coordinates of ith surface in grid  surf[:,k] = [d, z, y, x]978 INTEGER(iwp), DIMENSION(:,:), POINTER :: surf l !< coordinates of ith local surface in local grid  surfl[:,k] = [d, z, y, x]979 INTEGER(iwp), DIMENSION(: ,:), POINTER :: surf !< coordinates of ith surface in grid  surf[:,k] = [d, z, y, x]980 INTEGER(iwp), DIMENSION(:,:), POINTER :: surfl !< coordinates of ith local surface in local grid  surfl[:,k] = [d, z, y, x, m] 981 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfl_linear !< dtto (linearly allocated array) 982 INTEGER(iwp), DIMENSION(:,:), POINTER :: surf !< coordinates of ith surface in grid  surf[:,k] = [d, z, y, x, m] 983 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surf_linear !< dtto (linearly allocated array) 980 984 INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor 981 985 INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: nsurfs !< array of number of all surfaces in individual processors … … 2557 2561 ! Calculate initial values of current (cosine of) the zenith angle and 2558 2562 ! whether the sun is up 2559 CALL calc_zenith 2560 ! readjust date and time to its initial value 2563 CALL calc_zenith 2564 ! 2565 ! readjust date and time to its initial value 2561 2566 CALL init_date_and_time 2562 2567 ! … … 2566 2571 ! 2567 2572 ! Horizontally aligned natural and urban surfaces 2568 CALL calc_albedo( surf_lsm_h 2569 CALL calc_albedo( surf_usm_h 2573 CALL calc_albedo( surf_lsm_h ) 2574 CALL calc_albedo( surf_usm_h ) 2570 2575 ! 2571 2576 ! Vertically aligned natural and urban surfaces … … 2652 2657 IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN 2653 2658 ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2654 rad_lw_in 2659 rad_lw_in = 0.0_wp 2655 2660 ENDIF 2656 2661 … … 2661 2666 IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN 2662 2667 ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2663 rad_lw_out 2668 rad_lw_out = 0.0_wp 2664 2669 ENDIF 2665 2670 … … 2917 2922 surf%rad_sw_out = albedo_urb * surf%rad_sw_in 2918 2923 2919 surf%rad_lw_in = 0.8_wp* sigma_sb * (pt1 * exner(k+1))**42924 surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4 2920 2925 2921 2926 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 2976 2981 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 2977 2982 pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i) 2978 surf%rad_lw_in(m) = 0.8_wp* sigma_sb * (pt1 * exner(k))**42983 surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4 2979 2984 ELSE 2980 surf%rad_lw_in(m) = 0.8_wp* sigma_sb * (pt(k,j,i) * exner(k))**42985 surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt(k,j,i) * exner(k))**4 2981 2986 ENDIF 2982 2987 … … 3093 3098 surf%rad_net = net_radiation 3094 3099 3095 surf%rad_lw_in = 0.8_wp* sigma_sb * (pt1 * exner(nz_urban_t+1))**43100 surf%rad_lw_in = emissivity_atm_clsky * sigma_sb * (pt1 * exner(nz_urban_t+1))**4 3096 3101 3097 3102 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 & … … 3129 3134 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 3130 3135 pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i) 3131 surf%rad_lw_in(m) = 0.8_wp* sigma_sb * (pt1 * exner(k))**43136 surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4 3132 3137 ELSE 3133 surf%rad_lw_in(m) = 0.8_wp* sigma_sb * &3138 surf%rad_lw_in(m) = emissivity_atm_clsky * sigma_sb * & 3134 3139 ( pt(k,j,i) * exner(k) )**4 3135 3140 ENDIF … … 4988 4993 REAL(wp), DIMENSION(0:nsurf_type) :: costheta !< direct irradiance factor of solar angle 4989 4994 REAL(wp), DIMENSION(nz_urban_b:nz_urban_t) :: pchf_prep !< precalculated factor for canopy temperature tendency 4990 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4995 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: synchronize with rotation_angle 4996 !< from netcdf_data_input_mod) 4991 4997 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff 4992 4998 REAL(wp) :: asrc !< area of source face … … 5437 5443 ! Transfer radiation arrays required for energy balance to the respective data types 5438 5444 DO i = 1, nsurfl 5439 m = surfl( 5,i)5445 m = surfl(im,i) 5440 5446 ! 5441 5447 ! (1) Urban surfaces … … 6096 6102 ENDIF 6097 6103 6098 ! fill surfl (the ordering of local surfaces given by the following 6104 ! 6105 ! Fill surfl (the ordering of local surfaces given by the following 6099 6106 ! cycles must not be altered, certain file input routines may depend 6100 ! on it) 6101 ALLOCATE(surfl_l(5*nsurfl)) ! is it necessary to allocate it with (5,nsurfl)? 6102 surfl(1:5,1:nsurfl) => surfl_l(1:5*nsurfl) 6107 ! on it). 6108 ! 6109 ! We allocate the array as linear and then use a twodimensional pointer 6110 ! into it, because some MPI implementations crash with 2Dallocated arrays. 6111 ALLOCATE(surfl_linear(nidx_surf*nsurfl)) 6112 surfl(1:nidx_surf,1:nsurfl) => surfl_linear(1:nidx_surf*nsurfl) 6103 6113 isurf = 0 6104 6114 IF ( rad_angular_discretization ) THEN … … 6356 6366 surfstart(numprocs) = k 6357 6367 nsurf = k 6358 ALLOCATE(surf_l(5*nsurf)) 6359 surf(1:5,1:nsurf) => surf_l(1:5*nsurf) 6368 ! 6369 ! We allocate the array as linear and then use a twodimensional pointer 6370 ! into it, because some MPI implementations crash with 2Dallocated arrays. 6371 ALLOCATE(surf_linear(nidx_surf*nsurf)) 6372 surf(1:nidx_surf,1:nsurf) => surf_linear(1:nidx_surf*nsurf) 6360 6373 6361 6374 #if defined( __parallel ) 6362 CALL MPI_AllGatherv(surfl_l, nsurfl*5, MPI_INTEGER, surf_l, nsurfs*5, & 6363 surfstart(0:numprocs1)*5, MPI_INTEGER, comm2d, ierr) 6375 CALL MPI_AllGatherv(surfl_linear, nsurfl*nidx_surf, MPI_INTEGER, & 6376 surf_linear, nsurfs*nidx_surf, & 6377 surfstart(0:numprocs1)*nidx_surf, MPI_INTEGER, & 6378 comm2d, ierr) 6364 6379 IF ( ierr /= 0 ) THEN 6365 WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_l), nsurfl*5, & 6366 SIZE(surf_l), nsurfs*5, surfstart(0:numprocs1)*5 6380 WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_linear), & 6381 nsurfl*nidx_surf, SIZE(surf_linear), nsurfs*nidx_surf, & 6382 surfstart(0:numprocs1)*nidx_surf 6367 6383 FLUSH(9) 6368 6384 ENDIF … … 7378 7394 icsf = 1 !< reading index 7379 7395 kcsf = 1 !< writing index 7380 DO while(icsf < npcsfl)7396 DO WHILE (icsf < npcsfl) 7381 7397 ! here kpcsf(kcsf) already has values from kpcsf(icsf) 7382 7398 IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1) .AND. & … … 8670 8686 8671 8687 8688 !! 8689 ! 8690 ! Description: 8691 !  8692 !> Grows the CSF array exponentially after it is full. During that, the ray 8693 !> canopy sink factors with common source face and target plant canopy grid 8694 !> cell are merged together so that the size doesn't grow out of control. 8695 !! 8672 8696 SUBROUTINE merge_and_grow_csf(newsize) 8673 8697 INTEGER(iwp), INTENT(in) :: newsize !< new array size after grow, must be >= ncsfl
Note: See TracChangeset
for help on using the changeset viewer.