Changeset 3014 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- May 9, 2018 8:42:38 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (59 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r3004 r3014 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Introduced plant canopy height similar to urban canopy height to limit 31 ! the memory requirement to allocate lad. 32 ! Deactivated automatic setting of minimum raytracing distance. 33 ! 34 ! 3004 2018-04-27 12:33:25Z Giersch 30 35 ! Further allocation checks implemented (averaged data will be assigned to fill 31 36 ! values if no allocation happened so far) … … 329 334 330 335 USE cloud_parameters, & 331 ONLY: cp, l_d_cp, r_d, rho_l336 ONLY: cp, l_d_cp, l_v, r_d, rho_l 332 337 333 338 USE constants, & … … 372 377 373 378 USE plant_canopy_model_mod, & 374 ONLY: pc_heating_rate, lad_s379 ONLY: lad_s, pc_heating_rate, pc_transpiration_rate 375 380 376 381 USE pegrid … … 646 651 !-- Parameters of urban and land surface models 647 652 INTEGER(iwp) :: nzu !< number of layers of urban surface (will be calculated) 653 INTEGER(iwp) :: nzp !< number of layers of plant canopy (will be calculated) 648 654 INTEGER(iwp) :: nzub,nzut !< bottom and top layer of urban surface (will be calculated) 655 INTEGER(iwp) :: nzpt !< top layer of plant canopy (will be calculated) 649 656 !-- parameters of urban and land surface models 650 657 INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers above top of of topography … … 4324 4331 REAL(wp), DIMENSION(3) :: sunorig_grid !< grid squashed solar direction unit vector (zyx) 4325 4332 REAL(wp), DIMENSION(0:nsurf_type) :: costheta !< direct irradiance factor of solar angle 4326 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temp tendency 4333 REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep !< precalculated factor for canopy temperature tendency 4334 REAL(wp), DIMENSION(nzub:nzut) :: pctf_prep !< precalculated factor for canopy transpiration tendency 4327 4335 REAL(wp), PARAMETER :: alpha = 0._wp !< grid rotation (TODO: add to namelist or remove) 4328 4336 REAL(wp) :: pc_box_area, pc_abs_frac, pc_abs_eff … … 4345 4353 REAL(wp) :: area_surf !< total area of surfaces in all processor 4346 4354 4355 4356 4347 4357 #if ! defined( __nopointer ) 4348 4358 IF ( plant_canopy ) THEN 4349 4359 pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 4350 4360 / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T) 4361 pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp & 4362 / (l_v * hyp(nzub:nzut) * dx*dy*dz) 4351 4363 ENDIF 4352 4364 #endif … … 4619 4631 !-- push heat flux absorbed by plant canopy to respective 3D arrays 4620 4632 IF ( npcbl > 0 ) THEN 4621 pc_heating_rate(:,:,:) = 0._wp 4633 pc_heating_rate(:,:,:) = 0.0_wp 4634 pc_transpiration_rate(:,:,:) = 0.0_wp 4622 4635 DO ipcgb = 1, npcbl 4623 4636 … … 4630 4643 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 4631 4644 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt 4645 4646 ! pc_transpiration_rate(kk,j,i) = 0.75_wp* (pcbinsw(ipcgb) + pcbinlw(ipcgb)) & 4647 ! * pctf_prep(k) * pt(k, j, i) !-- = dq/dt 4648 4632 4649 ENDDO 4633 4650 ENDIF … … 4995 5012 INTEGER(iwp) :: k_topo !< vertical index indicating topography top for given (j,i) 4996 5013 INTEGER(iwp) :: k_topo2 !< vertical index indicating topography top for given (j,i) 4997 INTEGER(iwp) :: nz ubl, nzutl, isurf, ipcgb5014 INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb 4998 5015 INTEGER(iwp) :: procid 4999 5016 REAL(wp) :: mrl … … 5040 5057 5041 5058 nzutl = MAX( nzutl, MAXVAL( pct ) ) 5059 nzptl = MAXVAL( pct ) 5042 5060 !-- code of plant canopy model uses parameter pch_index 5043 5061 !-- we need to setup it here to right value … … 5058 5076 CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr ) 5059 5077 CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5078 CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr ) 5060 5079 #else 5061 5080 nzub = nzubl 5062 5081 nzut = nzutl 5082 nzpt = nzptl 5063 5083 #endif 5064 5084 ! 5065 !-- global number of urban layers5085 !-- global number of urban and plant layers 5066 5086 nzu = nzut - nzub + 1 5087 nzp = nzpt - nzub + 1 5067 5088 ! 5068 5089 !-- check max_raytracing_dist relative to urban surface layer height 5069 mrl = 2.0_wp * nzu * dz5070 IF ( max_raytracing_dist <= mrl ) THEN5071 IF ( max_raytracing_dist /= -999.0_wp ) THEN5072 ! -- max_raytracing_dist too low5073 WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist too low, ' &5074 // 'override to value ', mrl5075 CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0)5076 ENDIF5077 max_raytracing_dist = mrl5078 ENDIF5090 ! mrl = 2.0_wp * nzu * dz 5091 ! IF ( max_raytracing_dist <= mrl ) THEN 5092 ! IF ( max_raytracing_dist /= -999.0_wp ) THEN 5093 ! !-- max_raytracing_dist too low 5094 ! WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist too low, ' & 5095 ! // 'override to value ', mrl 5096 ! CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0) 5097 ! ENDIF 5098 ! max_raytracing_dist = mrl 5099 ! ENDIF 5079 5100 ! 5080 5101 !-- allocate urban surfaces grid … … 5103 5124 IF ( npcbl > 0 ) THEN 5104 5125 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 5105 ALLOCATE( gridpcbl(nzub:nz ut,nys:nyn,nxl:nxr) )5126 ALLOCATE( gridpcbl(nzub:nzpt,nys:nyn,nxl:nxr) ) 5106 5127 pcbl = -1 5107 5128 gridpcbl(:,:,:) = 0 … … 5388 5409 IF ( plant_canopy ) THEN 5389 5410 ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) ) 5390 maxboxesg = nx + ny + nz u+ 15411 maxboxesg = nx + ny + nzp + 1 5391 5412 max_track_len = nx + ny + 1 5392 5413 !-- temporary arrays storing values for csf calculation during raytracing … … 5416 5437 !-- optimization of memory should be done 5417 5438 !-- Argument X of function STORAGE_SIZE(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now 5418 size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nz u5439 size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nzp 5419 5440 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, & 5420 5441 lad_s_rma_p, win_lad, ierr) 5421 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nz u, nny, nnx /))5442 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp, nny, nnx /)) 5422 5443 sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:) 5423 5444 ELSE 5424 ALLOCATE(sub_lad(nzub:nz ut, nys:nyn, nxl:nxr))5445 ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr)) 5425 5446 ENDIF 5426 5447 #else 5427 5448 plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) ) 5428 ALLOCATE(sub_lad(nzub:nz ut, nys:nyn, nxl:nxr))5449 ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr)) 5429 5450 #endif 5430 5451 plantt_max = MAXVAL(plantt) … … 5437 5458 k = get_topography_top_index_ji( j, i, 's' ) 5438 5459 5439 sub_lad(k:nz ut, j, i) = lad_s(0:nzut-k, j, i)5460 sub_lad(k:nzpt, j, i) = lad_s(0:nzpt-k, j, i) 5440 5461 ENDDO 5441 5462 ENDDO … … 5446 5467 CALL MPI_Win_lock_all(0, win_lad, ierr) 5447 5468 ELSE 5448 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nz u-1) )5449 CALL MPI_AllGather( sub_lad, nnx*nny*nz u, MPI_REAL, &5450 sub_lad_g, nnx*nny*nz u, MPI_REAL, comm2d, ierr )5469 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) ) 5470 CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, & 5471 sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr ) 5451 5472 ENDIF 5452 5473 #endif … … 6142 6163 #if defined( __parallel ) 6143 6164 lad_ip(ncsb) = ip 6144 lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nz u) + (box(2)-py*nny)*nzu+ box(1)-nzub6165 lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzp) + (box(2)-py*nny)*nzp + box(1)-nzub 6145 6166 #endif 6146 6167 ENDIF … … 6185 6206 lad_s_target = lad_s_ray(i) 6186 6207 ELSE 6187 lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nz u+ lad_disp(i))6208 lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nzp + lad_disp(i)) 6188 6209 ENDIF 6189 6210 #else … … 6387 6408 ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny 6388 6409 IF ( plantt(ig) <= nzterr(ig) ) CYCLE 6389 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nz u) + (rt2_track(1,i)-py*nny)*nzu+ nzterr(ig)+1-nzub6410 wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + nzterr(ig)+1-nzub 6390 6411 wcount = plantt(ig)-nzterr(ig) 6391 6412 ! TODO send request ASAP - even during raytracing … … 6411 6432 py = rt2_track(1,i)/nny 6412 6433 ip = px*pdims(2)+py 6413 ig = ip*nnx*nny*nz u + (rt2_track(2,i)-px*nnx)*(nny*nzu) + (rt2_track(1,i)-py*nny)*nzu6434 ig = ip*nnx*nny*nzp + (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp 6414 6435 rt2_track_lad(nzub:plantt_max, i) = sub_lad_g(ig:ig+nly-1) 6415 6436 ENDDO … … 6429 6450 !--Assert that we have space allocated for CSFs 6430 6451 !-- 6431 maxboxes = (ntrack + MAX(origin(1) - nzub, nz ut - origin(1))) * SIZE(zdirs, 1)6452 maxboxes = (ntrack + MAX(origin(1) - nzub, nzpt - origin(1))) * SIZE(zdirs, 1) 6432 6453 IF ( ncsfl + maxboxes > ncsfla ) THEN 6433 6454 !-- use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1) … … 6814 6835 ENDIF 6815 6836 ! 6816 !-- Close binary file 6837 !-- Close binary file 6817 6838 CALL close_file( fsvf ) 6818 6839 … … 7617 7638 !------------------------------------------------------------------------------! 7618 7639 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, & 7619 local_pf, two_d )7640 local_pf, two_d, nzb_do, nzt_do ) 7620 7641 7621 7642 USE indices … … 7635 7656 INTEGER(iwp) :: k !< 7636 7657 INTEGER(iwp) :: m !< index of surface element at grid point (j,i) 7658 INTEGER(iwp) :: nzb_do !< 7659 INTEGER(iwp) :: nzt_do !< 7637 7660 7638 7661 LOGICAL :: found !< … … 7641 7664 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 7642 7665 7643 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb :nzt+1) :: local_pf !<7666 REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 7644 7667 7645 7668 found = .TRUE. … … 7685 7708 DO i = nxl, nxr 7686 7709 DO j = nys, nyn 7687 DO k = nzb , nzt+17710 DO k = nzb_do, nzt_do 7688 7711 local_pf(i,j,k) = rad_lw_in(k,j,i) 7689 7712 ENDDO … … 7697 7720 DO i = nxl, nxr 7698 7721 DO j = nys, nyn 7699 DO k = nzb , nzt+17722 DO k = nzb_do, nzt_do 7700 7723 local_pf(i,j,k) = rad_lw_in_av(k,j,i) 7701 7724 ENDDO … … 7709 7732 DO i = nxl, nxr 7710 7733 DO j = nys, nyn 7711 DO k = nzb , nzt+17734 DO k = nzb_do, nzt_do 7712 7735 local_pf(i,j,k) = rad_lw_out(k,j,i) 7713 7736 ENDDO … … 7721 7744 DO i = nxl, nxr 7722 7745 DO j = nys, nyn 7723 DO k = nzb , nzt+17746 DO k = nzb_do, nzt_do 7724 7747 local_pf(i,j,k) = rad_lw_out_av(k,j,i) 7725 7748 ENDDO … … 7733 7756 DO i = nxl, nxr 7734 7757 DO j = nys, nyn 7735 DO k = nzb , nzt+17758 DO k = nzb_do, nzt_do 7736 7759 local_pf(i,j,k) = rad_lw_cs_hr(k,j,i) 7737 7760 ENDDO … … 7745 7768 DO i = nxl, nxr 7746 7769 DO j = nys, nyn 7747 DO k = nzb , nzt+17770 DO k = nzb_do, nzt_do 7748 7771 local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i) 7749 7772 ENDDO … … 7757 7780 DO i = nxl, nxr 7758 7781 DO j = nys, nyn 7759 DO k = nzb , nzt+17782 DO k = nzb_do, nzt_do 7760 7783 local_pf(i,j,k) = rad_lw_hr(k,j,i) 7761 7784 ENDDO … … 7769 7792 DO i = nxl, nxr 7770 7793 DO j = nys, nyn 7771 DO k = nzb , nzt+17794 DO k = nzb_do, nzt_do 7772 7795 local_pf(i,j,k) = rad_lw_hr_av(k,j,i) 7773 7796 ENDDO … … 7781 7804 DO i = nxl, nxr 7782 7805 DO j = nys, nyn 7783 DO k = nzb , nzt+17806 DO k = nzb_do, nzt_do 7784 7807 local_pf(i,j,k) = rad_sw_in(k,j,i) 7785 7808 ENDDO … … 7793 7816 DO i = nxl, nxr 7794 7817 DO j = nys, nyn 7795 DO k = nzb , nzt+17818 DO k = nzb_do, nzt_do 7796 7819 local_pf(i,j,k) = rad_sw_in_av(k,j,i) 7797 7820 ENDDO … … 7805 7828 DO i = nxl, nxr 7806 7829 DO j = nys, nyn 7807 DO k = nzb , nzt+17830 DO k = nzb_do, nzt_do 7808 7831 local_pf(i,j,k) = rad_sw_out(k,j,i) 7809 7832 ENDDO … … 7829 7852 DO i = nxl, nxr 7830 7853 DO j = nys, nyn 7831 DO k = nzb , nzt+17854 DO k = nzb_do, nzt_do 7832 7855 local_pf(i,j,k) = rad_sw_cs_hr(k,j,i) 7833 7856 ENDDO … … 7841 7864 DO i = nxl, nxr 7842 7865 DO j = nys, nyn 7843 DO k = nzb , nzt+17866 DO k = nzb_do, nzt_do 7844 7867 local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i) 7845 7868 ENDDO … … 7853 7876 DO i = nxl, nxr 7854 7877 DO j = nys, nyn 7855 DO k = nzb , nzt+17878 DO k = nzb_do, nzt_do 7856 7879 local_pf(i,j,k) = rad_sw_hr(k,j,i) 7857 7880 ENDDO … … 7865 7888 DO i = nxl, nxr 7866 7889 DO j = nys, nyn 7867 DO k = nzb , nzt+17890 DO k = nzb_do, nzt_do 7868 7891 local_pf(i,j,k) = rad_sw_hr_av(k,j,i) 7869 7892 ENDDO … … 7888 7911 !> Subroutine defining 3D output variables 7889 7912 !------------------------------------------------------------------------------! 7890 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )7913 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 7891 7914 7892 7915 … … 7904 7927 INTEGER(iwp) :: j !< 7905 7928 INTEGER(iwp) :: k !< 7929 INTEGER(iwp) :: nzb_do !< 7930 INTEGER(iwp) :: nzt_do !< 7906 7931 7907 7932 LOGICAL :: found !< … … 7909 7934 REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute 7910 7935 7911 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb :nzt+1) :: local_pf !<7936 REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< 7912 7937 7913 7938 … … 7921 7946 DO i = nxl, nxr 7922 7947 DO j = nys, nyn 7923 DO k = nzb , nzt+17948 DO k = nzb_do, nzt_do 7924 7949 local_pf(i,j,k) = rad_sw_in(k,j,i) 7925 7950 ENDDO … … 7933 7958 DO i = nxl, nxr 7934 7959 DO j = nys, nyn 7935 DO k = nzb , nzt+17960 DO k = nzb_do, nzt_do 7936 7961 local_pf(i,j,k) = rad_sw_in_av(k,j,i) 7937 7962 ENDDO … … 7944 7969 DO i = nxl, nxr 7945 7970 DO j = nys, nyn 7946 DO k = nzb , nzt+17971 DO k = nzb_do, nzt_do 7947 7972 local_pf(i,j,k) = rad_sw_out(k,j,i) 7948 7973 ENDDO … … 7956 7981 DO i = nxl, nxr 7957 7982 DO j = nys, nyn 7958 DO k = nzb , nzt+17983 DO k = nzb_do, nzt_do 7959 7984 local_pf(i,j,k) = rad_sw_out_av(k,j,i) 7960 7985 ENDDO … … 7967 7992 DO i = nxl, nxr 7968 7993 DO j = nys, nyn 7969 DO k = nzb , nzt+17994 DO k = nzb_do, nzt_do 7970 7995 local_pf(i,j,k) = rad_sw_cs_hr(k,j,i) 7971 7996 ENDDO … … 7979 8004 DO i = nxl, nxr 7980 8005 DO j = nys, nyn 7981 DO k = nzb , nzt+18006 DO k = nzb_do, nzt_do 7982 8007 local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i) 7983 8008 ENDDO … … 7990 8015 DO i = nxl, nxr 7991 8016 DO j = nys, nyn 7992 DO k = nzb , nzt+18017 DO k = nzb_do, nzt_do 7993 8018 local_pf(i,j,k) = rad_sw_hr(k,j,i) 7994 8019 ENDDO … … 8002 8027 DO i = nxl, nxr 8003 8028 DO j = nys, nyn 8004 DO k = nzb , nzt+18029 DO k = nzb_do, nzt_do 8005 8030 local_pf(i,j,k) = rad_sw_hr_av(k,j,i) 8006 8031 ENDDO … … 8013 8038 DO i = nxl, nxr 8014 8039 DO j = nys, nyn 8015 DO k = nzb , nzt+18040 DO k = nzb_do, nzt_do 8016 8041 local_pf(i,j,k) = rad_lw_in(k,j,i) 8017 8042 ENDDO … … 8025 8050 DO i = nxl, nxr 8026 8051 DO j = nys, nyn 8027 DO k = nzb , nzt+18052 DO k = nzb_do, nzt_do 8028 8053 local_pf(i,j,k) = rad_lw_in_av(k,j,i) 8029 8054 ENDDO … … 8036 8061 DO i = nxl, nxr 8037 8062 DO j = nys, nyn 8038 DO k = nzb , nzt+18063 DO k = nzb_do, nzt_do 8039 8064 local_pf(i,j,k) = rad_lw_out(k,j,i) 8040 8065 ENDDO … … 8048 8073 DO i = nxl, nxr 8049 8074 DO j = nys, nyn 8050 DO k = nzb , nzt+18075 DO k = nzb_do, nzt_do 8051 8076 local_pf(i,j,k) = rad_lw_out_av(k,j,i) 8052 8077 ENDDO … … 8059 8084 DO i = nxl, nxr 8060 8085 DO j = nys, nyn 8061 DO k = nzb , nzt+18086 DO k = nzb_do, nzt_do 8062 8087 local_pf(i,j,k) = rad_lw_cs_hr(k,j,i) 8063 8088 ENDDO … … 8071 8096 DO i = nxl, nxr 8072 8097 DO j = nys, nyn 8073 DO k = nzb , nzt+18098 DO k = nzb_do, nzt_do 8074 8099 local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i) 8075 8100 ENDDO … … 8082 8107 DO i = nxl, nxr 8083 8108 DO j = nys, nyn 8084 DO k = nzb , nzt+18109 DO k = nzb_do, nzt_do 8085 8110 local_pf(i,j,k) = rad_lw_hr(k,j,i) 8086 8111 ENDDO … … 8094 8119 DO i = nxl, nxr 8095 8120 DO j = nys, nyn 8096 DO k = nzb , nzt+18121 DO k = nzb_do, nzt_do 8097 8122 local_pf(i,j,k) = rad_lw_hr_av(k,j,i) 8098 8123 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
