Changeset 3843

Timestamp:

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

Author:

pavelkrc

Message:

Code review for radiation_model_mod.f90

File:

• palm/trunk/SOURCE/radiation_model_mod.f90 (modified) (18 diffs)

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -23 +23 @@
 ! Current revisions:
 ! ------------------
-! 
+! Code review.
 ! 
 ! Former revisions:
@@ -761 +761 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_in_xy_av  !< average of incoming shortwave radiation at surface
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
+
+    REAL(wp), PARAMETER :: emissivity_atm_clsky = 0.8_wp       !< emissivity of the clear-sky atmosphere
 !
 !-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)     
@@ -934 +936 @@
     INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
     INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
+    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
 
     INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
@@ -974 +978 @@
 
 !-- indices and sizes of urban and land surface models
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_l          !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_l           !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
-    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
-    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_linear     !< dtto (linearly allocated array)
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x, m]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_linear      !< dtto (linearly allocated array)
     INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
     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 
 !--       whether the sun is up
-          CALL calc_zenith     
-          ! readjust date and time to its initial value
+          CALL calc_zenith
+!
+!--       readjust date and time to its initial value
           CALL init_date_and_time
 !
@@ -2566 +2571 @@
 !
 !--          Horizontally aligned natural and urban surfaces
-             CALL calc_albedo( surf_lsm_h    )
-             CALL calc_albedo( surf_usm_h    )
+             CALL calc_albedo( surf_lsm_h )
+             CALL calc_albedo( surf_usm_h )
 !
 !--          Vertically aligned natural and urban surfaces
@@ -2652 +2657 @@
           IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
              ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-             rad_lw_in     = 0.0_wp
+             rad_lw_in = 0.0_wp
           ENDIF
 
@@ -2661 +2666 @@
           IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
              ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-            rad_lw_out    = 0.0_wp
+            rad_lw_out = 0.0_wp
           ENDIF
 
@@ -2917 +2922 @@
                 surf%rad_sw_out = albedo_urb * surf%rad_sw_in
                 
-                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(k+1))**4
+                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k+1))**4
 
                 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
@@ -2976 +2981 @@
                    IF ( bulk_cloud_model  .OR.  cloud_droplets  )  THEN
                       pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
-                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
+                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
                    ELSE
-                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt(k,j,i) * exner(k))**4
+                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt(k,j,i) * exner(k))**4
                    ENDIF
 
@@ -3093 +3098 @@
                 surf%rad_net = net_radiation
 
-                surf%rad_lw_in  = 0.8_wp * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
+                surf%rad_lw_in  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
 
                 surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4   &
@@ -3129 +3134 @@
                    IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
                       pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
-                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
+                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb * (pt1 * exner(k))**4
                    ELSE
-                      surf%rad_lw_in(m)  = 0.8_wp * sigma_sb *                 &
+                      surf%rad_lw_in(m)  = emissivity_atm_clsky * sigma_sb *                 &
                                              ( pt(k,j,i) * exner(k) )**4
                    ENDIF
@@ -4988 +4993 @@
      REAL(wp), DIMENSION(0:nsurf_type) :: costheta           !< direct irradiance factor of solar angle
      REAL(wp), DIMENSION(nz_urban_b:nz_urban_t)    :: pchf_prep          !< precalculated factor for canopy temperature tendency
-     REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: add to namelist or remove)
+     REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: synchronize with rotation_angle
+                                                             !< from netcdf_data_input_mod)
      REAL(wp)                          :: pc_box_area, pc_abs_frac, pc_abs_eff
      REAL(wp)                          :: asrc               !< area of source face
@@ -5437 +5443 @@
 !--     Transfer radiation arrays required for energy balance to the respective data types
      DO  i = 1, nsurfl
-        m  = surfl(5,i)
+        m  = surfl(im,i)
 !
 !--     (1) Urban surfaces
@@ -6096 +6102 @@
        ENDIF
 
-!--    fill surfl (the ordering of local surfaces given by the following
+!
+!--    Fill surfl (the ordering of local surfaces given by the following
 !--    cycles must not be altered, certain file input routines may depend
-!--    on it)
-       ALLOCATE(surfl_l(5*nsurfl))  ! is it necessary to allocate it with (5,nsurfl)?
-       surfl(1:5,1:nsurfl) => surfl_l(1:5*nsurfl)
+!--    on it).
+!
+!--    We allocate the array as linear and then use a two-dimensional pointer
+!--    into it, because some MPI implementations crash with 2D-allocated arrays.
+       ALLOCATE(surfl_linear(nidx_surf*nsurfl))
+       surfl(1:nidx_surf,1:nsurfl) => surfl_linear(1:nidx_surf*nsurfl)
        isurf = 0
        IF ( rad_angular_discretization )  THEN
@@ -6356 +6366 @@
        surfstart(numprocs) = k
        nsurf = k
-       ALLOCATE(surf_l(5*nsurf))
-       surf(1:5,1:nsurf) => surf_l(1:5*nsurf)
+!
+!--    We allocate the array as linear and then use a two-dimensional pointer
+!--    into it, because some MPI implementations crash with 2D-allocated arrays.
+       ALLOCATE(surf_linear(nidx_surf*nsurf))
+       surf(1:nidx_surf,1:nsurf) => surf_linear(1:nidx_surf*nsurf)
 
 #if defined( __parallel )
-       CALL MPI_AllGatherv(surfl_l, nsurfl*5, MPI_INTEGER, surf_l, nsurfs*5, &
-           surfstart(0:numprocs-1)*5, MPI_INTEGER, comm2d, ierr)
+       CALL MPI_AllGatherv(surfl_linear, nsurfl*nidx_surf, MPI_INTEGER,    &
+                           surf_linear, nsurfs*nidx_surf,                  &
+                           surfstart(0:numprocs-1)*nidx_surf, MPI_INTEGER, &
+                           comm2d, ierr)
        IF ( ierr /= 0 ) THEN
-           WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_l), nsurfl*5, &
-                      SIZE(surf_l), nsurfs*5, surfstart(0:numprocs-1)*5
+           WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_linear),    &
+                      nsurfl*nidx_surf, SIZE(surf_linear), nsurfs*nidx_surf, &
+                      surfstart(0:numprocs-1)*nidx_surf
            FLUSH(9)
        ENDIF
@@ -7378 +7394 @@
                 icsf = 1 !< reading index
                 kcsf = 1 !< writing index
-                DO while (icsf < npcsfl)
+                DO WHILE (icsf < npcsfl)
 !--                 here kpcsf(kcsf) already has values from kpcsf(icsf)
                     IF ( kpcsflt(3,icsf) == kpcsflt(3,icsf+1)  .AND.  &
@@ -8670 +8686 @@
 
     
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Grows the CSF array exponentially after it is full. During that, the ray
+!> canopy sink factors with common source face and target plant canopy grid
+!> cell are merged together so that the size doesn't grow out of control.
+!------------------------------------------------------------------------------!
     SUBROUTINE merge_and_grow_csf(newsize)
         INTEGER(iwp), INTENT(in)                :: newsize  !< new array size after grow, must be >= ncsfl