Changeset 2977 for palm/trunk/SOURCE/radiation_model_mod.f90
- Timestamp:
- Apr 17, 2018 10:27:57 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/radiation_model_mod.f90
r2967 r2977 22 22 ! 23 23 ! Current revisions: 24 ! ----------------- 24 ! ------------------ 25 25 ! 26 26 ! … … 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Implement changes from branch radiation (r2948-2971) with minor modifications, 31 ! plus some formatting. 32 ! (moh.hefny): 33 ! - replaced plant_canopy by npcbl to check tree existence to avoid weird 34 ! allocation of related arrays (after domain decomposition some domains 35 ! contains no trees although plant_canopy (global parameter) is still TRUE). 36 ! - added a namelist parameter to force RTM settings 37 ! - enabled the option to switch radiation reflections off 38 ! - renamed surf_reflections to surface_reflections 39 ! - removed average_radiation flag from the namelist (now it is implicitly set 40 ! in init_3d_model according to RTM) 41 ! - edited read and write sky view factors and CSF routines to account for 42 ! the sub-domains which may not contain any of them 43 ! 44 ! 2967 2018-04-13 11:22:08Z raasch 30 45 ! bugfix: missing parallel cpp-directives added 31 46 ! … … 432 447 sun_direction = .FALSE., & !< flag parameter indicating whether solar direction shall be calculated 433 448 average_radiation = .FALSE., & !< flag to set the calculation of radiation averaging for the domain 434 radiation_interactions = .FALSE., & !< flag to control if radiation interactions via sky-view factors shall be considered435 surf _reflections = .TRUE.!< flag to switch the calculation of radiation interaction between surfaces.449 radiation_interactions = .FALSE., & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist) 450 surface_reflections = .TRUE., & !< flag to switch the calculation of radiation interaction between surfaces. 436 451 !< When it switched off, only the effect of buildings and trees shadow will 437 452 !< will be considered. However fewer SVFs are expected. 453 radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees 438 454 439 455 … … 680 696 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy 681 697 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy 682 INTEGER(iwp) :: npcbl 698 INTEGER(iwp) :: npcbl = 0 !< number of the plant canopy gridboxes in local processor 683 699 INTEGER(wp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i] 684 700 REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box … … 955 971 surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart, & 956 972 surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw, & 957 pcbl, npcbl,iup_u, inorth_u, isouth_u, ieast_u, iwest_u, &973 iup_u, inorth_u, isouth_u, ieast_u, iwest_u, & 958 974 iup_l, inorth_l, isouth_l, ieast_l, iwest_l, & 959 975 nsurf_type, nzub, nzut, nzu, pch, nsurf, & … … 961 977 idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, & 962 978 radiation_interactions, startwall, startland, endland, endwall, & 963 skyvf, skyvft 979 skyvf, skyvft, radiation_interactions_on, average_radiation 964 980 965 981 #if defined ( __rrtmg ) … … 1309 1325 CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 ) 1310 1326 ENDIF 1311 ENDIF1312 1313 !1314 !-- Radiation interactions1315 IF ( nrefsteps < 1 .AND. radiation_interactions ) THEN1316 message_string = 'nrefsteps must be > 0 when using LSM/USM to' // &1317 'account for surface outgoing SW flux.' // &1318 'You may set surf_reflections = .FALSE. to ' // &1319 'diable surface reflections instead.'1320 CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )1321 1327 ENDIF 1322 1328 … … 2634 2640 nrefsteps, mrt_factors, rma_lad_raytrace, & 2635 2641 dist_max_svf, & 2636 average_radiation,&2637 surf_reflections, svfnorm_report_thresh2642 surface_reflections, svfnorm_report_thresh, & 2643 radiation_interactions_on 2638 2644 2639 2645 NAMELIST /radiation_parameters/ albedo, albedo_type, albedo_lw_dir, & … … 2647 2653 nrefsteps, mrt_factors, rma_lad_raytrace, & 2648 2654 dist_max_svf, & 2649 average_radiation,&2650 surf_reflections, svfnorm_report_thresh2655 surface_reflections, svfnorm_report_thresh, & 2656 radiation_interactions_on 2651 2657 2652 2658 line = ' ' … … 2693 2699 radiation = .TRUE. 2694 2700 2701 IF ( .NOT. radiation_interactions_on .AND. surface_reflections ) THEN 2702 message_string = 'surface_reflections is allowed only when ' // & 2703 'radiation_interactions_on is set to TRUE' 2704 CALL message( 'radiation_parin', 'PA0487',1, 2, 0, 6, 0 ) 2705 ENDIF 2706 2695 2707 12 CONTINUE 2696 2697 2698 2699 !-- Set radiation_interactions flag according to urban_ and land_surface flag2700 IF ( urban_surface .OR. land_surface ) radiation_interactions = .TRUE.2701 2708 2702 2709 END SUBROUTINE radiation_parin … … 4354 4361 sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2)) 4355 4362 4356 IF ( plant_canopy) THEN4363 IF ( npcbl > 0 ) THEN 4357 4364 !-- precompute effective box depth with prototype Leaf Area Density 4358 4365 pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1 … … 4465 4472 #endif 4466 4473 4467 DO isvf = 1, nsvfl 4468 isurf = svfsurf(1, isvf) 4469 k = surfl(iz, isurf) 4470 j = surfl(iy, isurf) 4471 i = surfl(ix, isurf) 4472 isurfsrc = svfsurf(2, isvf) 4473 4474 !-- for surface-to-surface factors we calculate thermal radiation in 1st pass 4475 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4476 ENDDO 4474 IF ( surface_reflections) THEN 4475 DO isvf = 1, nsvfl 4476 isurf = svfsurf(1, isvf) 4477 k = surfl(iz, isurf) 4478 j = surfl(iy, isurf) 4479 i = surfl(ix, isurf) 4480 isurfsrc = svfsurf(2, isvf) 4481 ! 4482 !-- For surface-to-surface factors we calculate thermal radiation in 1st pass 4483 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc) 4484 ENDDO 4485 ENDIF 4477 4486 4478 4487 !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors … … 4494 4503 ENDIF 4495 4504 4496 IF ( plant_canopy) THEN4505 IF ( npcbl > 0 ) THEN 4497 4506 4498 4507 pcbinswdir(:) = 0._wp … … 4536 4545 ! surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw 4537 4546 4547 IF ( .NOT. surface_reflections ) THEN 4548 ! 4549 !-- Set nrefsteps to 0 to disable reflections 4550 nrefsteps = 0 4551 surfoutsl = albedo_surf * surfins 4552 surfoutll = (1._wp - emiss_surf) * surfinl 4553 surfoutsw = surfoutsw + surfoutsl 4554 surfoutlw = surfoutlw + surfoutll 4555 ENDIF 4556 4538 4557 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 4539 4558 !-- Next passes - reflections … … 4585 4604 4586 4605 !-- push heat flux absorbed by plant canopy to respective 3D arrays 4587 IF ( plant_canopy) THEN4606 IF ( npcbl > 0 ) THEN 4588 4607 pc_heating_rate(:,:,:) = 0._wp 4589 4608 DO ipcgb = 1, npcbl … … 4733 4752 !-- domain when using average_radiation in the respective radiation model 4734 4753 4735 IF ( average_radiation ) THEN 4736 4737 !-- precalculate face areas for all face directions using normal vector 4738 DO d = 0, nsurf_type 4739 facearea(d) = 1._wp 4740 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 4741 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 4742 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 4743 ENDDO 4744 ! 4745 !-- absorbed/received SW & LW and emitted LW energy of all physical 4746 !-- surfaces (land and urban) in local processor 4747 pinswl = 0._wp 4748 pinlwl = 0._wp 4749 pabsswl = 0._wp 4750 pabslwl = 0._wp 4751 pemitlwl = 0._wp 4752 emiss_sum_surfl = 0._wp 4753 area_surfl = 0._wp 4754 DO i = 1, nsurfl 4755 d = surfl(id, i) 4756 !-- received SW & LW 4757 pinswl = pinswl + surfinsw(i) * facearea(d) 4758 pinlwl = pinlwl + surfinlw(i) * facearea(d) 4759 !-- absorbed SW & LW 4760 pabsswl = pabsswl + (1._wp - albedo_surf(i)) * & 4761 surfinsw(i) * facearea(d) 4762 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d) 4763 !-- emitted LW 4764 pemitlwl = pemitlwl + surfoutlw(i) * facearea(d) 4765 !-- emissivity and area sum 4766 emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d) 4767 area_surfl = area_surfl + facearea(d) 4768 END DO 4769 ! 4770 !-- add the absorbed SW energy by plant canopy 4771 IF ( plant_canopy ) THEN 4772 pabsswl = pabsswl + SUM(pcbinsw) 4773 pabslwl = pabslwl + SUM(pcbinlw) 4774 ENDIF 4775 ! 4776 !-- gather all rad flux energy in all processors 4754 !-- Precalculate face areas for all face directions using normal vector 4755 DO d = 0, nsurf_type 4756 facearea(d) = 1._wp 4757 IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx 4758 IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy 4759 IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz 4760 ENDDO 4761 ! 4762 !-- absorbed/received SW & LW and emitted LW energy of all physical 4763 !-- surfaces (land and urban) in local processor 4764 pinswl = 0._wp 4765 pinlwl = 0._wp 4766 pabsswl = 0._wp 4767 pabslwl = 0._wp 4768 pemitlwl = 0._wp 4769 emiss_sum_surfl = 0._wp 4770 area_surfl = 0._wp 4771 DO i = 1, nsurfl 4772 d = surfl(id, i) 4773 !-- received SW & LW 4774 pinswl = pinswl + surfinsw(i) * facearea(d) 4775 pinlwl = pinlwl + surfinlw(i) * facearea(d) 4776 !-- absorbed SW & LW 4777 pabsswl = pabsswl + (1._wp - albedo_surf(i)) * & 4778 surfinsw(i) * facearea(d) 4779 pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d) 4780 !-- emitted LW 4781 pemitlwl = pemitlwl + surfoutlw(i) * facearea(d) 4782 !-- emissivity and area sum 4783 emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d) 4784 area_surfl = area_surfl + facearea(d) 4785 END DO 4786 ! 4787 !-- add the absorbed SW energy by plant canopy 4788 IF ( npcbl > 0 ) THEN 4789 pabsswl = pabsswl + SUM(pcbinsw) 4790 pabslwl = pabslwl + SUM(pcbinlw) 4791 ENDIF 4792 ! 4793 !-- gather all rad flux energy in all processors 4777 4794 #if defined( __parallel ) 4778 4779 4780 4781 4782 4783 4784 4795 CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 4796 CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr) 4797 CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4798 CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4799 CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4800 CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4801 CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 4785 4802 #else 4786 4787 4788 4789 4790 4791 4792 4803 pinsw = pinswl 4804 pinlw = pinlwl 4805 pabssw = pabsswl 4806 pabslwl = pabslw 4807 pemitlwl = pemitlw 4808 emiss_sum_surf = emiss_sum_surfl 4809 area_surf = area_surfl 4793 4810 #endif 4794 4811 4795 !-- (1) albedo 4796 IF ( pinsw /= 0.0_wp ) albedo_urb = 1._wp - pabssw / pinsw 4797 4798 !-- (2) average emmsivity 4799 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 4800 4801 !-- (3) temperature 4802 t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp 4803 4804 ENDIF 4812 !-- (1) albedo 4813 IF ( pinsw /= 0.0_wp ) albedo_urb = 1._wp - pabssw / pinsw 4814 4815 !-- (2) average emmsivity 4816 IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf 4817 4818 !-- (3) temperature 4819 t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp 4820 4805 4821 4806 4822 CONTAINS … … 5071 5087 5072 5088 !-- fill gridpcbl and pcbl 5073 IF ( plant_canopy) THEN5089 IF ( npcbl > 0 ) THEN 5074 5090 ALLOCATE( pcbl(iz:ix, 1:npcbl) ) 5075 5091 ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) ) … … 5593 5609 5594 5610 DEALLOCATE ( zdirs, zbdry, vffrac, ztransp ) 5595 5596 DO isurfs = 1, nsurf 5597 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 5598 surfl(iz, isurflt), surfl(id, isurflt), & 5599 surf(ix, isurfs), surf(iy, isurfs), & 5600 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5601 CYCLE 5602 ENDIF 5611 ! 5612 !-- Following calculations only required for surface_reflections 5613 IF ( surface_reflections ) THEN 5614 5615 DO isurfs = 1, nsurf 5616 IF ( .NOT. surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), & 5617 surfl(iz, isurflt), surfl(id, isurflt), & 5618 surf(ix, isurfs), surf(iy, isurfs), & 5619 surf(iz, isurfs), surf(id, isurfs)) ) THEN 5620 CYCLE 5621 ENDIF 5603 5622 5604 sd = surf(id, isurfs)5605 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &5606 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &5607 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)5608 5609 !-- unit vector source -> target5610 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)5611 sqdist = SUM(uv(:)**2)5612 uv = uv / SQRT(sqdist)5613 5614 !-- reject raytracing above max distance5615 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN5616 ray_skip_maxdist = ray_skip_maxdist + 15617 CYCLE5618 ENDIF5619 5620 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area5621 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction5622 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction5623 / (pi * sqdist) & ! square of distance between centers5624 * facearea(sd)5625 5626 !-- reject raytracing for potentially too small view factor values5627 IF ( rirrf < min_irrf_value ) THEN5628 ray_skip_minval = ray_skip_minval + 15629 CYCLE5630 ENDIF5631 5632 !-- raytrace + process plant canopy sinks within5633 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &5634 visible, transparency, win_lad)5635 5636 IF ( .NOT. visible ) CYCLE5637 ! rsvf = rirrf * transparency5638 5639 !-- write to the svf array5640 nsvfl = nsvfl + 15641 !-- check dimmension of asvf array and enlarge it if needed5642 IF ( nsvfla < nsvfl ) THEN5643 k = nsvfla * 25644 IF ( msvf == 0 ) THEN5623 sd = surf(id, isurfs) 5624 sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), & 5625 REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), & 5626 REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /) 5627 5628 !-- unit vector source -> target 5629 uv = (/ (ta(1)-sa(1))*dz, (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /) 5630 sqdist = SUM(uv(:)**2) 5631 uv = uv / SQRT(sqdist) 5632 5633 !-- reject raytracing above max distance 5634 IF ( SQRT(sqdist) > max_raytracing_dist ) THEN 5635 ray_skip_maxdist = ray_skip_maxdist + 1 5636 CYCLE 5637 ENDIF 5638 5639 !-- irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area 5640 rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction 5641 * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) & ! cosine of target normal and reverse direction 5642 / (pi * sqdist) & ! square of distance between centers 5643 * facearea(sd) 5644 5645 !-- reject raytracing for potentially too small view factor values 5646 IF ( rirrf < min_irrf_value ) THEN 5647 ray_skip_minval = ray_skip_minval + 1 5648 CYCLE 5649 ENDIF 5650 5651 !-- raytrace + process plant canopy sinks within 5652 CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., & 5653 visible, transparency, win_lad) 5654 5655 IF ( .NOT. visible ) CYCLE 5656 ! rsvf = rirrf * transparency 5657 5658 !-- write to the svf array 5659 nsvfl = nsvfl + 1 5660 !-- check dimmension of asvf array and enlarge it if needed 5661 IF ( nsvfla < nsvfl ) THEN 5662 k = nsvfla * 2 5663 IF ( msvf == 0 ) THEN 5645 5664 msvf = 1 5646 5665 ALLOCATE( asvf1(k) ) … … 5648 5667 asvf1(1:nsvfla) = asvf2 5649 5668 DEALLOCATE( asvf2 ) 5650 ELSE5669 ELSE 5651 5670 msvf = 0 5652 5671 ALLOCATE( asvf2(k) ) … … 5654 5673 asvf2(1:nsvfla) = asvf1 5655 5674 DEALLOCATE( asvf1 ) 5656 ENDIF 5657 5658 ! WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 5659 ! CALL radiation_write_debug_log( msg ) 5660 5661 nsvfla = k 5662 ENDIF 5663 !-- write svf values into the array 5664 asvf(nsvfl)%isurflt = isurflt 5665 asvf(nsvfl)%isurfs = isurfs 5666 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency 5667 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor 5668 ENDDO 5675 ENDIF 5676 5677 ! WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k 5678 ! CALL radiation_write_debug_log( msg ) 5679 5680 nsvfla = k 5681 ENDIF 5682 !-- write svf values into the array 5683 asvf(nsvfl)%isurflt = isurflt 5684 asvf(nsvfl)%isurfs = isurfs 5685 asvf(nsvfl)%rsvf = rirrf !we postopne multiplication by transparency 5686 asvf(nsvfl)%rtransp = transparency !a.k.a. Direct Irradiance Factor 5687 ENDDO 5688 ENDIF 5669 5689 ENDDO 5670 5690 … … 6733 6753 !-- read nsvfl, ncsfl 6734 6754 READ ( fsvf ) nsvfl, ncsfl, nsurfl_from_file 6735 IF ( nsvfl < =0 .OR. ncsfl < 0 ) THEN6755 IF ( nsvfl < 0 .OR. ncsfl < 0 ) THEN 6736 6756 WRITE( message_string, * ) 'Wrong number of SVF or CSF' 6737 6757 CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 ) … … 6747 6767 CALL message( 'radiation_read_svf', 'PA0490', 1, 2, 0, 6, 0 ) 6748 6768 ENDIF 6749 6750 IF ( .NOT. ALLOCATED( skyvf ) ) ALLOCATE( skyvf(nsurfl) ) 6751 IF ( .NOT. ALLOCATED( skyvft ) ) ALLOCATE( skyvft(nsurfl) ) 6752 IF ( .NOT. ALLOCATED( svf ) ) ALLOCATE( svf(ndsvf,nsvfl) ) 6753 IF ( .NOT. ALLOCATED( svfsurf ) ) ALLOCATE( svfsurf(idsvf,nsvfl) ) 6754 READ(fsvf) skyvf 6755 READ(fsvf) skyvft 6756 READ(fsvf) svf 6757 READ(fsvf) svfsurf 6758 IF ( plant_canopy ) THEN 6759 IF ( .NOT. ALLOCATED( csf ) ) ALLOCATE( csf(ndcsf,ncsfl) ) 6760 IF ( .NOT. ALLOCATED( csfsurf ) ) ALLOCATE( csfsurf(idcsf,ncsfl) ) 6761 READ(fsvf) csf 6762 READ(fsvf) csfsurf 6769 IF ( nsurfl > 0 ) THEN 6770 IF ( .NOT. ALLOCATED( skyvf ) ) ALLOCATE( skyvf(nsurfl) ) 6771 IF ( .NOT. ALLOCATED( skyvft ) ) ALLOCATE( skyvft(nsurfl) ) 6772 READ(fsvf) skyvf 6773 READ(fsvf) skyvft 6774 ENDIF 6775 6776 IF ( nsvfl > 0 ) THEN 6777 IF ( .NOT. ALLOCATED( svf ) ) ALLOCATE( svf(ndsvf,nsvfl) ) 6778 IF ( .NOT. ALLOCATED( svfsurf ) ) ALLOCATE( svfsurf(idsvf,nsvfl) ) 6779 READ(fsvf) svf 6780 READ(fsvf) svfsurf 6781 ENDIF 6782 6783 IF ( plant_canopy .AND. ncsfl > 0 ) THEN 6784 IF ( .NOT. ALLOCATED( csf ) ) ALLOCATE( csf(ndcsf,ncsfl) ) 6785 IF ( .NOT. ALLOCATED( csfsurf ) ) ALLOCATE( csfsurf(idcsf,ncsfl) ) 6786 READ(fsvf) csf 6787 READ(fsvf) csfsurf 6763 6788 ENDIF 6764 6789 READ ( fsvf ) svf_code_field 6765 6790 6766 6791 IF ( TRIM(svf_code_field) /= TRIM(svf_code) ) THEN 6767 6768 6792 WRITE( message_string, * ) 'Wrong structure of binary svf file' 6793 CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 ) 6769 6794 ENDIF 6770 6795 ! 6771 !-- Close binary file6772 CALL close_file( fsvf )6796 !-- Close binary file 6797 CALL close_file( fsvf ) 6773 6798 6774 6799 ENDIF 6775 6800 #if defined( __parallel ) 6776 6801 CALL MPI_BARRIER( comm2d, ierr ) 6777 6802 #endif 6778 6803 ENDDO … … 6801 6826 WRITE ( fsvf ) rad_version 6802 6827 WRITE ( fsvf ) nsvfl, ncsfl, nsurfl 6803 WRITE ( fsvf ) skyvf 6804 WRITE ( fsvf ) skyvft 6805 WRITE ( fsvf ) svf 6806 WRITE ( fsvf ) svfsurf 6807 IF ( plant_canopy ) THEN 6828 IF ( nsurfl > 0 ) THEN 6829 WRITE ( fsvf ) skyvf 6830 WRITE ( fsvf ) skyvft 6831 ENDIF 6832 IF ( nsvfl > 0 ) THEN 6833 WRITE ( fsvf ) svf 6834 WRITE ( fsvf ) svfsurf 6835 ENDIF 6836 IF ( plant_canopy .AND. ncsfl > 0 ) THEN 6808 6837 WRITE ( fsvf ) csf 6809 6838 WRITE ( fsvf ) csfsurf
Note: See TracChangeset
for help on using the changeset viewer.
