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

Timestamp:

Jul 11, 2018 9:59:11 AM (6 years ago)

Author:

maronga

Message:

bugfixes in RRTMG interface

File:

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

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Bugfix: water vapor was not transfered to RRTMG when cloud_physics = .F.
+! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
+! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
+! 
+! 3116 2018-07-10 14:31:58Z suehring
 ! Output of long/shortwave radiation at surface
 ! 
@@ -370 +375 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_droplets, cloud_physics, coupling_char, dz, g,            &
+               humidity,                                                       &
                initializing_actions, io_blocks, io_group,                      &
                latitude, longitude, large_scale_forcing, lsf_surf,             &
@@ -618 +624 @@
 
     REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
-                                           q_snd,       & !< mixing ratio from sounding data (kg/kg) - dummy at the moment
                                            rrtm_tsfc,   & !< dummy array for storing surface temperature
                                            t_snd          !< actual temperature from sounding data (hPa)
@@ -829 +834 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf           !< array of total radiation flux incoming to minus outgoing from local surface
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
 
 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
@@ -2850 +2856 @@
           CALL calc_mean_profile( pt, 4 )
           pt_av = hom(:, 1, 4, 0)
-
+          
+          IF ( humidity )  THEN
+             CALL calc_mean_profile( q, 41 )
+             q_av  = hom(:, 1, 41, 0)
+          ENDIF
 !
 !--       Prepare profiles of temperature and H2O volume mixing ratio
@@ -2856 +2866 @@
 
           IF ( cloud_physics )  THEN
-             CALL calc_mean_profile( q, 41 )
-             ! average  q is now in hom(:, 1, 41, 0)
-             q_av = hom(:, 1, 41, 0)
+
              CALL calc_mean_profile( ql, 54 )
              ! average ql is now in hom(:, 1, 54, 0)
@@ -2872 +2880 @@
                 rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
                                  )**.286_wp
-                rrtm_h2ovmr(0,k) = 0._wp
-              ENDDO
+             ENDDO
+
+             IF ( humidity )  THEN
+                DO k = nzb+1, nzt+1
+                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q_av(k)
+                ENDDO
+             ELSE
+                rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
+             ENDIF
           ENDIF
 
@@ -2943 +2958 @@
 !--       Set surface temperature
           rrtm_tsfc = t_rad_urb
-
+          
           IF ( lw_radiation )  THEN
              CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
@@ -3081 +3096 @@
                       rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp   &
                                        )**0.286_wp
-                      rrtm_h2ovmr(0,k) = 0.0_wp
                    ENDDO
+
+                   IF ( humidity )  THEN
+                      DO k = nzb+1, nzt+1
+                         rrtm_h2ovmr(0,k) =  mol_mass_air_d_wv * q(k,j,i)
+                      ENDDO   
+                   ELSE
+                      rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
+                   ENDIF
                 ENDIF
 
@@ -3742 +3764 @@
           DEALLOCATE( hyp_snd )
           DEALLOCATE( t_snd )
-          DEALLOCATE( q_snd  )
           DEALLOCATE ( rrtm_play )
           DEALLOCATE ( rrtm_plev )
@@ -3748 +3769 @@
           DEALLOCATE ( rrtm_tlev )
 
-          DEALLOCATE ( rrtm_h2ovmr )
           DEALLOCATE ( rrtm_cicewp )
           DEALLOCATE ( rrtm_cldfr )
@@ -3844 +3864 @@
 
 !
-!--    Save data above LES domain in hyp_snd, t_snd and q_snd
-!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
+!--    Save data above LES domain in hyp_snd, t_snd
        ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
        ALLOCATE( t_snd(nzb+1:nzt_rad)   )
-       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
        hyp_snd = 0.0_wp
        t_snd = 0.0_wp
-       q_snd = 0.0_wp
 
        hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
@@ -3891 +3908 @@
        ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
        ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
-       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
 
        DO k = nzt+8, nzt_rad
           rrtm_tlay(0,k)   = t_snd(k)
-          rrtm_h2ovmr(0,k) = q_snd(k)
        ENDDO
        rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
@@ -3905 +3920 @@
                              * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
        ENDDO
-       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
 
        rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
@@ -3994 +4008 @@
        IMPLICIT NONE
 
