Changeset 2213 for palm

Timestamp:

Apr 24, 2017 3:10:35 PM (8 years ago)

Author:

kanani

Message:

bugfix in PCM output, minor formatting and clean-up

Location:

palm/trunk/SOURCE

Files:

• plant_canopy_model_mod.f90 (modified) (8 diffs)
• urban_surface_mod.f90 (modified) (4 diffs)

palm/trunk/SOURCE/plant_canopy_model_mod.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: exchange of ghost points in array pc_heating_rate needed for output 
+! of pcm_heatrate, onetime ghost point exchange of lad_s after initialization.
+! Formatting and clean-up of subroutine pcm_read_plant_canopy_3d,
+! minor re-organization of canopy-heating initialization.
 ! 
 ! Former revisions:
@@ -346 +349 @@
 
       CASE ( 'pcm_heatrate' )
+         CALL exchange_horiz( pc_heating_rate, nbgp )
          IF ( av == 0 )  THEN
             DO  i = nxlg, nxrg
                DO  j = nysg, nyng
                   DO  k = nzb_s_inner(j,i), nz_do3d
-                     local_pf(i,j,k) = pc_heating_rate(k,j,i)
+                     local_pf(i,j,k) = pc_heating_rate(k-nzb_s_inner(j,i),j,i)
                   ENDDO
                ENDDO
@@ -358 +362 @@
     
       CASE ( 'pcm_lad' )
-
          IF ( av == 0 )  THEN
             DO  i = nxlg, nxrg
                DO  j = nysg, nyng
                   DO  k = nzb_s_inner(j,i), nz_do3d
-                     local_pf(i,j,k) = lad_s(k,j,i)
+                     local_pf(i,j,k) = lad_s(k-nzb_s_inner(j,i),j,i)
                   ENDDO
                ENDDO
@@ -542 +545 @@
 
        USE control_parameters,                                                 &
-           ONLY: coupling_char, dz, humidity, io_blocks, io_group,             &
-                 message_string, ocean, passive_scalar 
-
+           ONLY:  coupling_char, dz, humidity, io_blocks, io_group,            &
+                  message_string, ocean, passive_scalar 
+
+       USE control_parameters,                                                 &
+           ONLY:  urban_surface
 
        IMPLICIT NONE
@@ -661 +666 @@
 
 !
-!--    Allocate 3D-array for the leaf area density (lad_s). In case of a 
-!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
-!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
+!--    Allocate 3D-array for the leaf area density (lad_s).
        ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
-       IF ( cthf /= 0.0_wp )  THEN
-          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
-                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       ENDIF
 
 !
@@ -723 +721 @@
 !--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3, 
 !--    73–96). This approach has been applied e.g. by Shaw & Schumann (1992;
-!--    Bound.-Layer Meteorol. 61, 47–64)
-       IF ( cthf /= 0.0_wp )  THEN
+!--    Bound.-Layer Meteorol. 61, 47–64).
+!--    When using the urban surface model (USM), canopy heating (pc_heating_rate)
+!--    by radiation is calculated in the USM.
+       IF ( cthf /= 0.0_wp  .AND. .NOT.  urban_surface)  THEN
+
+          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
+                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
 !--       Piecewise calculation of the cumulative leaf area index by vertical
@@ -852 +855 @@
 !>
 !> dtype=1: leaf area density (lad_s)
-!> dtype=2: canopy drag coefficient (cdc)
-!> dtype=3: leaf scalar exchange coefficient (lsec)
-!> dtype=4: leaf surface concentration (lsc)
+!> dtype=2....n: some additional plant canopy input data quantity
 !>
 !> Zeros are added automatically above num_levels until top of domain.  Any
