Changeset 3117


Ignore:
Timestamp:
Jul 11, 2018 9:59:11 AM (4 years ago)
Author:
maronga
Message:

bugfixes in RRTMG interface

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3116 r3117  
    2828! -----------------
    2929! $Id$
     30! Bugfix: water vapor was not transfered to RRTMG when cloud_physics = .F.
     31! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
     32! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
     33!
     34! 3116 2018-07-10 14:31:58Z suehring
    3035! Output of long/shortwave radiation at surface
    3136!
     
    370375    USE control_parameters,                                                    &
    371376        ONLY:  cloud_droplets, cloud_physics, coupling_char, dz, g,            &
     377               humidity,                                                       &
    372378               initializing_actions, io_blocks, io_group,                      &
    373379               latitude, longitude, large_scale_forcing, lsf_surf,             &
     
    618624
    619625    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
    620                                            q_snd,       & !< mixing ratio from sounding data (kg/kg) - dummy at the moment
    621626                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
    622627                                           t_snd          !< actual temperature from sounding data (hPa)
     
    829834    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    830835    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf           !< array of total radiation flux incoming to minus outgoing from local surface
     836    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfemitlwl      !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
    831837
    832838!-- block variables needed for calculation of the plant canopy model inside the urban surface model
     
    28502856          CALL calc_mean_profile( pt, 4 )
    28512857          pt_av = hom(:, 1, 4, 0)
    2852 
     2858         
     2859          IF ( humidity )  THEN
     2860             CALL calc_mean_profile( q, 41 )
     2861             q_av  = hom(:, 1, 41, 0)
     2862          ENDIF
    28532863!
    28542864!--       Prepare profiles of temperature and H2O volume mixing ratio
     
    28562866
    28572867          IF ( cloud_physics )  THEN
    2858              CALL calc_mean_profile( q, 41 )
    2859              ! average  q is now in hom(:, 1, 41, 0)
    2860              q_av = hom(:, 1, 41, 0)
     2868
    28612869             CALL calc_mean_profile( ql, 54 )
    28622870             ! average ql is now in hom(:, 1, 54, 0)
     
    28722880                rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp       &
    28732881                                 )**.286_wp
    2874                 rrtm_h2ovmr(0,k) = 0._wp
    2875               ENDDO
     2882             ENDDO
     2883
     2884             IF ( humidity )  THEN
     2885                DO k = nzb+1, nzt+1
     2886                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q_av(k)
     2887                ENDDO
     2888             ELSE
     2889                rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
     2890             ENDIF
    28762891          ENDIF
    28772892
     
    29432958!--       Set surface temperature
    29442959          rrtm_tsfc = t_rad_urb
    2945 
     2960         
    29462961          IF ( lw_radiation )  THEN
    29472962             CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
     
    30813096                      rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp   &
    30823097                                       )**0.286_wp
    3083                       rrtm_h2ovmr(0,k) = 0.0_wp
    30843098                   ENDDO
     3099
     3100                   IF ( humidity )  THEN
     3101                      DO k = nzb+1, nzt+1
     3102                         rrtm_h2ovmr(0,k) =  mol_mass_air_d_wv * q(k,j,i)
     3103                      ENDDO   
     3104                   ELSE
     3105                      rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
     3106                   ENDIF
    30853107                ENDIF
    30863108
     
    37423764          DEALLOCATE( hyp_snd )
    37433765          DEALLOCATE( t_snd )
    3744           DEALLOCATE( q_snd  )
    37453766          DEALLOCATE ( rrtm_play )
    37463767          DEALLOCATE ( rrtm_plev )
     
    37483769          DEALLOCATE ( rrtm_tlev )
    37493770
    3750           DEALLOCATE ( rrtm_h2ovmr )
    37513771          DEALLOCATE ( rrtm_cicewp )
    37523772          DEALLOCATE ( rrtm_cldfr )
     
    38443864
    38453865!
    3846 !--    Save data above LES domain in hyp_snd, t_snd and q_snd
    3847 !--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
     3866!--    Save data above LES domain in hyp_snd, t_snd
    38483867       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
    38493868       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
    3850        ALLOCATE( q_snd(nzb+1:nzt_rad)   )
    38513869       hyp_snd = 0.0_wp
    38523870       t_snd = 0.0_wp
    3853        q_snd = 0.0_wp
    38543871
    38553872       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
     
    38913908       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
    38923909       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
    3893        ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
    38943910
    38953911       DO k = nzt+8, nzt_rad
    38963912          rrtm_tlay(0,k)   = t_snd(k)
    3897           rrtm_h2ovmr(0,k) = q_snd(k)
    38983913       ENDDO
    38993914       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
     
    39053920                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
    39063921       ENDDO
    3907        rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
    39083922
    39093923       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
     
    39944008       IMPLICIT NONE
    39954009
    3996        INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
     4010       INTEGER(iwp), PARAMETER :: num_trace_gases = 10 !< number of trace gases (absorbers)
    39974011
    39984012       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
    39994013           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
    4000                            'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
     4014                           'CFC11', 'CFC12', 'CFC22', 'CCL4 ', 'H2O'/)
    40014015
    40024016       INTEGER(iwp) :: id,     & !< NetCDF id
     
    40364050          DEALLOCATE ( rrtm_cfc22vmr )
    40374051          DEALLOCATE ( rrtm_ccl4vmr  )
     4052          DEALLOCATE ( rrtm_h2ovmr  )     
    40384053       ENDIF
    40394054
     
    40494064       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
    40504065       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
     4066       ALLOCATE ( rrtm_h2ovmr(0:0,1:nzt_rad+1)  )
    40514067
    40524068!
     
    42264242                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
    42274243
     4244             CASE ( 'H2O' )
     4245
     4246                rrtm_h2ovmr(0,:) = trace_path_tmp(:)
     4247               
    42284248             CASE DEFAULT
    42294249
     
    43884408     REAL(wp)                          :: area_surfl         !< total area of surfaces in local processor
    43894409     REAL(wp)                          :: area_surf          !< total area of surfaces in all processor
     4410     REAL(wp)                          :: area_horl          !< total horizontal area of domain in local processor
     4411     REAL(wp)                          :: area_hor           !< total horizontal area of domain in all processor
    43904412
    43914413
     
    46054627     surfoutsw = 0.0_wp
    46064628     surfoutlw = surfoutll
     4629     surfemitlwl = surfoutll
    46074630!        surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
    46084631
     
    48264849        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
    48274850     ENDDO
     4851!--  calculate horizontal area
     4852! !!! ATTENTION!!! uniform grid is assumed here
     4853     area_hor = (nx+1) * (ny+1) * dx * dy
    48284854!
    48294855!--  absorbed/received SW & LW and emitted LW energy of all physical
     
    48394865        d = surfl(id, i)
    48404866!--  received SW & LW
    4841         pinswl = pinswl + surfinsw(i) * facearea(d)
    4842         pinlwl = pinlwl + surfinlw(i) * facearea(d)
     4867        pinswl = pinswl + (surfinswdir(i) + surfinswdif(i)) * facearea(d)
     4868        pinlwl = pinlwl + surfinlwdif(i) * facearea(d)
    48434869!--   absorbed SW & LW
    48444870        pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
     
    48464872        pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
    48474873!--   emitted LW
    4848         pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
     4874        pemitlwl = pemitlwl + surfemitlwl(i) * facearea(d)
    48494875!--   emissivity and area sum
    48504876        emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
     
    48564882        pabsswl = pabsswl + SUM(pcbinsw)
    48574883        pabslwl = pabslwl + SUM(pcbinlw)
     4884        pinswl = pinswl + SUM(pcbinswdir) + SUM(pcbinswdif)
    48584885     ENDIF
    48594886!
     
    48714898     pinlw = pinlwl
    48724899     pabssw = pabsswl
    4873      pabslwl = pabslw
    4874      pemitlwl = pemitlw
     4900     pabslw = pabslwl
     4901     pemitlw = pemitlwl
    48754902     emiss_sum_surf = emiss_sum_surfl
    48764903     area_surf = area_surfl
     
    48784905
    48794906!--  (1) albedo
    4880      IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
    4881 
     4907     IF ( pinsw /= 0.0_wp )  &
     4908          albedo_urb = (pinsw - pabssw) / pinsw
    48824909!--  (2) average emmsivity
    4883      IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
    4884 
     4910     IF ( area_surf /= 0.0_wp ) &
     4911          emissivity_urb = emiss_sum_surf / area_surf
    48854912!--  (3) temperature
    4886      t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
     4913     t_rad_urb = ( (pemitlw - pabslw + emissivity_urb*pinlw) / &
     4914          (emissivity_urb*sigma_sb * area_hor) )**0.25_wp
    48874915     
    4888 
    48894916    CONTAINS
    48904917
     
    53565383       ALLOCATE( skyvf(nsurfl) )
    53575384       ALLOCATE( skyvft(nsurfl) )
     5385       ALLOCATE( surfemitlwl(nsurfl) )
    53585386
    53595387!
     
    53625390!--    For now set an arbitrary initial value.
    53635391       IF ( average_radiation )  THEN
    5364           albedo_urb = 0.5_wp
    5365           emissivity_urb = 0.5_wp
     5392          albedo_urb = 0.1_wp
     5393          emissivity_urb = 0.9_wp
    53665394          t_rad_urb = pt_surface
    53675395       ENDIF
Note: See TracChangeset for help on using the changeset viewer.