Changeset 4632 for palm/trunk/SOURCE/radiation_model_mod.f90

Timestamp:

Aug 5, 2020 10:51:52 AM (4 years ago)

Author:

pavelkrc

Message:

Fix average transparency for MRT factors, fix rtm_svf_, rtm_dif_ outputs

File:

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

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! - Bugfix: rtm_svf_, rtm_dif_ outputs
+! - Bugfix: correct average transparency for MRT factors
+!
+! 4624 2020-07-24 09:53:17Z raasch
 ! File re-formatted to follow the PALM coding standard
 !
@@ -1396 +1400 @@
 
     IF ( directional )  THEN
-       SELECT CASE ( var )
-          CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw', 'rtm_rad_inswdir',                 &
-                 'rtm_rad_inswdif', 'rtm_rad_inswref', 'rtm_rad_inlwdif', 'rtm_rad_inlwref',       &
-                 'rtm_rad_outsw', 'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
-             IF ( .NOT.  radiation )  THEN
-                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = 'W/m2'
-
-          CASE ( 'rtm_svf', 'rtm_dif', 'rtm_skyvf', 'rtm_skyvft', 'rtm_surfalb', 'rtm_surfemis' )
-             IF ( .NOT.  radiation )  THEN
-                message_string = 'output of "' // var // '" requires radiation = .TRUE.'
-                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
-             ENDIF
-             unit = '1'
-
-          CASE DEFAULT
-             unit = 'illegal'
-       END SELECT
-
+       IF ( var(1:8) == 'rtm_svf_'  .OR.  var(1:8) == 'rtm_dif_' )  THEN
+          IF ( .NOT.  radiation ) THEN
+             message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+             CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+          ENDIF
+          unit = '1'
+       ELSE
+          SELECT CASE ( var )
+             CASE ( 'rtm_rad_net', 'rtm_rad_insw', 'rtm_rad_inlw', 'rtm_rad_inswdir',              &
+                    'rtm_rad_inswdif', 'rtm_rad_inswref', 'rtm_rad_inlwdif', 'rtm_rad_inlwref',    &
+                    'rtm_rad_outsw', 'rtm_rad_outlw', 'rtm_rad_ressw', 'rtm_rad_reslw' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = 'W/m2'
+
+             CASE ( 'rtm_skyvf', 'rtm_skyvft', 'rtm_surfalb', 'rtm_surfemis' )
+                IF ( .NOT.  radiation ) THEN
+                   message_string = 'output of "' // var // '" requires radiation = .TRUE.'
+                   CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = '1'
+
+             CASE DEFAULT
+                unit = 'illegal'
+          END SELECT
+       ENDIF
     ELSE
        SELECT CASE ( var )
@@ -7921 +7932 @@
           CALL quicksort_itarget( itarget, vffrac, ztransp, 1, nzn * naz )
 !
+!--       For aggregation, we need fractions multiplied by transmissivities
+          ztransp(:) = vffrac(:) * ztransp(:)
+!
 !--       Find the first valid position
           itarg0 = 1
@@ -7934 +7948 @@
                 IF ( itarget(i+1) == itarget(i) )  THEN
                    vffrac(i+1) = vffrac(i+1) + vffrac(i)
+                   ztransp(i+1) = ztransp(i+1) + ztransp(i)
                    CYCLE
                 ENDIF
@@ -7968 +7983 @@
              amrtf(nmrtf)%isurfs = itarget(i)
              amrtf(nmrtf)%rsvf = vffrac(i)
-             amrtf(nmrtf)%rtransp = ztransp(i)
+             amrtf(nmrtf)%rtransp = ztransp(i) / vffrac(i)
           ENDDO ! itarg
 
@@ -11680 +11695 @@
             isurfs = svfsurf(2, isvf)
 
-            IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND. surf(iz,isurfs) == ks  &
-               .AND.  (surf(id,isurfs) == idsint_u .OR. surfl(id,isurfs) == idsint_l ) )  THEN
+            IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND.  surf(iz,isurfs) == ks &
+                 .AND.  ( surfl(id,isurflt) == idsint_u .OR. surfl(id,isurflt) == idsint_l ) ) THEN
 !
 !--            Correct source surface