@@ -860 +861 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE pcm_read_plant_canopy_3d
-        USE control_parameters, &
-            ONLY: passive_scalar, message_string
-        IMPLICIT NONE
-
-        INTEGER(iwp)                            :: i, j, dtype, nzp, nzpltop, nzpl, kk
-        REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
-
-        lad_s = 0.0_wp
-        OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
-                action='READ', status='OLD', form='FORMATTED', err=515)
-        READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
-        ALLOCATE(col(1:nzp))
-        nzpltop = MIN(nzt+1, nzb+nzp-1)
-        nzpl = nzpltop - nzb + 1    !< no. of layers to assign
-
-        DO
-            READ(152, *, err=516, end=517) dtype, i, j, col(:)
-            IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
-
-            SELECT CASE (dtype)
-              CASE( 1 ) !< leaf area density
-                !-- only lad_s has flat z-coordinate, others have regular
-                kk = nzb_s_inner(j, i)
-                lad_s(nzb:nzpltop-kk, j, i) = col(1+kk:nzpl)
-!               CASE( 2 ) !< canopy drag coefficient
-!                 cdc(nzb:nzpltop, j, i) = col(1:nzpl)
-!               CASE( 3 ) !< leaf scalar exchange coefficient
-!                 lsec(nzb:nzpltop, j, i) = col(1:nzpl)
-!               CASE( 4 ) !< leaf surface concentration
-!                 lsc(nzb:nzpltop, j, i) = col(1:nzpl)
-              CASE DEFAULT
-                write(message_string, '(a,i2,a)')   &
-                    'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
-                CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
-            END SELECT
-        ENDDO
-
-515     message_string = 'error opening file PLANT_CANOPY_DATA_3D'
-        CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
-
-516     message_string = 'error reading file PLANT_CANOPY_DATA_3D'
-        CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
-
-517     CLOSE(152)
-        DEALLOCATE(col)
+    
+       USE control_parameters,                                                 &
+           ONLY:  message_string, passive_scalar
+
+       USE indices,                                                            &
+           ONLY:  nbgp
+           
+       IMPLICIT NONE
+
+       INTEGER(iwp)                        ::  dtype     !< type of input data (1=lad)
+       INTEGER(iwp)                        ::  i, j      !< running index
+       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
+       INTEGER(iwp)                        ::  nzpltop   !<
+       INTEGER(iwp)                        ::  nzpl      !<
+       
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
+
+!
+!--    Initialize lad_s array
+       lad_s = 0.0_wp
+       
+!
+!--    Open and read plant canopy input data
+       OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
+                 action='READ', status='OLD', form='FORMATTED', err=515)
+       READ(152, *, err=516, end=517) nzp   !< read first line = number of vertical layers
+       
+       ALLOCATE(col(0:nzp-1))
+
+       DO
+          READ(152, *, err=516, end=517) dtype, i, j, col(:)
+          IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
+
+             SELECT CASE (dtype)
+                CASE( 1 )   !< leaf area density
+!
+!--                 This is just the pure canopy layer assumed to be grounded to
+!--                 a flat domain surface. At locations where plant canopy sits 
+!--                 on top of any kind of topography, the vertical plant column 
+!--                 must be "lifted", which is done in SUBROUTINE pcm_tendency.
+                    lad_s(0:nzp-1, j, i) = col(0:nzp-1)
+                    
+                CASE DEFAULT
+                   write(message_string, '(a,i2,a)')   &
+                      'Unknown record type in file PLANT_CANOPY_DATA_3D: "', dtype, '"'
+                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0530', 1, 2, 0, 6, 0 )
+             END SELECT
+       ENDDO
+
+515    message_string = 'error opening file PLANT_CANOPY_DATA_3D'
+       CALL message( 'pcm_read_plant_canopy_3d', 'PA0531', 1, 2, 0, 6, 0 )
+
+516    message_string = 'error reading file PLANT_CANOPY_DATA_3D'
+       CALL message( 'pcm_read_plant_canopy_3d', 'PA0532', 1, 2, 0, 6, 0 )
+
+517    CLOSE(152)
+       DEALLOCATE(col)
+       
+       CALL exchange_horiz( lad_s, nbgp )
         
     END SUBROUTINE pcm_read_plant_canopy_3d

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -21 +21 @@
 ! Current revisions:
 ! ------------------
-! 
+! Removal of output quantities usm_lad and usm_canopy_hr
 ! 
 ! Former revisions:
@@ -2005 +2005 @@
                   var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis')  THEN
             unit = '1'
-        ELSE IF ( plant_canopy  .AND.  var(1:7) == 'usm_lad' )  THEN
-            unit = 'm2/m3'
-        ELSE IF ( plant_canopy  .AND.  var(1:13) == 'usm_canopy_hr' )  THEN
-            unit = 'K/s'
         ELSE
             unit = 'illegal'
@@ -2401 +2397 @@
                  ENDIF
               ENDDO
-
-          CASE ( 'usm_lad' )
-!--           leaf area density
-              DO i = nxl, nxr
-                 DO j = nys, nyn
-                     DO k = nzb_s_inner(j,i), nzut
-                         temp_pf(k,j,i) = lad_s(k-nzb_s_inner(j,i),j,i)
-                     ENDDO
-                 ENDDO
-              ENDDO
-              
-          CASE ( 'usm_canopy_hr' )
-!--           canopy heating rate
-              DO i = nxl, nxr
-                 DO j = nys, nyn
-                     DO k = nzb_s_inner(j,i), nzut
-                         temp_pf(k,j,i) = pc_heating_rate(k-nzb_s_inner(j,i),j,i)
-                     ENDDO
-                 ENDDO
-              ENDDO
              
           CASE DEFAULT
@@ -2471 +2447 @@
              var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.                     & 
              var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.                  &
-             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis'  .OR.            &
-             var(1:7) == 'usm_lad'  .OR.  var(1:13) == 'usm_canopy_hr' )  THEN
+             var(1:11) == 'usm_surfalb'  .OR.  var(1:12) == 'usm_surfemis' )  THEN
 
             found = .TRUE.