-       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
+       INTEGER(iwp), PARAMETER :: num_trace_gases = 10 !< number of trace gases (absorbers)
 
        CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
            trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
-                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
+                           'CFC11', 'CFC12', 'CFC22', 'CCL4 ', 'H2O'/)
 
        INTEGER(iwp) :: id,     & !< NetCDF id
@@ -4036 +4050 @@
           DEALLOCATE ( rrtm_cfc22vmr )
           DEALLOCATE ( rrtm_ccl4vmr  )
+          DEALLOCATE ( rrtm_h2ovmr  )     
        ENDIF
 
@@ -4049 +4064 @@
        ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
        ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
+       ALLOCATE ( rrtm_h2ovmr(0:0,1:nzt_rad+1)  )
 
 !
@@ -4226 +4242 @@
                 rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
 
+             CASE ( 'H2O' )
+
+                rrtm_h2ovmr(0,:) = trace_path_tmp(:)
+                
              CASE DEFAULT
 
@@ -4388 +4408 @@
      REAL(wp)                          :: area_surfl         !< total area of surfaces in local processor
      REAL(wp)                          :: area_surf          !< total area of surfaces in all processor
+     REAL(wp)                          :: area_horl          !< total horizontal area of domain in local processor
+     REAL(wp)                          :: area_hor           !< total horizontal area of domain in all processor
 
 
@@ -4605 +4627 @@
      surfoutsw = 0.0_wp
      surfoutlw = surfoutll
+     surfemitlwl = surfoutll
 !        surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
 
@@ -4826 +4849 @@
         IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
      ENDDO
+!--  calculate horizontal area
+! !!! ATTENTION!!! uniform grid is assumed here
+     area_hor = (nx+1) * (ny+1) * dx * dy
 !
 !--  absorbed/received SW & LW and emitted LW energy of all physical
@@ -4839 +4865 @@
         d = surfl(id, i)
 !--  received SW & LW
-        pinswl = pinswl + surfinsw(i) * facearea(d)
-        pinlwl = pinlwl + surfinlw(i) * facearea(d)
+        pinswl = pinswl + (surfinswdir(i) + surfinswdif(i)) * facearea(d)
+        pinlwl = pinlwl + surfinlwdif(i) * facearea(d)
 !--   absorbed SW & LW
         pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
@@ -4846 +4872 @@
         pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
 !--   emitted LW
-        pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
+        pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
 !--   emissivity and area sum
         emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
@@ -4856 +4882 @@
         pabsswl = pabsswl + SUM(pcbinsw)
         pabslwl = pabslwl + SUM(pcbinlw)
+        pinswl = pinswl + SUM(pcbinswdir) + SUM(pcbinswdif)
      ENDIF
 !
@@ -4871 +4898 @@
      pinlw = pinlwl
      pabssw = pabsswl
-     pabslwl = pabslw
-     pemitlwl = pemitlw
+     pabslw = pabslwl
+     pemitlw = pemitlwl
      emiss_sum_surf = emiss_sum_surfl
      area_surf = area_surfl
@@ -4878 +4905 @@
 
 !--  (1) albedo
-     IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
-
+     IF ( pinsw /= 0.0_wp )  &
+          albedo_urb = (pinsw - pabssw) / pinsw
 !--  (2) average emmsivity
-     IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
-
+     IF ( area_surf /= 0.0_wp ) &
+          emissivity_urb = emiss_sum_surf / area_surf
 !--  (3) temperature
-     t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
+     t_rad_urb = ( (pemitlw - pabslw + emissivity_urb*pinlw) / &
+          (emissivity_urb*sigma_sb * area_hor) )**0.25_wp
      
-
     CONTAINS
 
@@ -5356 +5383 @@
        ALLOCATE( skyvf(nsurfl) )
        ALLOCATE( skyvft(nsurfl) )
+       ALLOCATE( surfemitlwl(nsurfl) )
 
 !
@@ -5362 +5390 @@
 !--    For now set an arbitrary initial value.
        IF ( average_radiation )  THEN
-          albedo_urb = 0.5_wp
-          emissivity_urb = 0.5_wp
+          albedo_urb = 0.1_wp
+          emissivity_urb = 0.9_wp
           t_rad_urb = pt_surface
        ENDIF