Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 1585 for palm


Ignore:
Timestamp:
Apr 30, 2015 7:05:52 AM (10 years ago)
Author:
maronga
Message:

Added support for RRTMG radiation code

Location:
palm/trunk
Files:
79 added
21 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r1576 r1585  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Added user_init_radiation.f90
    2323#
    2424# Former revisions:
    … …  
    223223        init_pegrid.f90 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \
    224224        interaction_droplets_ptq.f90 land_surface_model.f90 local_flush.f90 \
    225    local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
     225        local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
    226226        local_tremain_ini.f90 lpm.f90 lpm_advec.f90 lpm_boundary_conds.f90 \
    227227        lpm_calc_liquid_water_content.f90 lpm_collision_kernels.f90 \
    … …  
    253253        user_dvrp_coltab.f90 user_header.f90 user_init.f90 \
    254254        user_init_3d_model.f90 user_init_grid.f90 user_init_land_surface.f90 \
    255         user_init_plant_canopy.f90 user_last_actions.f90 user_lpm_advec.f90 \
     255        user_init_plant_canopy.f90 user_init_radiation.f90 \
     256        user_last_actions.f90 user_lpm_advec.f90 \
    256257        user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
    257258        user_parin.f90 user_read_restart_data.f90 \
    258259        user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
    259260        write_3d_binary.f90 write_var_list.f90
     261
    260262
    261263OBJS=$(SOURCES:.f90=.o)
    … …  
    427429        nudging.o plant_canopy_model.o production_e.o subsidence.o user_actions.o
    428430progress_bar.o: modules.o mod_kinds.o
    429 radiation_model.o : modules.o
     431radiation_model.o : modules.o mod_particle_attributes.o
    430432random_function.o: mod_kinds.o
    431433random_gauss.o: mod_kinds.o random_function.o random_generator_parallel.o
    … …  
    473475user_init_land_surface.o: modules.o mod_kinds.o user_module.o land_surface_model.o
    474476user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o
     477user_init_radiation.o: modules.o mod_kinds.o user_module.o radiation_model.o
    475478user_last_actions.o: modules.o mod_kinds.o user_module.o
    476479user_lpm_advec.o: modules.o mod_kinds.o user_module.o

  • palm/trunk/SOURCE/average_3d_data.f90

    r1556 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    8888
    8989    USE radiation_model_mod,                                                   &
    90         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
     90        ONLY:  rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,       &
     91               rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,             &
     92               rad_sw_out_av
    9193
    9294    IMPLICIT NONE
    … …  
    349351             ENDDO
    350352
    351          CASE ( 'rad_sw_in*' )
    352              DO  i = nxlg, nxrg
    353                 DO  j = nysg, nyng
    354                    rad_sw_in_av(j,i) = rad_sw_in_av(j,i) / REAL( average_count_3d, KIND=wp )
    355                 ENDDO
    356              ENDDO
    357 
    358353         CASE ( 'rad_net*' )
    359354             DO  i = nxlg, nxrg
    360355                DO  j = nysg, nyng
    361356                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
     357                ENDDO
     358             ENDDO
     359
     360          CASE ( 'rad_lw_in' )
     361             DO  i = nxlg, nxrg
     362                DO  j = nysg, nyng
     363                   DO  k = nzb, nzt+1
     364                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     365                   ENDDO
     366                ENDDO
     367             ENDDO
     368
     369          CASE ( 'rad_lw_out' )
     370             DO  i = nxlg, nxrg
     371                DO  j = nysg, nyng
     372                   DO  k = nzb, nzt+1
     373                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     374                   ENDDO
     375                ENDDO
     376             ENDDO
     377
     378          CASE ( 'rad_sw_in' )
     379             DO  i = nxlg, nxrg
     380                DO  j = nysg, nyng
     381                   DO  k = nzb, nzt+1
     382                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     383                   ENDDO
     384                ENDDO
     385             ENDDO
     386
     387          CASE ( 'rad_sw_out' )
     388             DO  i = nxlg, nxrg
     389                DO  j = nysg, nyng
     390                   DO  k = nzb, nzt+1
     391                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     392                   ENDDO
    362393                ENDDO
    363394             ENDDO

  • palm/trunk/SOURCE/check_parameters.f90

    r1576 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for RRTMG
    2323!
    2424! Former revisions:
    … …  
    11431143
    11441144    IF ( radiation )  THEN
    1145        IF ( radiation_scheme == 'constant' )  THEN
    1146           irad_scheme = 0
    1147        ELSEIF ( radiation_scheme == 'clear-sky' )  THEN
    1148           irad_scheme = 1
    1149        ELSEIF ( radiation_scheme == 'rrtm' )  THEN
    1150           irad_scheme = 2
    1151        ELSE
     1145       IF ( radiation_scheme /= 'constant'  .AND.                              &
     1146            radiation_scheme /= 'clear-sky' .AND.                              &
     1147            radiation_scheme /= 'rrtmg' )  THEN
    11521148          message_string = 'unknown radiation_scheme = '//                     &
    11531149                           TRIM( radiation_scheme )
    11541150          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
     1151       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     1152#if ! defined ( __rrtmg )
     1153          message_string = 'radiation_scheme = "rrtmg" requires ' //           &
     1154                           'compilation of PALM with pre-processor ' //        &
     1155                           'directive -D__rrtmg'
     1156          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
     1157#endif
     1158       ENDIF
     1159       IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND.                &
     1160            radiation_scheme == 'clear-sky')  THEN
     1161          message_string = 'radiation_scheme = "clear-sky" in combination' //  &
     1162                           'with albedo_type = 0 requires setting of albedo'// &
     1163                           ' /= 9999999.9'
     1164          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
    11551165       ENDIF
    11561166    ENDIF
    … …  
    30583068             ENDIF
    30593069
     3070          CASE ( 'rad_lw_in' )
     3071             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3072                message_string = 'data_output_pr = ' //                        &
     3073                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     3074                                 'lable for radiation = .FALSE. or ' //        &
     3075                                 'radiation_scheme = "constant"'
     3076                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     3077             ELSE
     3078                dopr_index(i) = 102
     3079                dopr_unit(i)  = 'W/m2'
     3080                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
     3081             ENDIF
     3082
     3083          CASE ( 'rad_lw_out' )
     3084             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3085                message_string = 'data_output_pr = ' //                        &
     3086                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     3087                                 'lable for radiation = .FALSE. or ' //        &
     3088                                 'radiation_scheme = "constant"'
     3089                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     3090             ELSE
     3091                dopr_index(i) = 103
     3092                dopr_unit(i)  = 'W/m2'
     3093                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
     3094             ENDIF
     3095
     3096          CASE ( 'rad_sw_in' )
     3097             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3098                message_string = 'data_output_pr = ' //                        &
     3099                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     3100                                 'lable for radiation = .FALSE. or ' //        &
     3101                                 'radiation_scheme = "constant"'
     3102                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     3103             ELSE
     3104                dopr_index(i) = 104
     3105                dopr_unit(i)  = 'W/m2'
     3106                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
     3107             ENDIF
     3108
     3109          CASE ( 'rad_sw_out')
     3110             IF ( (.NOT. radiation) .OR. radiation_scheme == 'constant' )  THEN
     3111                message_string = 'data_output_pr = ' //                        &
     3112                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
     3113                                 'lable for radiation = .FALSE. or ' //        &
     3114                                 'radiation_scheme = "constant"'
     3115                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
     3116             ELSE
     3117                dopr_index(i) = 105
     3118                dopr_unit(i)  = 'W/m2'
     3119                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
     3120             ENDIF
    30603121
    30613122          CASE DEFAULT
    … …  
    32573318             unit = 'kg/kg'
    32583319
     3320
     3321          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
     3322             IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' )  THEN
     3323                message_string = '"output of "' // TRIM( var ) // '" requi' //  &
     3324                                 'res radiation = .TRUE. and ' //              &
     3325                                 'radiation_scheme = "rrtmg"'
     3326                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
     3327             ENDIF
     3328             unit = 'W/m2'
     3329
    32593330          CASE ( 'rho' )
    32603331             IF ( .NOT. ocean )  THEN
    … …  
    32933364                 'm_liq_eb*', 'pra*', 'prr*', 'qsws*', 'qsws_eb*',             &
    32943365                 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*', 'rad_net*',  &
    3295                  'rad_sw_in*', 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*',        &
    3296                  'u*', 'z0*', 'z0h*' )
     3366                 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*', 'rrtm_asdir*',   &
     3367                 'r_a*', 'r_s*', 'shf*', 'shf_eb*', 't*', 'u*', 'z0*', 'z0h*' )
    32973368             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    32983369                message_string = 'illegal value for data_output: "' //         &
    … …  
    33013372                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    33023373             ENDIF
     3374             IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" )  THEN
     3375                IF ( TRIM( var ) == 'rrtm_aldif*' .OR.                         &
     3376                     TRIM( var ) == 'rrtm_aldir*' .OR.                         &
     3377                     TRIM( var ) == 'rrtm_asdif*' .OR.                         &
     3378                     TRIM( var ) == 'rrtm_asdir*'      )                       &
     3379                THEN
     3380                   message_string = 'output of "' // TRIM( var ) // '" require'&
     3381                                    // 's radiation = .TRUE. and radiation_sch'&
     3382                                    // 'eme = "rrtmg"'
     3383                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
     3384                ENDIF
     3385             ENDIF
     3386
    33033387             IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT. land_surface )  THEN
    33043388                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    … …  
    34053489             IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
    34063490             IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
    3407              IF ( TRIM( var ) == 'rad_net*')       unit = 'W/m2'     
    3408              IF ( TRIM( var ) == 'rad_sw_in*')     unit = 'W/m2'
     3491             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
     3492             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
     3493             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = ''
     3494             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = ''
     3495             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
    34093496             IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
    34103497             IF ( TRIM( var ) == 'r_s*')     unit = 's/m'

  • palm/trunk/SOURCE/data_output_2d.f90

    r1556 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for RRTMG
    2323!
    2424! Former revisions:
    … …  
    154154
    155155    USE radiation_model_mod,                                                   &
    156         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
     156        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
     157               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
     158               rad_lw_out_av
    157159
    158160    IMPLICIT NONE
    … …  
    913915                level_z(nzb+1) = zu(nzb+1)
    914916
    915              CASE ( 'rad_sw_in*_xy' )        ! 2d-array
    916                 IF ( av == 0 ) THEN
    917                    DO  i = nxlg, nxrg
    918                       DO  j = nysg, nyng
    919                          local_pf(i,j,nzb+1) =  rad_sw_in(j,i)
    920                       ENDDO
    921                    ENDDO
    922                 ELSE
    923                    DO  i = nxlg, nxrg
    924                       DO  j = nysg, nyng
    925                          local_pf(i,j,nzb+1) =  rad_sw_in_av(j,i)
    926                       ENDDO
    927                    ENDDO
    928                 ENDIF
    929                 resorted = .TRUE.
    930                 two_d = .TRUE.
    931                 level_z(nzb+1) = zu(nzb+1)
     917
     918             CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
     919                IF ( av == 0 )  THEN
     920                   to_be_resorted => rad_lw_in
     921                ELSE
     922                   to_be_resorted => rad_lw_in_av
     923                ENDIF
     924
     925             CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
     926                IF ( av == 0 )  THEN
     927                   to_be_resorted => rad_lw_out
     928                ELSE
     929                   to_be_resorted => rad_lw_out_av
     930                ENDIF
     931
     932             CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
     933                IF ( av == 0 )  THEN
     934                   to_be_resorted => rad_sw_in
     935                ELSE
     936                   to_be_resorted => rad_sw_in_av
     937                ENDIF
     938
     939             CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
     940                IF ( av == 0 )  THEN
     941                   to_be_resorted => rad_sw_out
     942                ELSE
     943                   to_be_resorted => rad_sw_out_av
     944                ENDIF
    932945
    933946             CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )

  • palm/trunk/SOURCE/data_output_3d.f90

    r1552 r1585  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Added support for RRTMG
    2323!
    2424! Former revisions:
    … …  
    128128       
    129129    USE pegrid
     130
     131    USE radiation_model_mod,                                                   &
     132        ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
     133               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
     134
    130135
    131136    IMPLICIT NONE
    … …  
    467472             ELSE
    468473                to_be_resorted => qv_av
     474             ENDIF
     475
     476          CASE ( 'rad_sw_in' )
     477             IF ( av == 0 )  THEN
     478                to_be_resorted => rad_sw_in
     479             ELSE
     480                to_be_resorted => rad_sw_in_av
     481             ENDIF
     482
     483          CASE ( 'rad_sw_out' )
     484             IF ( av == 0 )  THEN
     485                to_be_resorted => rad_sw_out
     486             ELSE
     487                to_be_resorted => rad_sw_out_av
     488             ENDIF
     489
     490          CASE ( 'rad_lw_in' )
     491             IF ( av == 0 )  THEN
     492                to_be_resorted => rad_lw_in
     493             ELSE
     494                to_be_resorted => rad_lw_in_av
     495             ENDIF
     496
     497          CASE ( 'rad_lw_out' )
     498             IF ( av == 0 )  THEN
     499                to_be_resorted => rad_lw_out
     500             ELSE
     501                to_be_resorted => rad_lw_out_av
    469502             ENDIF
    470503

  • palm/trunk/SOURCE/data_output_mask.f90

    r1439 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for RRTMG
    2323!
    2424! Former revisions:
    … …  
    108108   
    109109    USE pegrid
     110
     111    USE radiation_model_mod,                                                   &
     112        ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
     113               rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
    110114
    111115    IMPLICIT NONE
    … …  
    393397             ENDIF
    394398
     399          CASE ( 'rad_sw_in' )
     400             IF ( av == 0 )  THEN
     401                to_be_resorted => rad_sw_in
     402             ELSE
     403                to_be_resorted => rad_sw_in_av
     404             ENDIF
     405
     406          CASE ( 'rad_sw_out' )
     407             IF ( av == 0 )  THEN
     408                to_be_resorted => rad_sw_out
     409             ELSE
     410                to_be_resorted => rad_sw_out_av
     411             ENDIF
     412
     413          CASE ( 'rad_lw_in' )
     414             IF ( av == 0 )  THEN
     415                to_be_resorted => rad_lw_in
     416             ELSE
     417                to_be_resorted => rad_lw_in_av
     418             ENDIF
     419
     420          CASE ( 'rad_lw_out' )
     421             IF ( av == 0 )  THEN
     422                to_be_resorted => rad_lw_out
     423             ELSE
     424                to_be_resorted => rad_lw_out_av
     425             ENDIF
     426
    395427          CASE ( 'rho' )
    396428             IF ( av == 0 )  THEN

  • palm/trunk/SOURCE/flow_statistics.f90

    r1572 r1585  
    2121! Current revisions:
    2222! -----------------
    23 !
     23! Added output of timeseries and profiles for RRTMG
    2424!
    2525! Former revisions:
    … …  
    178178
    179179    USE radiation_model_mod,                                                   &
    180         ONLY :  dots_rad, rad_net, rad_sw_in, radiation
    181    
     180        ONLY:  dots_rad, radiation, radiation_scheme, rad_net,                 &
     181               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
     182
     183#if defined ( __rrtmg )
     184    USE radiation_model_mod,                                                   &
     185        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
     186#endif
     187 
    182188    USE statistics
    183189
    … …  
    741747             IF ( radiation )  THEN
    742748                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
    743                 sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_sw_in(j,i)
     749                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_lw_in(nzb,j,i)
     750                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_lw_out(nzb,j,i)
     751                sums_l(nzb,104,tn)  = sums_l(nzb,104,tn)  + rad_sw_in(nzb,j,i)
     752                sums_l(nzb,105,tn)  = sums_l(nzb,105,tn)  + rad_sw_out(nzb,j,i)
     753
     754#if defined ( __rrtmg )
     755                IF ( radiation_scheme == 'rrtmg' )  THEN
     756                   sums_l(nzb,106,tn)  = sums_l(nzb,106,tn)  + rrtm_aldif(0,j,i)
     757                   sums_l(nzb,107,tn)  = sums_l(nzb,107,tn)  + rrtm_aldir(0,j,i)
     758                   sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_asdif(0,j,i)
     759                   sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_asdir(0,j,i)
     760                ENDIF
     761#endif
     762
    744763             ENDIF
    745764
    … …  
    11261145       ENDIF
    11271146       
    1128 
     1147       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
     1148          !$OMP DO
     1149          DO  i = nxl, nxr
     1150             DO  j =  nys, nyn
     1151                DO  k = nzb_s_inner(j,i)+1, nzt+1
     1152                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
     1153                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
     1154                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
     1155                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
     1156                ENDDO
     1157             ENDDO
     1158          ENDDO
     1159       ENDIF
    11291160!
    11301161!--    Calculate the user-defined profiles
    … …  
    11861217          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
    11871218          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
    1188           sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
    1189           sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     1219          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
     1220          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    11901221       ENDDO
    11911222
    … …  
    13191350
    13201351       IF ( radiation )  THEN
    1321           hom(:,1,101 ,sr) = sums(:,101)            ! rad_net
    1322           hom(:,1,102,sr)  = sums(:,102)            ! rad_sw_in
     1352          hom(:,1,101,sr) = sums(:,101)            ! rad_net
     1353          hom(:,1,102,sr) = sums(:,102)            ! rad_lw_in
     1354          hom(:,1,103,sr) = sums(:,103)            ! rad_lw_out
     1355          hom(:,1,104,sr) = sums(:,104)            ! rad_sw_in
     1356          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
     1357
     1358#if defined ( __rrtmg )
     1359          IF ( radiation_scheme == 'rrtmg' )  THEN
     1360             hom(:,1,106,sr) = sums(:,106)            ! rrtm_aldif
     1361             hom(:,1,107,sr) = sums(:,107)            ! rrtm_aldir
     1362             hom(:,1,108,sr) = sums(:,108)            ! rrtm_asdif
     1363             hom(:,1,109,sr) = sums(:,109)            ! rrtm_asdir
     1364          ENDIF
     1365#endif
     1366
    13231367       ENDIF
    13241368
    … …  
    14831527!--    Collect radiation model timeseries
    14841528       IF ( radiation )  THEN
    1485           ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)           ! rad_net
    1486           ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
     1529          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
     1530          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
     1531          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
     1532          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_lw_in
     1533          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_lw_out
     1534
     1535#if defined ( __rrtmg )
     1536          IF ( radiation_scheme == 'rrtmg' )  THEN
     1537             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
     1538             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
     1539             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
     1540             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
     1541          ENDIF
     1542#endif
     1543
    14871544       ENDIF
    14881545
    … …  
    29513008
    29523009
     3010       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
     3011          !$OMP DO
     3012          DO  i = nxl, nxr
     3013             DO  j =  nys, nyn
     3014                DO  k = nzb_s_inner(j,i)+1, nzt+1
     3015                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
     3016                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
     3017                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
     3018                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
     3019                ENDDO
     3020             ENDDO
     3021          ENDDO
     3022       ENDIF
     3023
    29533024!
    29543025!--    Calculate the user-defined profiles
    … …  
    30123083          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
    30133084          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
    3014           sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
    3015           sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
     3085          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
     3086          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    30163087       ENDDO
    30173088
    … …  
    32883359!--    Collect radiation model timeseries
    32893360       IF ( radiation )  THEN
    3290           ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)           ! rad_net
    3291           ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
     3361          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
     3362          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
     3363          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
     3364          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_lw_in
     3365          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_lw_out
     3366
     3367#if defined ( __rrtmg )
     3368          IF ( radiation_scheme == 'rrtmg' )  THEN
     3369             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
     3370             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
     3371             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
     3372             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
     3373          ENDIF
     3374#endif
     3375
    32923376       ENDIF
    32933377

  • palm/trunk/SOURCE/header.f90

    r1576 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Further output for radiation model(s).
    2323!
    2424! Former revisions:
    … …  
    245245
    246246    USE radiation_model_mod,                                                   &
    247         ONLY:  albedo, day_init, dt_radiation, lambda, net_radiation,          &
    248                radiation, radiation_scheme, time_utc_init
     247        ONLY:  albedo, albedo_type, albedo_type_name, constant_albedo,         &
     248               day_init, dt_radiation, lambda, lw_radiation, net_radiation,    &
     249               radiation, radiation_scheme, sw_radiation, time_utc_init
    249250   
    250251    USE spectrum,                                                              &
    … …  
    927928    IF ( radiation )  THEN
    928929!
    929 !--    Write land surface model header
     930!--    Write radiation model header
    930931       WRITE( io, 444 )
    931932
    … …  
    934935       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
    935936          WRITE( io, 446 )
    936        ELSE
    937           WRITE( io, 447 ) radiation_scheme
    938        ENDIF
    939 
    940        WRITE( io, 448 ) albedo
     937       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
     938          WRITE( io, 447 )
     939          IF ( .NOT. lw_radiation )  WRITE( io, 458 )
     940          IF ( .NOT. sw_radiation )  WRITE( io, 459 )
     941       ENDIF
     942
     943       IF ( albedo_type == 0 )  THEN
     944          WRITE( io, 448 ) albedo
     945       ELSE
     946          WRITE( io, 456 ) albedo_type_name(albedo_type)
     947       ENDIF
     948       IF ( constant_albedo )  THEN
     949          WRITE( io, 457 )
     950       ENDIF
    941951       WRITE( io, 449 ) dt_radiation
    942 
    943952    ENDIF
    944953
    … …  
    22622271446 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',  &
    22632272                   ' default)')
    2264 447 FORMAT (' --> Radiation scheme:', A)
    2265 448 FORMAT (/'    Surface albedo: albedo = ', F5.3)
    2266 449 FORMAT  ('    Timestep: dt_radiation = ', F5.2, '  s')
     2273447 FORMAT (' --> RRTMG scheme is used')
     2274448 FORMAT (/'     User-specific surface albedo: albedo = ', F5.3)
     2275449 FORMAT  ('     Timestep: dt_radiation = ', F5.2, '  s')
    22672276
    22682277450 FORMAT (//' LES / Turbulence quantities:'/ &
    … …  
    22732282454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
    22742283455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
     2284456 FORMAT (/'    Albedo is set for land surface type: ', A)
     2285457 FORMAT (/'    --> Albedo is fixed during the run')
     2286458 FORMAT (/'    --> Longwave radiation is disabled')
     2287459 FORMAT (/'    --> Shortwave radiation is disabled.')
    22752288470 FORMAT (//' Actions during the simulation:'/ &
    22762289              ' -----------------------------'/)

  • palm/trunk/SOURCE/init_3d_model.f90

    r1576 r1585  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Initialization of radiation code is now done after LSM initializtion
    2323!
    2424! Former revisions:
    … …  
    16061606    IF ( particle_advection )  CALL lpm_init
    16071607
     1608!
     1609!-- If required, initialize quantities needed for the LSM
     1610    IF ( land_surface )  THEN
     1611       CALL location_message( 'initializing land surface model', .FALSE. )
     1612       CALL init_lsm
     1613       CALL location_message( 'finished', .TRUE. )
     1614    ENDIF
    16081615
    16091616!
    16101617!-- If required, initialize radiation model
    16111618    IF ( radiation )  THEN
     1619       CALL location_message( 'initializing radiation model', .FALSE. )
    16121620       CALL init_radiation
    1613     ENDIF
    1614 
    1615 !
    1616 !-- If required, initialize quantities needed for the LSM
    1617     IF ( land_surface )  THEN
    1618        CALL init_lsm
     1621       CALL location_message( 'finished', .TRUE. )
    16191622    ENDIF
    16201623

  • palm/trunk/SOURCE/land_surface_model.f90

    r1572 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Modifications for RRTMG. Changed tables to PARAMETER type.
    2323!
    2424! Former revisions:
    … …  
    8080!------------------------------------------------------------------------------!
    8181     USE arrays_3d,                                                            &
    82          ONLY:  pt, pt_p, q, q_p, qsws, rif, shf, ts, us, z0, z0h
     82         ONLY:  pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h
    8383
    8484     USE cloud_parameters,                                                     &
    … …  
    9090                max_masks, precipitation, pt_surface, rho_surface,             &
    9191                roughness_length, surface_pressure, timestep_scheme, tsc,      &
    92                 z0h_factor
     92                z0h_factor, time_since_reference_point
    9393
    9494     USE indices,                                                              &
    95          ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner
     95         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
    9696
    9797     USE kinds
    … …  
    101101
    102102     USE radiation_model_mod,                                                  &
    103          ONLY:  irad_scheme, rad_net, rad_sw_in, sigma_sb
     103         ONLY:  radiation_scheme, rad_net, rad_sw_in, sigma_sb
    104104
    105105     USE statistics,                                                           &
    … …  
    243243              qsws_veg_eb,      & !: surface flux of latent heat (vegetation portion)
    244244              qsws_veg_eb_av,   & !: average of qsws_veg_eb
     245              rad_net_l,        & !: local copy of rad_net (net radiation at surface)
    245246              r_a,              & !: aerodynamic resistance
    246247              r_a_av,           & !: avergae of r_a
    … …  
    276277              t_soil_av, t_soil_1, t_soil_2,                                   &
    277278              m_soil_av, m_soil_1, m_soil_2
    278 
    279 
    280279#endif
    281280
    … …  
    289288!
    290289!-- Predefined Land surface classes (veg_type)
    291     CHARACTER(26), DIMENSION(0:19) :: veg_type_name = (/          &
    292                                    'user defined',                & ! 0
    293                                    'crops, mixed farming',        & !  1
    294                                    'short grass',                 & !  2
    295                                    'evergreen needleleaf trees',  & !  3
    296                                    'deciduous needleleaf trees',  & !  4
    297                                    'evergreen broadleaf trees' ,  & !  5
    298                                    'deciduous broadleaf trees',   & !  6
    299                                    'tall grass',                  & !  7
    300                                    'desert',                      & !  8
    301                                    'tundra',                      & !  9
    302                                    'irrigated crops',             & ! 10
    303                                    'semidesert',                  & ! 11
    304                                    'ice caps and glaciers' ,      & ! 12
    305                                    'bogs and marshes',            & ! 13
    306                                    'inland water',                & ! 14
    307                                    'ocean',                       & ! 15
    308                                    'evergreen shrubs',            & ! 16
    309                                    'deciduous shrubs',            & ! 17
    310                                    'mixed forest/woodland',       & ! 18
    311                                    'interrupted forest'           & ! 19
    312                                                        /)
     290    CHARACTER(26), DIMENSION(0:19), PARAMETER :: veg_type_name = (/ &
     291                                   'user defined',                  & ! 0
     292                                   'crops, mixed farming',          & !  1
     293                                   'short grass',                   & !  2
     294                                   'evergreen needleleaf trees',    & !  3
     295                                   'deciduous needleleaf trees',    & !  4
     296                                   'evergreen broadleaf trees' ,    & !  5
     297                                   'deciduous broadleaf trees',     & !  6
     298                                   'tall grass',                    & !  7
     299                                   'desert',                        & !  8
     300                                   'tundra',                        & !  9
     301                                   'irrigated crops',               & ! 10
     302                                   'semidesert',                    & ! 11
     303                                   'ice caps and glaciers' ,        & ! 12
     304                                   'bogs and marshes',              & ! 13
     305                                   'inland water',                  & ! 14
     306                                   'ocean',                         & ! 15
     307                                   'evergreen shrubs',              & ! 16
     308                                   'deciduous shrubs',              & ! 17
     309                                   'mixed forest/woodland',         & ! 18
     310                                   'interrupted forest'             & ! 19
     311                                                                 /)
    313312
    314313!
    315314!-- Soil model classes (soil_type)
    316     CHARACTER(12), DIMENSION(0:7) :: soil_type_name = (/ &
    317                                    'user defined',        & ! 0
    318                                    'coarse',              & ! 1
    319                                    'medium',              & ! 2
    320                                    'medium-fine',         & ! 3
    321                                    'fine',                & ! 4
    322                                    'very fine' ,          & ! 5
    323                                    'organic',             & ! 6
    324                                    'loamy (CH)'           & ! 7
    325                                                         /)
     315    CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
     316                                   'user defined',                  & ! 0
     317                                   'coarse',                        & ! 1
     318                                   'medium',                        & ! 2
     319                                   'medium-fine',                   & ! 3
     320                                   'fine',                          & ! 4
     321                                   'very fine' ,                    & ! 5
     322                                   'organic',                       & ! 6
     323                                   'loamy (CH)'                     & ! 7
     324                                                                 /)
    326325!
    327326!-- Land surface parameters according to the respective classes (veg_type)
    … …  
    330329!-- Land surface parameters I
    331330!--                          r_canopy_min,     lai,   c_veg,     g_d
    332     REAL(wp), DIMENSION(0:3,1:19) :: veg_pars = RESHAPE( (/          &
    333                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & !  1
    334                                  110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & !  2
    335                                  500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  3
    336                                  500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  4
    337                                  175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & !  5
    338                                  240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp, & !  6
    339                                  100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp, & !  7
    340                                  250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp, & !  8
    341                                   80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp, & !  9
    342                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 10
    343                                  150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp, & ! 11
    344                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
    345                                  240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp, & ! 13
    346                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
    347                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
    348                                  225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp, & ! 16
    349                                  225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17
    350                                  250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18
    351                                  175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp  & ! 19
     331    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: veg_pars = RESHAPE( (/ &
     332                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
     333                                 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
     334                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
     335                                 500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
     336                                 175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
     337                                 240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
     338                                 100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
     339                                 250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
     340                                  80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
     341                                 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
     342                                 150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
     343                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
     344                                 240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
     345                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
     346                                   0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
     347                                 225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
     348                                 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
     349                                 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
     350                                 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp   & ! 19
    352351                                 /), (/ 4, 19 /) )
    353352
    354353!
    355354!-- Land surface parameters II          z0,         z0h
    356     REAL(wp), DIMENSION(0:1,1:19) :: roughness_par = RESHAPE( (/ &
    357                                    0.25_wp,  0.25E-2_wp,         & !  1
    358                                    0.20_wp,  0.20E-2_wp,         & !  2
    359                                    2.00_wp,     2.00_wp,         & !  3
    360                                    2.00_wp,     2.00_wp,         & !  4
    361                                    2.00_wp,     2.00_wp,         & !  5
    362                                    2.00_wp,     2.00_wp,         & !  6
    363                                    0.47_wp,  0.47E-2_wp,         & !  7
    364                                   0.013_wp, 0.013E-2_wp,         & !  8
    365                                   0.034_wp, 0.034E-2_wp,         & !  9
    366                                     0.5_wp,  0.50E-2_wp,         & ! 10
    367                                    0.17_wp,  0.17E-2_wp,         & ! 11
    368                                  1.3E-3_wp,   1.3E-4_wp,         & ! 12
    369                                    0.83_wp,  0.83E-2_wp,         & ! 13
    370                                    0.00_wp,  0.00E-2_wp,         & ! 14
    371                                    0.00_wp,  0.00E-2_wp,         & ! 15
    372                                    0.10_wp,  0.10E-2_wp,         & ! 16
    373                                    0.25_wp,  0.25E-2_wp,         & ! 17
    374                                    2.00_wp,  2.00E-2_wp,         & ! 18
    375                                    1.10_wp,  1.10E-2_wp          & ! 19
     355    REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ &
     356                                   0.25_wp,  0.25E-2_wp,                    & !  1
     357                                   0.20_wp,  0.20E-2_wp,                    & !  2
     358                                   2.00_wp,     2.00_wp,                    & !  3
     359                                   2.00_wp,     2.00_wp,                    & !  4
     360                                   2.00_wp,     2.00_wp,                    & !  5
     361                                   2.00_wp,     2.00_wp,                    & !  6
     362                                   0.47_wp,  0.47E-2_wp,                    & !  7
     363                                  0.013_wp, 0.013E-2_wp,                    & !  8
     364                                  0.034_wp, 0.034E-2_wp,                    & !  9
     365                                    0.5_wp,  0.50E-2_wp,                    & ! 10
     366                                   0.17_wp,  0.17E-2_wp,                    & ! 11
     367                                 1.3E-3_wp,   1.3E-4_wp,                    & ! 12
     368                                   0.83_wp,  0.83E-2_wp,                    & ! 13
     369                                   0.00_wp,  0.00E-2_wp,                    & ! 14
     370                                   0.00_wp,  0.00E-2_wp,                    & ! 15
     371                                   0.10_wp,  0.10E-2_wp,                    & ! 16
     372                                   0.25_wp,  0.25E-2_wp,                    & ! 17
     373                                   2.00_wp,  2.00E-2_wp,                    & ! 18
     374                                   1.10_wp,  1.10E-2_wp                     & ! 19
    376375                                 /), (/ 2, 19 /) )
    377376
    378377!
    379378!-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
    380     REAL(wp), DIMENSION(0:2,1:19) :: surface_pars = RESHAPE( (/          &
    381                                       10.0_wp,       10.0_wp, 0.05_wp, & !  1
    382                                       10.0_wp,       10.0_wp, 0.05_wp, & !  2
    383                                       20.0_wp,       15.0_wp, 0.03_wp, & !  3
    384                                       20.0_wp,       15.0_wp, 0.03_wp, & !  4
    385                                       20.0_wp,       15.0_wp, 0.03_wp, & !  5
    386                                       20.0_wp,       15.0_wp, 0.03_wp, & !  6
    387                                       10.0_wp,       10.0_wp, 0.05_wp, & !  7
    388                                       15.0_wp,       15.0_wp, 0.00_wp, & !  8
    389                                       10.0_wp,       10.0_wp, 0.05_wp, & !  9
    390                                       10.0_wp,       10.0_wp, 0.05_wp, & ! 10
    391                                       10.0_wp,       10.0_wp, 0.05_wp, & ! 11
    392                                       58.0_wp,       58.0_wp, 0.00_wp, & ! 12
    393                                       10.0_wp,       10.0_wp, 0.05_wp, & ! 13
    394                                     1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 14
    395                                     1.0E20_wp,     1.0E20_wp, 0.00_wp, & ! 15
    396                                       10.0_wp,       10.0_wp, 0.05_wp, & ! 16
    397                                       10.0_wp,       10.0_wp, 0.05_wp, & ! 17
    398                                       20.0_wp,       15.0_wp, 0.03_wp, & ! 18
    399                                       20.0_wp,       15.0_wp, 0.03_wp  & ! 19
     379    REAL(wp), DIMENSION(0:2,1:19), PARAMETER :: surface_pars = RESHAPE( (/ &
     380                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  1
     381                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  2
     382                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  3
     383                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  4
     384                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  5
     385                                      20.0_wp,       15.0_wp, 0.03_wp,     & !  6
     386                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  7
     387                                      15.0_wp,       15.0_wp, 0.00_wp,     & !  8
     388                                      10.0_wp,       10.0_wp, 0.05_wp,     & !  9
     389                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
     390                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
     391                                      58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
     392                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
     393                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 14
     394                                    1.0E20_wp,     1.0E20_wp, 0.00_wp,     & ! 15
     395                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
     396                                      10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
     397                                      20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
     398                                      20.0_wp,       15.0_wp, 0.03_wp      & ! 19
    400399                                      /), (/ 3, 19 /) )
    401400
    402401!
    403402!-- Root distribution (sum = 1)  level 1, level 2, level 3, level 4,
    404     REAL(wp), DIMENSION(0:3,1:19) :: root_distribution = RESHAPE( (/ &
    405                                  0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & !  1
    406                                  0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & !  2
    407                                  0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & !  3
    408                                  0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & !  4
    409                                  0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & !  5
    410                                  0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & !  6
    411                                  0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & !  7
    412                                  1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  8
    413                                  0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & !  9
    414                                  0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 10
    415                                  0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 11
    416                                  0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 12
    417                                  0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 13
    418                                  0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 14
    419                                  0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 15
    420                                  0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
    421                                  0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17
    422                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18
    423                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp  & ! 19
     403    REAL(wp), DIMENSION(0:3,1:19), PARAMETER :: root_distribution = RESHAPE( (/ &
     404                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  1
     405                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  2
     406                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  3
     407                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  4
     408                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  5
     409                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  6
     410                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  7
     411                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  8
     412                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & !  9
     413                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 10
     414                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 11
     415                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 12
     416                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 13
     417                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 14
     418                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 15
     419                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
     420                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 17
     421                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 18
     422                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 19
    424423                                 /), (/ 4, 19 /) )
    425424
    … …  
    429428!
    430429!-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
    431     REAL(wp), DIMENSION(0:3,1:7) :: soil_pars = RESHAPE( (/                &
     430    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
    432431                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
    433432                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
    … …  
    441440!
    442441!-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
    443     REAL(wp), DIMENSION(0:3,1:7) :: m_soil_pars = RESHAPE( (/            &
     442    REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
    444443                                 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
    445444                                 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
    … …  
    490489!-- Public prognostic variables
    491490    PUBLIC m_liq_eb, m_liq_eb_av, m_soil, m_soil_av, t_soil, t_soil_av
    492 
    493491
    494492    INTERFACE init_lsm
    … …  
    573571       ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
    574572       ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
     573       ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
    575574       ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
    576575       ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
    … …  
    583582#if ! defined( __nopointer )
    584583!
    585 !-- Initial assignment of the pointers
     584!--    Initial assignment of the pointers
    586585       t_soil    => t_soil_1;    t_soil_p    => t_soil_2
    587586       t_surface => t_surface_1; t_surface_p => t_surface_2
    … …  
    890889!
    891890!--    Add timeseries for land surface model
     891
    892892       dots_label(dots_num+1) = "ghf_eb"
    893893       dots_label(dots_num+2) = "shf_eb"
    … …  
    904904       dots_soil = dots_num + 1
    905905       dots_num  = dots_num + 8
    906 
    907906
    908907       RETURN
    … …  
    957956       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    958957
     958       
    959959
    960960       DO i = nxlg, nxrg
    … …  
    984984!--          f1: correction for incoming shortwave radiation (stomata close at
    985985!--          night)
    986              IF ( irad_scheme /= 0 )  THEN
    987                 f1 = MIN(1.0_wp, ( 0.004_wp * rad_sw_in(j,i) + 0.05_wp ) /     &
    988                               (0.81_wp * (0.004_wp * rad_sw_in(j,i) + 1.0_wp) ))
     986             IF ( radiation_scheme /= 'constant' )  THEN
     987                f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /&
     988                              (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)        &
     989                               + 1.0_wp)) )
    989990             ELSE
    990991                f1 = 1.0_wp
    … …  
    11121113!
    11131114!--          Add LW up so that it can be removed in prognostic equation
    1114              rad_net(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
     1115             rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4
    11151116
    11161117             IF ( humidity )  THEN
    … …  
    11181119!
    11191120!--             Numerator of the prognostic equation
    1120                 coef_1 = rad_net(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
     1121                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
    11211122                         + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s   &
    11221123                         + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
    … …  
    11321133!
    11331134!--             Numerator of the prognostic equation
    1134                 coef_1 = rad_net(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
     1135                coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4&
    11351136                         + f_shf / exn * T_1 + lambda_surface                  &
    11361137                         * t_soil(nzb_soil,j,i)
    … …  
    11751176!
    11761177!--          Calculate fluxes
    1177              rad_net(j,i)   = rad_net(j,i) + 3.0_wp * sigma_sb                 &
     1178             rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb                 &
    11781179                              * t_surface(j,i)**4 - 4.0_wp * sigma_sb          &
    11791180                              * t_surface(j,i)**3 * t_surface_p(j,i)
    … …  
    16211622          CASE ( 0 )
    16221623
    1623              t_surface = t_surface_p
    1624              t_soil    = t_soil_p
    1625              IF ( humidity )  THEN
    1626                 m_soil    = m_soil_p
    1627                 m_liq_eb  = m_liq_eb_p
    1628              ENDIF
    1629 
    1630 
    1631           CASE ( 1 )
    1632 
    16331624             t_surface  => t_surface_1; t_surface_p  => t_surface_2
    16341625             t_soil     => t_soil_1;    t_soil_p     => t_soil_2
    … …  
    16381629             ENDIF
    16391630
     1631
     1632          CASE ( 1 )
     1633
     1634             t_surface  => t_surface_2; t_surface_p  => t_surface_1
     1635             t_soil     => t_soil_2;    t_soil_p     => t_soil_1
     1636             IF ( humidity )  THEN
     1637                m_soil    => m_soil_2;   m_soil_p    => m_soil_1
     1638                m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
     1639             ENDIF
     1640
    16401641       END SELECT
    16411642#endif

  • palm/trunk/SOURCE/netcdf.f90

    r1552 r1585  
    148148    USE profil_parameter,                                                       &
    149149        ONLY:  crmax, cross_profiles, dopr_index,profile_columns, profile_rows
     150
     151    USE radiation_model_mod,                                                    &
     152        ONLY:  rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
    150153
    151154    USE spectrum,                                                               &
    … …  
    523526                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',&
    524527                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',    &
    525                        'rho', 's', 'sa', 'vpt' )
     528                       'rho', 's', 'sa', 'vpt', 'rad_lw_in', 'rad_lw_out', &
     529                       'rad_sw_in', 'rad_sw_out' )
    526530
    527531                   grid_x = 'x'
    … …  
    11031107                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
    11041108                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
    1105                        's', 'sa', 'vpt' )
     1109                       's', 'sa', 'vpt' , 'rad_lw_in', 'rad_lw_out',          &
     1110                       'rad_sw_in', 'rad_sw_out' )
    11061111
    11071112                   grid_x = 'x'
    … …  
    17561761                             'prr_xy', 'pt_xy', 'q_xy', 'qc_xy', 'ql_xy',        &
    17571762                             'ql_c_xy', 'ql_v_xy', 'ql_vp_xy', 'qr_xy', 'qv_xy', &
    1758                              'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
     1763                             'rho_xy', 's_xy', 'sa_xy', 'vpt_xy', 'rad_lw_in_xy',&
     1764                             'rad_lw_out_xy', 'rad_sw_in_xy', 'rad_sw_out_xy' )
    17591765
    17601766                         grid_x = 'x'
    … …  
    24392445                          'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
    24402446                          'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
    2441                           'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
     2447                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' , 'rad_lw_in_xz',&
     2448                             'rad_lw_out_xz', 'rad_sw_in_xz', 'rad_sw_out_xz' )
    24422449
    24432450                      grid_x = 'x'
    … …  
    31153122                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',        &
    31163123                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz', &
    3117                           'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
     3124                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz', 'rad_lw_in_yz',&
     3125                             'rad_lw_out_yz', 'rad_sw_in_yz', 'rad_sw_out_yz' )
    31183126
    31193127                      grid_x = 'x'

  • palm/trunk/SOURCE/package_parin.f90

    r1576 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added several radiation_par parameters
    2323!
    2424! Former revisions:
    … …  
    144144
    145145    USE radiation_model_mod,                                                   &
    146         ONLY: albedo, day_init, dt_radiation, lambda, net_radiation, radiation,&
    147               radiation_scheme, time_radiation, time_utc_init
     146        ONLY: albedo, albedo_type, albedo_lw_dif, albedo_lw_dir, albedo_sw_dif,&
     147              albedo_sw_dir, day_init, constant_albedo, dt_radiation, lambda,  &
     148              lw_radiation, net_radiation, radiation, radiation_scheme,        &
     149              sw_radiation, time_utc_init
    148150               
    149151
    … …  
    221223                                  write_particle_statistics
    222224
    223     NAMELIST /radiation_par/      albedo, day_init, dt_radiation, lambda,      &
    224                                   net_radiation, radiation_scheme, time_utc_init
     225    NAMELIST /radiation_par/      albedo, albedo_type, albedo_lw_dir,          &
     226                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
     227                                  constant_albedo, day_init, dt_radiation,     &
     228                                  lambda, lw_radiation, net_radiation,         &
     229                                  radiation_scheme, sw_radiation, time_utc_init
    225230
    226231    NAMELIST /spectra_par/        averaging_interval_sp, comp_spectra_level,   &

  • palm/trunk/SOURCE/prognostic_equations.f90

    r1518 r1585  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Added call for temperature tendency calculation due to radiative flux divergence
    2323!
    2424! Former revisions:
    … …  
    277277        ONLY:  production_e, production_e_acc
    278278
     279    USE radiation_model_mod,                                                   &
     280        ONLY:  radiation, radiation_scheme, radiation_tendency
     281
    279282    USE statistics,                                                            &
    280283        ONLY:  hom
    … …  
    593596                CALL subsidence( i, j, tend, pt, pt_init, 2 )
    594597             ENDIF
     598
     599!
     600!--          If required, add tendency due to radiative heating/cooling
     601             IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
     602                CALL radiation_tendency ( i, j, tend )
     603             ENDIF
     604
    595605
    596606             CALL user_actions( i, j, 'pt-tendency' )
    … …  
    12461256            .NOT. use_subsidence_tendencies )  THEN
    12471257          CALL subsidence( tend, pt, pt_init, 2 )
     1258       ENDIF
     1259
     1260!
     1261!--    If required, add tendency due to radiative heating/cooling
     1262       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
     1263          CALL radiation_tendency ( tend )
    12481264       ENDIF
    12491265

  • palm/trunk/SOURCE/radiation_model.f90

    r1572 r1585  
    1515! PALM. If not, see <http://www.gnu.org/licenses/>.
    1616!
    17 ! Copyright 1997-2014 Leibniz Universitaet Hannover
     17! Copyright 1997-2015 Leibniz Universitaet Hannover
    1818!--------------------------------------------------------------------------------!
    1919!
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Added support for RRTMG
    2323!
    2424! Former revisions:
    … …  
    4040! Description:
    4141! ------------
    42 ! Radiation model(s), to be used e.g. with the land surface scheme
     42! Radiation models and interfaces
     43! To do: move variable definitions used in init_radiation only to the subroutine
     44! as they are no longer required after initialization.
     45! Note that many variables have a leading dummy dimension (0:0) in order to
     46! match the assume-size shape expected by the RRTMG model
     47!
     48! To do:
     49! - Output of full column vertical profiles used in RRTMG
     50! - Output of other rrtm arrays (such as volume mixing ratios)
     51! - Adapt for use with topography
     52! - Remove 3D dummy arrays (such as clear-sky output)
    4353!------------------------------------------------------------------------------!
    4454
    4555    USE arrays_3d,                                                             &
    46         ONLY: pt
     56        ONLY:  hyp, pt, q, ql, zw
     57
     58    USE cloud_parameters,                                                      &
     59        ONLY:  cp, l_v, nc_const, rho_l, sigma_gc 
     60
     61    USE constants,                                                             &
     62        ONLY:  pi
    4763
    4864    USE control_parameters,                                                    &
    49         ONLY: phi, surface_pressure, time_since_reference_point
     65        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
     66               large_scale_forcing, lsf_surf, phi, pt_surface,                 &
     67               surface_pressure, time_since_reference_point
    5068
    5169    USE indices,                                                               &
    52         ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner
     70        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
    5371
    5472    USE kinds
     73
     74    USE netcdf
    5575
    5676    USE netcdf_control,                                                        &
    5777        ONLY:  dots_label, dots_num, dots_unit
    5878
     79#if defined ( __rrtmg )
     80    USE parrrsw,                                                               &
     81        ONLY:  naerec, nbndsw
     82
     83    USE parrrtm,                                                               &
     84        ONLY:  nbndlw
     85
     86    USE rrtmg_lw_init,                                                         &
     87        ONLY:  rrtmg_lw_ini
     88
     89    USE rrtmg_sw_init,                                                         &
     90        ONLY:  rrtmg_sw_ini
     91
     92    USE rrtmg_lw_rad,                                                          &
     93        ONLY:  rrtmg_lw
     94
     95    USE rrtmg_sw_rad,                                                          &
     96        ONLY:  rrtmg_sw
     97#endif
     98
     99
    59100
    60101    IMPLICIT NONE
    61102
    62     CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtm'
    63 
    64     INTEGER(iwp) :: i, j, k
    65 
    66 
    67     INTEGER(iwp) :: day_init     = 172,  & !: day of the year at model start (21/06)
    68                     dots_rad     = 0,    & !: starting index for timeseries output
    69                     irad_scheme  = 0
    70 
    71     LOGICAL :: radiation = .FALSE.  !: flag parameter indicating wheather the radiation model is used
    72 
    73     REAL(wp), PARAMETER :: sw_0 = 1368.0_wp, &    !: solar constant 
    74                            pi = 3.14159265358979323_wp, &
    75                            sigma_sb  = 5.67E-8_wp !: Stefan-Boltzmann constant
     103    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
     104
     105!
     106!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
     107    CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/    &
     108                                   'user defined',                         & !  0
     109                                   'ocean',                                & !  1
     110                                   'mixed farming, tall grassland',        & !  2
     111                                   'tall/medium grassland',                & !  3
     112                                   'evergreen shrubland',                  & !  4
     113                                   'short grassland/meadow/shrubland',     & !  5
     114                                   'evergreen needleleaf forest',          & !  6
     115                                   'mixed deciduous evergreen forest',     & !  7
     116                                   'deciduous forest',                     & !  8
     117                                   'tropical evergreen broadleaved forest',& !  9
     118                                   'medium/tall grassland/woodland',       & ! 10
     119                                   'desert, sandy',                        & ! 11
     120                                   'desert, rocky',                        & ! 12
     121                                   'tundra',                               & ! 13
     122                                   'landice',                              & ! 14
     123                                   'sea ice',                              & ! 15
     124                                   'snow'                                  & ! 16
     125                                                         /)
     126
     127    INTEGER(iwp) :: albedo_type  = 5,    & !: Albedo surface type (default: short grassland)
     128                    day,                 & !: current day of the year
     129                    day_init     = 172,  & !: day of the year at model start (21/06)
     130                    dots_rad     = 0       !: starting index for timeseries output
     131
     132
     133
     134
     135
     136
     137    LOGICAL ::  constant_albedo = .FALSE., & !: flag parameter indicating whether the albedo may change depending on zenith
     138                lw_radiation = .TRUE.,     & !: flag parameter indicing whether longwave radiation shall be calculated
     139                radiation = .FALSE.,       & !: flag parameter indicating whether the radiation model is used
     140                sun_up    = .TRUE.,        & !: flag parameter indicating whether the sun is up or down
     141                sw_radiation = .TRUE.        !: flag parameter indicing whether shortwave radiation shall be calculated
     142
     143    REAL(wp), PARAMETER :: sigma_sb       = 5.67E-8_wp, & !: Stefan-Boltzmann constant
     144                           solar_constant = 1368.0_wp     !: solar constant at top of atmosphere
    76145 
    77     REAL(wp) :: albedo = 0.2_wp,             & !: NAMELIST alpha
    78                 dt_radiation = 0.0_wp,       & !: radiation model timestep
    79                 exn,                         & !: Exner function
    80                 lon = 0.0_wp,                & !: longitude in radians
    81                 lat = 0.0_wp,                & !: latitude in radians
    82                 decl_1,                      & !: declination coef. 1
    83                 decl_2,                      & !: declination coef. 2
    84                 decl_3,                      & !: declination coef. 3
    85                 time_utc,                    & !: current time in UTC
    86                 time_utc_init = 43200.0_wp,  & !: UTC time at model start (noon)
    87                 day,                         & !: current day of the year
    88                 lambda = 0.0_wp,             & !: longitude in degrees
    89                 declination,                 & !: solar declination angle
    90                 net_radiation = 0.0_wp,      & !: net radiation at surface
    91                 hour_angle,                  & !: solar hour angle
    92                 time_radiation = 0.0_wp,     &
    93                 zenith,                      & !: solar zenith angle
    94                 sky_trans                      !: sky transmissivity
     146    REAL(wp) :: albedo = 9999999.9_wp,        & !: NAMELIST alpha
     147                albedo_lw_dif = 9999999.9_wp, & !: NAMELIST aldif
     148                albedo_lw_dir = 9999999.9_wp, & !: NAMELIST aldir
     149                albedo_sw_dif = 9999999.9_wp, & !: NAMELIST asdif
     150                albedo_sw_dir = 9999999.9_wp, & !: NAMELIST asdir
     151                dt_radiation = 0.0_wp,        & !: radiation model timestep
     152                emissivity = 0.95_wp,         & !: NAMELIST surface emissivity
     153                exn,                          & !: Exner function
     154                lon = 0.0_wp,                 & !: longitude in radians
     155                lat = 0.0_wp,                 & !: latitude in radians
     156                decl_1,                       & !: declination coef. 1
     157                decl_2,                       & !: declination coef. 2
     158                decl_3,                       & !: declination coef. 3
     159                time_utc,                     & !: current time in UTC
     160                time_utc_init = 43200.0_wp,   & !: UTC time at model start (noon)
     161                lambda = 0.0_wp,              & !: longitude in degrees
     162                net_radiation = 0.0_wp,       & !: net radiation at surface
     163                time_radiation = 0.0_wp,      & !: time since last call of radiation code
     164                sky_trans                       !: sky transmissivity
     165
     166
     167    REAL(wp), DIMENSION(0:0) ::  zenith        !: solar zenith angle
    95168
    96169    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
    97                 alpha,                       & !: surface albedo
     170                alpha,                       & !: surface broadband albedo (used for clear-sky scheme)
    98171                rad_net,                     & !: net radiation at the surface
    99                 rad_net_av,                  & !: average of rad_net
    100                 rad_lw_in,                   & !: incoming longwave radiation
    101                 rad_lw_out,                  & !: outgoing longwave radiation
    102                 rad_sw_in,                   & !: incoming shortwave radiation
    103                 rad_sw_in_av,                & !: average of rad_sw_in
    104                 rad_sw_out                     !: outgoing shortwave radiation
    105 
     172                rad_net_av                     !: average of rad_net
     173
     174!
     175!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
     176!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
     177    REAL(wp), DIMENSION(0:2,1:15), PARAMETER :: albedo_pars = RESHAPE( (/&
     178                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
     179                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
     180                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
     181                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
     182                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
     183                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
     184                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
     185                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
     186                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
     187                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
     188                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
     189                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
     190                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
     191                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
     192                                   0.95_wp, 0.70_wp, 0.82_wp             & ! 15
     193                                 /), (/ 3, 15 /) )
     194
     195    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
     196                        rad_sw_in,                     & !: incoming shortwave radiation (W/m2)
     197                        rad_sw_in_av,                  & !: average of rad_sw_in
     198                        rad_sw_out,                    & !: outgoing shortwave radiation (W/m2)
     199                        rad_sw_out_av,                 & !: average of rad_sw_out
     200                        rad_lw_in,                     & !: incoming longwave radiation (W/m2)
     201                        rad_lw_in_av,                  & !: average of rad_lw_in
     202                        rad_lw_out,                    & !: outgoing longwave radiation (W/m2)
     203                        rad_lw_out_av                    !: average of rad_lw_out
     204
     205!
     206!-- Variables and parameters used in RRTMG only
     207#if defined ( __rrtmg )
     208    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !: name of the NetCDF input file (sounding data)
     209
     210
     211!
     212!-- Flag parameters for RRTMGS (should not be changed)
     213    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !: flag for lw cloud optical properties (0,1,2)
     214                               rrtm_iceflglw = 0, & !: flag for lw ice particle specifications (0,1,2,3)
     215                               rrtm_liqflglw = 1, & !: flag for lw liquid droplet specifications
     216                               rrtm_inflgsw  = 2, & !: flag for sw cloud optical properties (0,1,2)
     217                               rrtm_iceflgsw = 0, & !: flag for sw ice particle specifications (0,1,2,3)
     218                               rrtm_liqflgsw = 1    !: flag for sw liquid droplet specifications
     219
     220!
     221!-- The following variables should be only changed with care, as this will
     222!-- require further setting of some variables, which is currently not
     223!-- implemented (aerosols, ice phase).
     224    INTEGER(iwp) :: nzt_rad,           & !: upper vertical limit for radiation calculations
     225                    rrtm_icld = 0,     & !: cloud flag (0: clear sky column, 1: cloudy column)
     226                    rrtm_iaer = 0,     & !: aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
     227                    rrtm_idrv = 0        !: longwave upward flux calculation option (0,1)
     228
     229    LOGICAL :: snd_exists = .FALSE.      !: flag parameter to check whether a user-defined input files exists
     230
     231    REAL(wp), PARAMETER :: mol_weight_air_d_wv = 1.607793_wp !: molecular weight dry air / water vapor
     232
     233    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !: hypostatic pressure from sounding data (hPa)
     234                                           q_snd,       & !: specific humidity from sounding data (kg/kg) - dummy at the moment
     235                                           rrtm_tsfc,   & !: dummy array for storing surface temperature
     236                                           t_snd          !: actual temperature from sounding data (hPa)
     237
     238    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,         & !: longwave diffuse albedo solar angle of 60°
     239                                             aldir,         & !: longwave direct albedo solar angle of 60°
     240                                             asdif,         & !: shortwave diffuse albedo solar angle of 60°
     241                                             asdir,         & !: shortwave direct albedo solar angle of 60°
     242                                             rrtm_ccl4vmr,  & !: CCL4 volume mixing ratio (g/mol)
     243                                             rrtm_cfc11vmr, & !: CFC11 volume mixing ratio (g/mol)
     244                                             rrtm_cfc12vmr, & !: CFC12 volume mixing ratio (g/mol)
     245                                             rrtm_cfc22vmr, & !: CFC22 volume mixing ratio (g/mol)
     246                                             rrtm_ch4vmr,   & !: CH4 volume mixing ratio
     247                                             rrtm_cicewp,   & !: in-cloud ice water path (g/m²)
     248                                             rrtm_cldfr,    & !: cloud fraction (0,1)
     249                                             rrtm_cliqwp,   & !: in-cloud liquid water path (g/m²)
     250                                             rrtm_co2vmr,   & !: CO2 volume mixing ratio (g/mol)
     251                                             rrtm_emis,     & !: surface emissivity (0-1)   
     252                                             rrtm_h2ovmr,   & !: H2O volume mixing ratio
     253                                             rrtm_n2ovmr,   & !: N2O volume mixing ratio
     254                                             rrtm_o2vmr,    & !: O2 volume mixing ratio
     255                                             rrtm_o3vmr,    & !: O3 volume mixing ratio
     256                                             rrtm_play,     & !: pressure layers (hPa, zu-grid)
     257                                             rrtm_plev,     & !: pressure layers (hPa, zw-grid)
     258                                             rrtm_reice,    & !: cloud ice effective radius (microns)
     259                                             rrtm_reliq,    & !: cloud water drop effective radius (microns)
     260                                             rrtm_tlay,     & !: actual temperature (K, zu-grid)
     261                                             rrtm_tlev,     & !: actual temperature (K, zw-grid)
     262                                             rrtm_lwdflx,   & !: RRTM output of incoming longwave radiation flux (W/m2)
     263                                             rrtm_lwuflx,   & !: RRTM output of outgoing longwave radiation flux (W/m2)
     264                                             rrtm_lwhr,     & !: RRTM output of longwave radiation heating rate (K/d)
     265                                             rrtm_lwuflxc,  & !: RRTM output of incoming clear sky longwave radiation flux (W/m2)
     266                                             rrtm_lwdflxc,  & !: RRTM output of outgoing clear sky longwave radiation flux (W/m2)
     267                                             rrtm_lwhrc,    & !: RRTM output of incoming longwave clear sky radiation heating rate (K/d)
     268                                             rrtm_swdflx,   & !: RRTM output of incoming shortwave radiation flux (W/m2)
     269                                             rrtm_swuflx,   & !: RRTM output of outgoing shortwave radiation flux (W/m2)
     270                                             rrtm_swhr,     & !: RRTM output of shortwave radiation heating rate (K/d)
     271                                             rrtm_swuflxc,  & !: RRTM output of incoming clear sky shortwave radiation flux (W/m2)
     272                                             rrtm_swdflxc,  & !: RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
     273                                             rrtm_swhrc       !: RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
     274
     275!
     276!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
     277    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !: incoming clear sky longwave radiation (W/m2) (not used)
     278                                                rad_lw_cs_out,  & !: outgoing clear sky longwave radiation (W/m2) (not used)
     279                                                rad_lw_cs_hr,   & !: longwave clear sky radiation heating rate (K/d) (not used)
     280                                                rad_lw_hr,      & !: longwave radiation heating rate (K/d)
     281                                                rad_sw_cs_in,   & !: incoming clear sky shortwave radiation (W/m2) (not used)
     282                                                rad_sw_cs_out,  & !: outgoing clear sky shortwave radiation (W/m2) (not used)
     283                                                rad_sw_cs_hr,   & !: shortwave clear sky radiation heating rate (K/d) (not used)
     284                                                rad_sw_hr,      & !: shortwave radiation heating rate (K/d)
     285                                                rrtm_aldif,     & !: surface albedo for longwave diffuse radiation
     286                                                rrtm_aldir,     & !: surface albedo for longwave direct radiation
     287                                                rrtm_asdif,     & !: surface albedo for shortwave diffuse radiation
     288                                                rrtm_asdir,     & !: surface albedo for shortwave direct radiation
     289                                                rrtm_lw_tauaer, & !: lw aerosol optical depth
     290                                                rrtm_lw_taucld, & !: lw in-cloud optical depth
     291                                                rrtm_sw_taucld, & !: sw in-cloud optical depth
     292                                                rrtm_sw_ssacld, & !: sw in-cloud single scattering albedo
     293                                                rrtm_sw_asmcld, & !: sw in-cloud asymmetry parameter
     294                                                rrtm_sw_fsfcld, & !: sw in-cloud forward scattering fraction
     295                                                rrtm_sw_tauaer, & !: sw aerosol optical depth
     296                                                rrtm_sw_ssaaer, & !: sw aerosol single scattering albedo
     297                                                rrtm_sw_asmaer, & !: sw aerosol asymmetry parameter
     298                                                rrtm_sw_ecaer     !: sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
     299#endif
    106300
    107301    INTERFACE init_radiation
    … …  
    113307    END INTERFACE radiation_clearsky
    114308
    115     INTERFACE radiation_constant
    116        MODULE PROCEDURE radiation_constant
    117     END INTERFACE radiation_constant
    118 
    119     INTERFACE radiation_rrtm
    120        MODULE PROCEDURE radiation_rrtm
    121     END INTERFACE radiation_rrtm
    122 
     309    INTERFACE radiation_rrtmg
     310       MODULE PROCEDURE radiation_rrtmg
     311    END INTERFACE radiation_rrtmg
     312
     313    INTERFACE radiation_tendency
     314       MODULE PROCEDURE radiation_tendency
     315       MODULE PROCEDURE radiation_tendency_ij
     316    END INTERFACE radiation_tendency
    123317
    124318    SAVE
    … …  
    126320    PRIVATE
    127321
    128     PUBLIC albedo, day_init, dots_rad, dt_radiation, init_radiation,           &
    129            irad_scheme, lambda, net_radiation, rad_net, rad_net_av, radiation, &
    130            radiation_clearsky, radiation_constant, radiation_rrtm,             &
    131            radiation_scheme, rad_sw_in, rad_sw_in_av, sigma_sb,                &
    132            time_radiation, time_utc_init
    133 
    134 
     322    PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
     323           albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad,  &
     324           dt_radiation, init_radiation, lambda, lw_radiation, net_radiation,  &
     325           rad_net, rad_net_av, radiation, radiation_clearsky,                 &
     326           radiation_rrtmg, radiation_scheme, radiation_tendency, rad_lw_in,   &
     327           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_sw_in, rad_sw_in_av,   &
     328           rad_sw_out, rad_sw_out_av, sigma_sb, sw_radiation, time_radiation,  &
     329           time_utc_init
     330
     331#if defined ( __rrtmg )
     332    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
     333#endif
    135334
    136335 CONTAINS
    … …  
    143342    SUBROUTINE init_radiation
    144343   
    145 
    146344       IMPLICIT NONE
    147345
    148        ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
    149        ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
    150        ALLOCATE ( rad_lw_in(nysg:nyng,nxlg:nxrg) )
    151        ALLOCATE ( rad_lw_out(nysg:nyng,nxlg:nxrg) )
    152        ALLOCATE ( rad_sw_in(nysg:nyng,nxlg:nxrg) )
    153        ALLOCATE ( rad_sw_out(nysg:nyng,nxlg:nxrg) )
    154 
    155        rad_sw_in  = 0.0_wp
    156        rad_sw_out = 0.0_wp
    157        rad_lw_in  = 0.0_wp
    158        rad_lw_out = 0.0_wp
    159        rad_net    = 0.0_wp
    160 
    161        alpha = albedo
     346!
     347!--    Allocate array for storing the surface net radiation
     348       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
     349          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
     350          rad_net = 0.0_wp
     351       ENDIF
    162352
    163353!
    164354!--    Fix net radiation in case of radiation_scheme = 'constant'
    165        IF ( irad_scheme == 0 )  THEN
     355       IF ( radiation_scheme == 'constant' )  THEN
    166356          rad_net = net_radiation
    167 !
    168 !--    Calculate radiation scheme constants
    169        ELSEIF ( irad_scheme == 1 .OR. irad_scheme == 2 )  THEN
     357          radiation = .FALSE.
     358!
     359!--    Calculate orbital constants
     360       ELSE
    170361          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
    171362          decl_2 = 2.0_wp * pi / 365.0_wp
    172363          decl_3 = decl_2 * 81.0_wp
    173 !
    174 !--       Calculate latitude and longitude angles (lat and lon, respectively)
    175           lat = phi * pi / 180.0_wp
    176           lon = lambda * pi / 180.0_wp
     364          lat    = phi * pi / 180.0_wp
     365          lon    = lambda * pi / 180.0_wp
    177366       ENDIF
     367
     368
     369       IF ( radiation_scheme == 'clear-sky' )  THEN
     370
     371          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
     372
     373          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
     374             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
     375          ENDIF
     376          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
     377             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
     378          ENDIF
     379
     380          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
     381             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     382          ENDIF
     383          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
     384             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     385          ENDIF
     386
     387          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
     388             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
     389          ENDIF
     390          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
     391             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
     392          ENDIF
     393
     394          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
     395             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     396          ENDIF
     397          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
     398             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     399          ENDIF
     400
     401          rad_sw_in  = 0.0_wp
     402          rad_sw_out = 0.0_wp
     403          rad_lw_in  = 0.0_wp
     404          rad_lw_out = 0.0_wp
     405
     406!
     407!--       Overwrite albedo if manually set in parameter file
     408          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
     409             albedo = albedo_pars(2,albedo_type)
     410          ENDIF
     411   
     412          alpha = albedo
     413 
     414!
     415!--    Initialization actions for RRTMG
     416       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     417#if defined ( __rrtmg )
     418!
     419!--       Allocate albedos
     420          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
     421          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
     422          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
     423          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
     424          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
     425          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
     426          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
     427          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
     428
     429          IF ( albedo_type /= 0 )  THEN
     430             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
     431                albedo_lw_dif = albedo_pars(0,albedo_type)
     432                albedo_lw_dir = albedo_lw_dif
     433             ENDIF
     434             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
     435                albedo_sw_dif = albedo_pars(1,albedo_type)
     436                albedo_sw_dir = albedo_sw_dif
     437             ENDIF
     438          ENDIF
     439
     440          aldif(:,:) = albedo_lw_dif
     441          aldir(:,:) = albedo_lw_dir
     442          asdif(:,:) = albedo_sw_dif
     443          asdir(:,:) = albedo_sw_dir
     444!
     445!--       Calculate initial values of current (cosine of) the zenith angle and
     446!--       whether the sun is up
     447          CALL calc_zenith     
     448!
     449!--       Calculate initial surface albedo
     450          IF ( .NOT. constant_albedo )  THEN
     451             CALL calc_albedo
     452          ELSE
     453             rrtm_aldif(0,:,:) = aldif(:,:)
     454             rrtm_aldir(0,:,:) = aldir(:,:)
     455             rrtm_asdif(0,:,:) = asdif(:,:)
     456             rrtm_asdir(0,:,:) = asdir(:,:)   
     457          ENDIF
     458
     459!
     460!--       Allocate surface emissivity
     461          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
     462          rrtm_emis = emissivity
     463
     464!
     465!--       Allocate 3d arrays of radiative fluxes and heating rates
     466          ALLOCATE ( rad_sw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     467          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     468          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     469          ALLOCATE ( rad_sw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     470
     471          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
     472             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     473             rad_sw_in = 0.0_wp
     474
     475          ENDIF
     476
     477          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
     478             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     479             rad_sw_out = 0.0_wp
     480          ENDIF
     481
     482          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
     483             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     484          ENDIF
     485
     486          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
     487             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     488          ENDIF
     489
     490          ALLOCATE ( rad_lw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     491          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     492          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     493          ALLOCATE ( rad_lw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     494
     495          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
     496             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     497             rad_lw_in     = 0.0_wp
     498          ENDIF
     499
     500          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
     501             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     502          ENDIF
     503
     504          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
     505             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     506            rad_lw_out    = 0.0_wp
     507          ENDIF
     508
     509          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
     510             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     511          ENDIF
     512
     513          rad_sw_hr     = 0.0_wp
     514          rad_sw_cs_in  = 0.0_wp
     515          rad_sw_cs_out = 0.0_wp
     516          rad_sw_cs_hr  = 0.0_wp
     517
     518          rad_lw_hr     = 0.0_wp
     519          rad_lw_cs_in  = 0.0_wp
     520          rad_lw_cs_out = 0.0_wp
     521          rad_lw_cs_hr  = 0.0_wp
     522
     523!
     524!--       Allocate dummy array for storing surface temperature
     525          ALLOCATE ( rrtm_tsfc(1) )
     526
     527!
     528!--       Initialize RRTMG
     529          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
     530          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
     531
     532!
     533!--       Set input files for RRTMG
     534          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
     535          IF ( .NOT. snd_exists )  THEN
     536             rrtm_input_file = "rrtmg_lw.nc"
     537          ENDIF
     538
     539!
     540!--       Read vertical layers for RRTMG from sounding data
     541!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
     542!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
     543!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
     544          CALL read_sounding_data
     545
     546!
     547!--       Read trace gas profiles from file. This routine provides
     548!--       the rrtm_ arrays (1:nzt_rad+1)
     549          CALL read_trace_gas_data
     550#endif
     551       ENDIF
     552
     553!
     554!--    Perform user actions if required
     555       CALL user_init_radiation
     556
     557
    178558!
    179559!--    Add timeseries for radiation model
     560       dots_rad = dots_num + 1
    180561       dots_label(dots_num+1) = "rad_net"
    181        dots_label(dots_num+2) = "rad_sw_in"
    182        dots_unit(dots_num+1:dots_num+2) = "W/m2"
    183 
    184        dots_rad  = dots_num + 1
    185        dots_num  = dots_num + 2
     562       dots_label(dots_num+2) = "rad_lw_in"
     563       dots_label(dots_num+3) = "rad_lw_out"
     564       dots_label(dots_num+4) = "rad_sw_in"
     565       dots_label(dots_num+5) = "rad_sw_out"
     566       dots_unit(dots_num+1:dots_num+5) = "W/m2"
     567       dots_num  = dots_num + 5
     568
     569!
     570!--    Output of albedos is only required for RRTMG
     571       IF ( radiation_scheme == 'rrtmg' )  THEN
     572          dots_label(dots_num+1) = "rrtm_aldif"
     573          dots_label(dots_num+2) = "rrtm_aldir"
     574          dots_label(dots_num+3) = "rrtm_asdif"
     575          dots_label(dots_num+4) = "rrtm_asdir"
     576          dots_unit(dots_num+1:dots_num+4) = ""
     577          dots_num  = dots_num + 4
     578       ENDIF
     579
     580!
     581!--    Calculate radiative fluxes at model start
     582       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     583          IF ( radiation_scheme == 'clear-sky' )  THEN
     584              CALL radiation_clearsky
     585          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     586             CALL radiation_rrtmg
     587          ENDIF
     588       ENDIF
    186589
    187590       RETURN
    … …  
    196599!------------------------------------------------------------------------------!
    197600    SUBROUTINE radiation_clearsky
    198    
     601
     602       USE indices,                                                            &
     603           ONLY:  nbgp
    199604
    200605       IMPLICIT NONE
    201606
     607       INTEGER(iwp) :: i, j, k
     608
     609!
     610!--    Calculate current zenith angle
     611       CALL calc_zenith
     612
     613!
     614!--    Calculate sky transmissivity
     615       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
     616
     617!
     618!--    Calculate value of the Exner function
     619       exn = (surface_pressure / 1000.0_wp )**0.286_wp
     620!
     621!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
     622!--    point
     623       DO i = nxl, nxr
     624          DO j = nys, nyn
     625             k = nzb_s_inner(j,i)
     626             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
     627             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
     628             rad_lw_out(0,j,i) = sigma_sb * (pt(k,j,i) * exn)**4
     629             rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
     630             rad_net(j,i)      = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)          &
     631                                 + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
     632
     633          ENDDO
     634       ENDDO
     635
     636       CALL exchange_horiz_2d( rad_lw_in,  nbgp )
     637       CALL exchange_horiz_2d( rad_lw_out, nbgp )
     638       CALL exchange_horiz_2d( rad_sw_in,  nbgp )
     639       CALL exchange_horiz_2d( rad_sw_out, nbgp )
     640       CALL exchange_horiz_2d( rad_net,    nbgp )
     641
     642       RETURN
     643
     644    END SUBROUTINE radiation_clearsky
     645
     646
     647!------------------------------------------------------------------------------!
     648! Description:
     649! ------------
     650!-- Implementation of the RRTMG radiation_scheme
     651!------------------------------------------------------------------------------!
     652    SUBROUTINE radiation_rrtmg
     653
     654       USE indices,                                                            &
     655           ONLY:  nbgp
     656
     657       USE particle_attributes,                                                &
     658           ONLY:  grid_particles, number_of_particles, particles,              &
     659                  particle_advection_start, prt_count
     660
     661       IMPLICIT NONE
     662
     663#if defined ( __rrtmg )
     664
     665       INTEGER(iwp) :: i, j, k, n          !:
     666
     667       REAL(wp)     ::  s_r2      !: weighted sum over all droplets with r^2
     668       REAL(wp)     ::  s_r3      !: weighted sum over all droplets with r^3
     669
     670!
     671!--    Calculate current (cosine of) zenith angle and whether the sun is up
     672       CALL calc_zenith     
     673!
     674!--    Calculate surface albedo
     675       IF ( .NOT. constant_albedo )  THEN
     676          CALL calc_albedo
     677       ENDIF
     678
     679!
     680!--    Prepare input data for RRTMG
     681
     682!
     683!--    In case of large scale forcing with surface data, calculate new pressure
     684!--    profile. nzt_rad might be modified by these calls and all required arrays
     685!--    will then be re-allocated
     686       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
     687          CALL read_sounding_data
     688          CALL read_trace_gas_data
     689       ENDIF
     690!
     691!--    Loop over all grid points
     692       DO i = nxl, nxr
     693          DO j = nys, nyn
     694
     695!
     696!--          Prepare profiles of temperature and H2O volume mixing ratio
     697             rrtm_tlev(0,nzb+1) = pt(nzb,j,i)
     698
     699             DO k = nzb+1, nzt+1
     700                rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
     701                                 )**0.286_wp + ( l_v / cp ) * ql(k,j,i)
     702                rrtm_h2ovmr(0,k) = mol_weight_air_d_wv * (q(k,j,i) - ql(k,j,i))
     703
     704             ENDDO
     705
     706!
     707!--          Avoid temperature/humidity jumps at the top of the LES domain by
     708!--          linear interpolation from nzt+2 to nzt+7
     709             DO k = nzt+2, nzt+7
     710                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
     711                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
     712                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
     713                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
     714
     715                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
     716                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
     717                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
     718                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
     719
     720             ENDDO
     721
     722!--          Linear interpolate to zw grid
     723             DO k = nzb+2, nzt+8
     724                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
     725                                   rrtm_tlay(0,k-1))                           &
     726                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
     727                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
     728             ENDDO
     729
     730
     731!
     732!--          Calculate liquid water path and cloud fraction for each column.
     733!--          Note that LWP is required in g/m² instead of kg/kg m.
     734             rrtm_cldfr  = 0.0_wp
     735             rrtm_reliq  = 0.0_wp
     736             rrtm_cliqwp = 0.0_wp
     737
     738             DO k = nzb+1, nzt+1
     739                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                      &
     740                                  (rrtm_plev(0,k) - rrtm_plev(0,k+1)) * 100.0_wp / g
     741
     742                IF ( rrtm_cliqwp(0,k) .GT. 0 )  THEN
     743                   rrtm_cldfr(0,k) = 1.0_wp
     744                   rrtm_icld = 1
     745
     746!
     747!--                Calculate cloud droplet effective radius
     748                   IF ( cloud_physics )  THEN
     749                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)          &
     750                                      / ( 4.0_wp * pi * nc_const * rho_l )     &
     751                                      )**0.33333333333333_wp                   &
     752                                     * EXP( LOG( sigma_gc )**2 )
     753
     754                   ELSEIF ( cloud_droplets )  THEN
     755                      number_of_particles = prt_count(k,j,i)
     756
     757                      IF (number_of_particles <= 0)  CYCLE
     758                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
     759                      s_r2 = 0.0_wp
     760                      s_r3 = 0.0_wp
     761
     762                      DO  n = 1, number_of_particles
     763                         IF ( particles(n)%particle_mask )  THEN
     764                            s_r2 = s_r2 + particles(n)%radius**2 * &
     765                                   particles(n)%weight_factor
     766                            s_r3 = s_r3 + particles(n)%radius**3 * &
     767                                   particles(n)%weight_factor
     768                         ENDIF
     769                      ENDDO
     770
     771                      IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
     772
     773                   ENDIF
     774
     775!
     776!--                Limit effective radius
     777                   IF ( rrtm_reliq(0,k) .GT. 0.0_wp )  THEN
     778                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
     779                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
     780                  ENDIF
     781                ELSE
     782                   rrtm_icld = 0
     783                ENDIF
     784             ENDDO
     785
     786!
     787!--          Set surface temperature
     788             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
     789
     790             IF ( lw_radiation )  THEN
     791               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
     792               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
     793               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
     794               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
     795               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
     796               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
     797               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,&
     798               rrtm_reliq      , rrtm_lw_tauaer,                               &
     799               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
     800               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc )
     801
     802                DO k = nzb, nzt+1
     803                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
     804                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
     805                ENDDO
     806
     807
     808             ENDIF
     809
     810             IF ( sw_radiation .AND. sun_up )  THEN
     811                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
     812               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
     813               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
     814               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
     815               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
     816               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
     817               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
     818               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
     819               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
     820               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
     821               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
     822               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
     823 
     824                DO k = nzb, nzt+1
     825                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
     826                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
     827                ENDDO
     828             ENDIF
     829
     830!
     831!--          Calculate surface net radiation
     832             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
     833                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
     834
     835          ENDDO
     836       ENDDO
     837
     838       CALL exchange_horiz( rad_lw_in,  nbgp )
     839       CALL exchange_horiz( rad_lw_out, nbgp )
     840       CALL exchange_horiz( rad_sw_in,  nbgp )
     841       CALL exchange_horiz( rad_sw_out, nbgp )
     842       CALL exchange_horiz_2d( rad_net, nbgp )
     843#endif
     844
     845    END SUBROUTINE radiation_rrtmg
     846
     847
     848!------------------------------------------------------------------------------!
     849! Description:
     850! ------------
     851!-- Calculate the cosine of the zenith angle (variable is called zenith)
     852!------------------------------------------------------------------------------!
     853    SUBROUTINE calc_zenith
     854
     855       IMPLICIT NONE
     856
     857       REAL(wp) ::  declination,  & !: solar declination angle
     858                    hour_angle      !: solar hour angle
    202859!
    203860!--    Calculate current day and time based on the initial values and simulation
    204861!--    time
    205        day = day_init + FLOOR( (time_utc_init + time_since_reference_point)    &
    206                                / 86400.0_wp )
     862       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
     863                               / 86400.0_wp ), KIND=iwp)
    207864       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
    208865
    … …  
    210867!
    211868!--    Calculate solar declination and hour angle   
    212        declination = ASIN( decl_1 * SIN(decl_2 * day - decl_3) )
     869       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
    213870       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
    214871
    215872!
    216873!--    Calculate zenith angle
    217        zenith = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)       &
     874       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
    218875                                            * COS(hour_angle)
    219        zenith = MAX(0.0_wp,zenith)
    220 
    221 !
    222 !--    Calculate sky transmissivity
    223        sky_trans = 0.6_wp + 0.2_wp * zenith
    224 
    225 !
    226 !--    Calculate value of the Exner function
    227        exn = (surface_pressure / 1000.0_wp )**0.286_wp
    228 
    229 !
    230 !--    Calculate radiation fluxes and net radiation (rad_net) for each grid
    231 !--    point
    232        DO i = nxlg, nxrg
    233           DO j = nysg, nyng
    234 
    235              k = nzb_s_inner(j,i)
    236              rad_sw_in(j,i)  = sw_0 * sky_trans * zenith
    237              rad_sw_out(j,i) = - alpha(j,i) * rad_sw_in(j,i)
    238              rad_lw_out(j,i) = - sigma_sb * (pt(k,j,i) * exn)**4
    239              rad_lw_in(j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
    240              rad_net(j,i)    = rad_sw_in(j,i) + rad_sw_out(j,i)                &
    241                                 + rad_lw_in(j,i) + rad_lw_out(j,i)
    242 
     876       zenith(0) = MAX(0.0_wp,zenith(0))
     877
     878!
     879!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
     880       IF ( zenith(0) .GT. 0.0_wp )  THEN
     881          sun_up = .TRUE.
     882       ELSE
     883          sun_up = .FALSE.
     884       END IF
     885
     886    END SUBROUTINE calc_zenith
     887
     888#if defined ( __rrtmg )
     889!------------------------------------------------------------------------------!
     890! Description:
     891! ------------
     892!-- Calculates surface albedo components based on Briegleb (1992) and
     893!-- Briegleb et al. (1986)
     894!------------------------------------------------------------------------------!
     895    SUBROUTINE calc_albedo
     896
     897        IMPLICIT NONE
     898
     899        IF ( sun_up )  THEN
     900!
     901!--        Ocean
     902           IF ( albedo_type == 1 )  THEN
     903              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
     904                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
     905                                            * ( zenith(0) - 0.5_wp )           &
     906                                            * ( zenith(0) - 1.0_wp )
     907              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
     908!
     909!--        Snow
     910           ELSEIF ( albedo_type == 16 )  THEN
     911              IF ( zenith(0) < 0.5_wp )  THEN
     912                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
     913                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
     914                                     * zenith(0))) - 1.0_wp
     915                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
     916                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
     917                                     * zenith(0))) - 1.0_wp
     918
     919                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
     920                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
     921              ELSE
     922                 rrtm_aldir(0,:,:) = aldif
     923                 rrtm_asdir(0,:,:) = asdif
     924              ENDIF
     925!
     926!--        Sea ice
     927           ELSEIF ( albedo_type == 15 )  THEN
     928                 rrtm_aldir(0,:,:) = aldif
     929                 rrtm_asdir(0,:,:) = asdif
     930!
     931!--        Land surfaces
     932           ELSE
     933              SELECT CASE ( albedo_type )
     934
     935!
     936!--              Surface types with strong zenith dependence
     937                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
     938                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
     939                                        (1.0_wp + 0.8_wp * zenith(0))
     940                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
     941                                        (1.0_wp + 0.8_wp * zenith(0))
     942!
     943!--              Surface types with weak zenith dependence
     944                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
     945                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
     946                                        (1.0_wp + 0.2_wp * zenith(0))
     947                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
     948                                        (1.0_wp + 0.2_wp * zenith(0))
     949
     950                 CASE DEFAULT
     951
     952              END SELECT
     953           ENDIF
     954!
     955!--        Diffusive albedo is taken from Table 2
     956           rrtm_aldif(0,:,:) = aldif
     957           rrtm_asdif(0,:,:) = asdif
     958
     959        ELSE
     960
     961           rrtm_aldir(0,:,:) = 0.0_wp
     962           rrtm_asdir(0,:,:) = 0.0_wp
     963           rrtm_aldif(0,:,:) = 0.0_wp
     964           rrtm_asdif(0,:,:) = 0.0_wp
     965        ENDIF
     966    END SUBROUTINE calc_albedo
     967
     968!------------------------------------------------------------------------------!
     969! Description:
     970! ------------
     971!-- Read sounding data (pressure and temperature) from RADIATION_DATA.
     972!------------------------------------------------------------------------------!
     973    SUBROUTINE read_sounding_data
     974
     975       USE netcdf_control
     976
     977       IMPLICIT NONE
     978
     979       INTEGER(iwp) :: id, id_dim_zrad, id_var, k, nz_snd, nz_snd_start, nz_snd_end
     980
     981       REAL(wp) :: t_surface
     982
     983       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, t_snd_tmp
     984
     985!
     986!--    In case of updates, deallocate arrays first (sufficient to check one
     987!--    array as the others are automatically allocated). This is required
     988!--    because nzt_rad might change during the update
     989       IF ( ALLOCATED ( hyp_snd ) )  THEN
     990          DEALLOCATE( hyp_snd )
     991          DEALLOCATE( t_snd )
     992          DEALLOCATE( q_snd  )
     993          DEALLOCATE ( rrtm_play )
     994          DEALLOCATE ( rrtm_plev )
     995          DEALLOCATE ( rrtm_tlay )
     996          DEALLOCATE ( rrtm_tlev )
     997          DEALLOCATE ( rrtm_h2ovmr )
     998          DEALLOCATE ( rrtm_cicewp )
     999          DEALLOCATE ( rrtm_cldfr )
     1000          DEALLOCATE ( rrtm_cliqwp )
     1001          DEALLOCATE ( rrtm_reice )
     1002          DEALLOCATE ( rrtm_reliq )
     1003          DEALLOCATE ( rrtm_lw_taucld )
     1004          DEALLOCATE ( rrtm_lw_tauaer )
     1005          DEALLOCATE ( rrtm_lwdflx  )
     1006          DEALLOCATE ( rrtm_lwuflx  )
     1007          DEALLOCATE ( rrtm_lwhr  )
     1008          DEALLOCATE ( rrtm_lwuflxc )
     1009          DEALLOCATE ( rrtm_lwdflxc )
     1010          DEALLOCATE ( rrtm_lwhrc )
     1011          DEALLOCATE ( rrtm_sw_taucld )
     1012          DEALLOCATE ( rrtm_sw_ssacld )
     1013          DEALLOCATE ( rrtm_sw_asmcld )
     1014          DEALLOCATE ( rrtm_sw_fsfcld )
     1015          DEALLOCATE ( rrtm_sw_tauaer )
     1016          DEALLOCATE ( rrtm_sw_ssaaer )
     1017          DEALLOCATE ( rrtm_sw_asmaer )
     1018          DEALLOCATE ( rrtm_sw_ecaer )   
     1019          DEALLOCATE ( rrtm_swdflx  )
     1020          DEALLOCATE ( rrtm_swuflx  )
     1021          DEALLOCATE ( rrtm_swhr  )
     1022          DEALLOCATE ( rrtm_swuflxc )
     1023          DEALLOCATE ( rrtm_swdflxc )
     1024          DEALLOCATE ( rrtm_swhrc )
     1025       ENDIF
     1026
     1027!
     1028!--    Open file for reading
     1029       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
     1030       CALL handle_netcdf_error( 'netcdf', 549 )
     1031
     1032!
     1033!--    Inquire dimension of z axis and save in nz_snd
     1034       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
     1035       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
     1036       CALL handle_netcdf_error( 'netcdf', 551 )
     1037
     1038!
     1039! !--    Allocate temporary array for storing pressure data
     1040       ALLOCATE( hyp_snd_tmp(nzb+1:nz_snd) )
     1041       hyp_snd_tmp = 0.0_wp
     1042
     1043
     1044!--    Read pressure from file
     1045       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
     1046       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),    &
     1047                               count = (/nz_snd/) )
     1048       CALL handle_netcdf_error( 'netcdf', 552 )
     1049
     1050!
     1051!--    Allocate temporary array for storing temperature data
     1052       ALLOCATE( t_snd_tmp(nzb+1:nz_snd) )
     1053       t_snd_tmp = 0.0_wp
     1054
     1055!
     1056!--    Read temperature from file
     1057       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
     1058       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),      &
     1059                               count = (/nz_snd/) )
     1060       CALL handle_netcdf_error( 'netcdf', 553 )
     1061
     1062!
     1063!--    Calculate start of sounding data
     1064       nz_snd_start = nz_snd + 1
     1065       nz_snd_end   = nz_snd_end
     1066
     1067!
     1068!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
     1069!--    in Pa, hyp_snd in hPa).
     1070       DO  k = 1, nz_snd
     1071          IF ( hyp_snd_tmp(k) .LT. (hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
     1072             nz_snd_start = k
     1073             EXIT
     1074          END IF
     1075       END DO
     1076
     1077       IF ( nz_snd_start .LE. nz_snd )  THEN
     1078          nz_snd_end = nz_snd - 1
     1079       END IF
     1080
     1081
     1082!
     1083!--    Calculate of total grid points for RRTMG calculations
     1084       nzt_rad = nzt + nz_snd_end - nz_snd_start + 2
     1085
     1086!
     1087!--    Save data above LES domain in hyp_snd, t_snd and q_snd
     1088!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
     1089       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
     1090       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
     1091       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
     1092       hyp_snd = 0.0_wp
     1093       t_snd = 0.0_wp
     1094       q_snd = 0.0_wp
     1095
     1096       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start:nz_snd_end)
     1097       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start:nz_snd_end)
     1098
     1099       DEALLOCATE ( hyp_snd_tmp )
     1100       DEALLOCATE ( t_snd_tmp )
     1101
     1102       nc_stat = NF90_CLOSE( id )
     1103
     1104!
     1105!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
     1106!--    top of the LES domain. This routine does not consider horizontal or
     1107!--    vertical variability of pressure and temperature
     1108       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
     1109       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
     1110
     1111       t_surface = pt_surface * (surface_pressure / 1000.0_wp )**0.286_wp
     1112       DO k = nzb+1, nzt+1
     1113          rrtm_play(0,k) = hyp(k) * 0.01_wp
     1114          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
     1115                         t_surface )**(1.0_wp/0.286_wp)
     1116       ENDDO
     1117
     1118       DO k = nzt+2, nzt_rad
     1119          rrtm_play(0,k) = hyp_snd(k)
     1120          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
     1121       ENDDO
     1122       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
     1123                                   1.5 * hyp_snd(nzt_rad)                      &
     1124                                 - 0.5 * hyp_snd(nzt_rad-1) )
     1125       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
     1126                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
     1127
     1128       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
     1129
     1130!
     1131!--    Calculate temperature/humidity levels at top of the LES domain.
     1132!--    Currently, the temperature is taken from sounding data (might lead to a
     1133!--    temperature jump at interface. To do: Humidity is currently not
     1134!--    calculated above the LES domain.
     1135       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
     1136       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
     1137       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
     1138
     1139       DO k = nzt+8, nzt_rad
     1140          rrtm_tlay(0,k)   = t_snd(k)
     1141          rrtm_h2ovmr(0,k) = q_snd(k)
     1142       ENDDO
     1143       rrtm_tlay(0,nzt_rad+1)   = 2.0_wp * rrtm_tlay(0,nzt_rad)                &
     1144                                  - rrtm_tlay(0,nzt_rad-1)
     1145       DO k = nzt+9, nzt_rad+1
     1146          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
     1147                             - rrtm_tlay(0,k-1))                               &
     1148                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
     1149                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
     1150       ENDDO
     1151       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
     1152
     1153       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
     1154                                  - rrtm_tlev(0,nzt_rad)
     1155!
     1156!--    Allocate remaining RRTMG arrays
     1157       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
     1158       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
     1159       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
     1160       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
     1161       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
     1162       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
     1163       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
     1164       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
     1165       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
     1166       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
     1167       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
     1168       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
     1169       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
     1170       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
     1171       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
     1172
     1173!
     1174!--    The ice phase is currently not considered in PALM
     1175       rrtm_cicewp = 0.0_wp
     1176       rrtm_reice  = 0.0_wp
     1177
     1178!
     1179!--    Set other parameters (move to NAMELIST parameters in the future)
     1180       rrtm_lw_tauaer = 0.0_wp
     1181       rrtm_lw_taucld = 0.0_wp
     1182       rrtm_sw_taucld = 0.0_wp
     1183       rrtm_sw_ssacld = 0.0_wp
     1184       rrtm_sw_asmcld = 0.0_wp
     1185       rrtm_sw_fsfcld = 0.0_wp
     1186       rrtm_sw_tauaer = 0.0_wp
     1187       rrtm_sw_ssaaer = 0.0_wp
     1188       rrtm_sw_asmaer = 0.0_wp
     1189       rrtm_sw_ecaer  = 0.0_wp
     1190
     1191
     1192       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
     1193       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
     1194       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
     1195       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
     1196       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
     1197       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
     1198
     1199       rrtm_swdflx  = 0.0_wp
     1200       rrtm_swuflx  = 0.0_wp
     1201       rrtm_swhr    = 0.0_wp 
     1202       rrtm_swuflxc = 0.0_wp
     1203       rrtm_swdflxc = 0.0_wp
     1204       rrtm_swhrc   = 0.0_wp
     1205
     1206       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
     1207       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
     1208       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
     1209       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
     1210       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
     1211       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
     1212
     1213       rrtm_lwdflx  = 0.0_wp
     1214       rrtm_lwuflx  = 0.0_wp
     1215       rrtm_lwhr    = 0.0_wp 
     1216       rrtm_lwuflxc = 0.0_wp
     1217       rrtm_lwdflxc = 0.0_wp
     1218       rrtm_lwhrc   = 0.0_wp
     1219
     1220
     1221    END SUBROUTINE read_sounding_data
     1222
     1223
     1224!------------------------------------------------------------------------------!
     1225! Description:
     1226! ------------
     1227!-- Read trace gas data from file
     1228!------------------------------------------------------------------------------!
     1229    SUBROUTINE read_trace_gas_data
     1230
     1231       USE netcdf_control
     1232       USE rrsw_ncpar
     1233
     1234       IMPLICIT NONE
     1235
     1236       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !: number of trace gases
     1237
     1238       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              &
     1239           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
     1240                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
     1241
     1242       INTEGER(iwp) :: id, k, m, n, nabs, np, id_abs, id_dim, id_var
     1243
     1244       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
     1245
     1246
     1247       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !: pressure levels for the absorbers
     1248                                                 rrtm_play_tmp, & !: temporary array for pressure zu-levels
     1249                                                 rrtm_plev_tmp, & !: temporary array for pressure zw-levels
     1250                                                 trace_path_tmp   !: temporary array for storing trace gas path data
     1251
     1252       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !: array for storing the absorber amounts
     1253                                                 trace_mls_path, & !: array for storing trace gas path data
     1254                                                 trace_mls_tmp     !: temporary array for storing trace gas data
     1255
     1256
     1257!
     1258!--    In case of updates, deallocate arrays first (sufficient to check one
     1259!--    array as the others are automatically allocated)
     1260       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
     1261          DEALLOCATE ( rrtm_o3vmr  )
     1262          DEALLOCATE ( rrtm_co2vmr )
     1263          DEALLOCATE ( rrtm_ch4vmr )
     1264          DEALLOCATE ( rrtm_n2ovmr )
     1265          DEALLOCATE ( rrtm_o2vmr  )
     1266          DEALLOCATE ( rrtm_cfc11vmr )
     1267          DEALLOCATE ( rrtm_cfc12vmr )
     1268          DEALLOCATE ( rrtm_cfc22vmr )
     1269          DEALLOCATE ( rrtm_ccl4vmr  )
     1270       ENDIF
     1271
     1272!
     1273!--    Allocate trace gas profiles
     1274       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
     1275       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
     1276       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
     1277       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
     1278       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
     1279       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
     1280       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
     1281       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
     1282       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
     1283
     1284!
     1285!--    Open file for reading
     1286       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
     1287       CALL handle_netcdf_error( 'netcdf', 549 )
     1288!
     1289!--    Inquire dimension ids and dimensions
     1290       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
     1291       CALL handle_netcdf_error( 'netcdf', 550 )
     1292       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
     1293       CALL handle_netcdf_error( 'netcdf', 550 )
     1294
     1295       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
     1296       CALL handle_netcdf_error( 'netcdf', 550 )
     1297       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
     1298       CALL handle_netcdf_error( 'netcdf', 550 )
     1299   
     1300
     1301!
     1302!--    Allocate pressure, and trace gas arrays     
     1303       ALLOCATE( p_mls(1:np) )
     1304       ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
     1305       ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
     1306
     1307
     1308       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
     1309       CALL handle_netcdf_error( 'netcdf', 550 )
     1310       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
     1311       CALL handle_netcdf_error( 'netcdf', 550 )
     1312
     1313       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
     1314       CALL handle_netcdf_error( 'netcdf', 550 )
     1315       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
     1316       CALL handle_netcdf_error( 'netcdf', 550 )
     1317
     1318
     1319!
     1320!--    Write absorber amounts (mls) to trace_mls
     1321       DO n = 1, num_trace_gases
     1322          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
     1323
     1324          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
     1325
     1326!
     1327!--       Replace missing values by zero
     1328          WHERE ( trace_mls(n,:) > 2.0_wp ) 
     1329             trace_mls(n,:) = 0.0_wp
     1330          END WHERE
     1331       END DO
     1332
     1333       DEALLOCATE ( trace_mls_tmp )
     1334
     1335       nc_stat = NF90_CLOSE( id )
     1336       CALL handle_netcdf_error( 'netcdf', 551 )
     1337
     1338!
     1339!--    Add extra pressure level for calculations of the trace gas paths
     1340       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
     1341       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
     1342
     1343       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad)
     1344       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
     1345       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
     1346       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
     1347                                         * rrtm_plev(0,nzt_rad+1) )
     1348 
     1349!
     1350!--    Calculate trace gas path (zero at surface) with interpolation to the
     1351!--    sounding levels
     1352       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
     1353
     1354       trace_mls_path(nzb+1,:) = 0.0_wp
     1355       
     1356       DO k = nzb+2, nzt_rad+2
     1357          DO m = 1, num_trace_gases
     1358             trace_mls_path(k,m) = trace_mls_path(k-1,m)
     1359
     1360!
     1361!--          When the pressure level is higher than the trace gas pressure
     1362!--          level, assume that
     1363             IF ( rrtm_plev_tmp(k-1) .GT. p_mls(1) )  THEN             
     1364               
     1365                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
     1366                                      * ( rrtm_plev_tmp(k-1)                   &
     1367                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
     1368                                        ) / g
     1369             ENDIF
     1370
     1371!
     1372!--          Integrate for each sounding level from the contributing p_mls
     1373!--          levels
     1374             DO n = 2, np
     1375!
     1376!--             Limit p_mls so that it is within the model level
     1377                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
     1378                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
     1379                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
     1380                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
     1381
     1382                IF ( p_mls_l .GT. p_mls_u )  THEN
     1383
     1384!
     1385!--                Calculate weights for interpolation
     1386                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
     1387                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
     1388                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
     1389
     1390!
     1391!--                Add level to trace gas path
     1392                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
     1393                                         +  ( p_wgt_u * trace_mls(m,n)         &
     1394                                            + p_wgt_l * trace_mls(m,n-1) )     &
     1395                                         * (p_mls_l + p_mls_u) / g
     1396                ENDIF
     1397             ENDDO
     1398
     1399             IF ( rrtm_plev_tmp(k) .LT. p_mls(np) )  THEN
     1400                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
     1401                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
     1402                                          - rrtm_plev_tmp(k)                   &
     1403                                        ) / g
     1404             ENDIF 
    2431405          ENDDO
    2441406       ENDDO
    2451407
    246        RETURN
    247 
    248     END SUBROUTINE radiation_clearsky
     1408
     1409!
     1410!--    Prepare trace gas path profiles
     1411       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
     1412
     1413       DO m = 1, num_trace_gases
     1414
     1415          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
     1416                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
     1417                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
     1418                                       - rrtm_plev_tmp(2:nzt_rad+2) )
     1419
     1420!
     1421!--       Save trace gas paths to the respective arrays
     1422          SELECT CASE ( TRIM( trace_names(m) ) )
     1423
     1424             CASE ( 'O3' )
     1425
     1426                rrtm_o3vmr(0,:) = trace_path_tmp(:)
     1427
     1428             CASE ( 'CO2' )
     1429
     1430                rrtm_co2vmr(0,:) = trace_path_tmp(:)
     1431
     1432             CASE ( 'CH4' )
     1433
     1434                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
     1435
     1436             CASE ( 'N2O' )
     1437
     1438                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
     1439
     1440             CASE ( 'O2' )
     1441
     1442                rrtm_o2vmr(0,:) = trace_path_tmp(:)
     1443
     1444             CASE ( 'CFC11' )
     1445
     1446                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
     1447
     1448             CASE ( 'CFC12' )
     1449
     1450                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
     1451
     1452             CASE ( 'CFC22' )
     1453
     1454                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
     1455
     1456             CASE ( 'CCL4' )
     1457
     1458                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
     1459
     1460             CASE DEFAULT
     1461
     1462          END SELECT
     1463
     1464       ENDDO
     1465
     1466       DEALLOCATE ( trace_path_tmp )
     1467       DEALLOCATE ( trace_mls_path )
     1468       DEALLOCATE ( rrtm_play_tmp )
     1469       DEALLOCATE ( rrtm_plev_tmp )
     1470       DEALLOCATE ( trace_mls )
     1471       DEALLOCATE ( p_mls )
     1472
     1473    END SUBROUTINE read_trace_gas_data
     1474
     1475#endif
     1476
    2491477
    2501478!------------------------------------------------------------------------------!
    2511479! Description:
    2521480! ------------
    253 !-- Prescribed net radiation
    254 !------------------------------------------------------------------------------!
    255     SUBROUTINE radiation_constant
    256 
    257        rad_net = net_radiation
    258 
    259     END SUBROUTINE radiation_constant
     1481!-- Calculate temperature tendency due to radiative cooling/heating.
     1482!-- Cache-optimized version.
     1483!------------------------------------------------------------------------------!
     1484    SUBROUTINE radiation_tendency_ij ( i, j, tend )
     1485
     1486       USE arrays_3d,                                                          &
     1487           ONLY:  dzw
     1488
     1489       USE cloud_parameters,                                                   &
     1490           ONLY:  pt_d_t, cp
     1491
     1492       USE control_parameters,                                                 &
     1493           ONLY:  rho_surface
     1494
     1495       IMPLICIT NONE
     1496
     1497       INTEGER(iwp) :: i, j, k
     1498
     1499       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
     1500
     1501#if defined ( __rrtmg )
     1502
     1503       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
     1504
     1505       rad_sw_net_l = 0.0_wp
     1506       rad_sw_net_u = 0.0_wp
     1507       rad_lw_net_l = 0.0_wp
     1508       rad_lw_net_u = 0.0_wp
     1509
     1510!
     1511!--    Calculate radiative flux divergence
     1512       DO k = nzb+1, nzt+1
     1513
     1514          rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
     1515          rad_sw_net_u = rad_sw_in(k,j,i)   - rad_sw_out(k,j,i)
     1516          rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
     1517          rad_lw_net_u = rad_lw_in(k,j,i)   - rad_lw_out(k,j,i)
     1518
     1519          tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l            &
     1520                                      + rad_lw_net_u - rad_lw_net_l ) /        &
     1521                                     ( dzw(k) * cp * rho_surface ) * pt_d_t(k)
     1522       ENDDO
     1523
     1524#endif
     1525
     1526    END SUBROUTINE radiation_tendency_ij
     1527
    2601528
    2611529!------------------------------------------------------------------------------!
    2621530! Description:
    2631531! ------------
    264 !-- Implementation of the RRTM radiation_scheme
    265 !------------------------------------------------------------------------------!
    266     SUBROUTINE radiation_rrtm
    267 
    268 
    269     END SUBROUTINE radiation_rrtm
     1532!-- Calculate temperature tendency due to radiative cooling/heating.
     1533!-- Vector-optimized version
     1534!------------------------------------------------------------------------------!
     1535    SUBROUTINE radiation_tendency ( tend )
     1536
     1537       USE arrays_3d,                                                          &
     1538           ONLY:  dzw
     1539
     1540       USE cloud_parameters,                                                   &
     1541           ONLY:  pt_d_t, cp
     1542
     1543       USE indices,                                                            &
     1544           ONLY:  nxl, nxr, nyn, nys
     1545
     1546       USE control_parameters,                                                 &
     1547           ONLY:  rho_surface
     1548
     1549       IMPLICIT NONE
     1550
     1551       INTEGER(iwp) :: i, j, k
     1552
     1553       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
     1554
     1555#if defined ( __rrtmg )
     1556
     1557       REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
     1558
     1559       rad_sw_net_l = 0.0_wp
     1560       rad_sw_net_u = 0.0_wp
     1561       rad_lw_net_l = 0.0_wp
     1562       rad_lw_net_u = 0.0_wp
     1563
     1564       DO  i = nxl, nxr
     1565          DO  j = nys, nyn
     1566
     1567!
     1568!--          Calculate radiative flux divergence
     1569             DO k = nzb+1, nzt+1
     1570
     1571                rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
     1572                rad_sw_net_u = rad_sw_in(k  ,j,i) - rad_sw_out(k  ,j,i)
     1573                rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
     1574                rad_lw_net_u = rad_lw_in(k  ,j,i) - rad_lw_out(k  ,j,i)
     1575
     1576                tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l      &
     1577                                      + rad_lw_net_u - rad_lw_net_l ) /        &
     1578                                      ( dzw(k) * cp * rho_surface ) * pt_d_t(k)
     1579             ENDDO
     1580         ENDDO
     1581       ENDDO
     1582#endif
     1583
     1584    END SUBROUTINE radiation_tendency
    2701585
    2711586 END MODULE radiation_model_mod

  • palm/trunk/SOURCE/read_3d_binary.f90

    r1552 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    111111
    112112    USE radiation_model_mod,                                                   &
    113         ONLY: rad_net_av, rad_sw_in_av
     113        ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
     114              rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av
    114115
    115116    USE random_function_mod,                                                   &
    … …  
    333334!--    First compare the version numbers
    334335       READ ( 13 )  version_on_file
    335        binary_version = '4.0'
     336       binary_version = '4.1'
    336337       IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
    337338          WRITE( message_string, * ) 'version mismatch concerning data ',      &
    … …  
    812813                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    813814
     815                CASE ( 'rad_net' )
     816                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
     817                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
     818                   ENDIF 
     819                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     820                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     821                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     822
    814823                CASE ( 'rad_net_av' )
    815824                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
    … …  
    819828                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    820829                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     830                CASE ( 'rad_lw_in' )
     831                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
     832                      ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     833                   ENDIF 
     834                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     835                   rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     836                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     837
     838                CASE ( 'rad_lw_in_av' )
     839                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
     840                      ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     841                   ENDIF 
     842                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     843                   rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     844                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     845
     846                CASE ( 'rad_lw_out' )
     847                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
     848                      ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     849                   ENDIF 
     850                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     851                   rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     852                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     853
     854                CASE ( 'rad_lw_out_av' )
     855                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
     856                      ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     857                   ENDIF 
     858                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     859                   rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     860                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     861
     862                CASE ( 'rad_sw_in' )
     863                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
     864                      ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     865                   ENDIF 
     866                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     867                   rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     868                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    821869
    822870                CASE ( 'rad_sw_in_av' )
    823871                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    824                       ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
    825                    ENDIF 
    826                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    827                    rad_sw_in_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    828                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     872                      ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     873                   ENDIF 
     874                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     875                   rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     876                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     877
     878                CASE ( 'rad_sw_out' )
     879                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
     880                      ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     881                   ENDIF 
     882                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     883                   rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     884                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     885
     886                CASE ( 'rad_sw_out_av' )
     887                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
     888                      ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     889                   ENDIF 
     890                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     891                   rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     892                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    829893
    830894                CASE ( 'random_iv' )  ! still unresolved issue

  • palm/trunk/SOURCE/read_var_list.f90

    r1561 r1585  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    162162
    163163    USE pegrid
     164
     165    USE radiation_model_mod,                                                   &
     166        ONLY: time_radiation
    164167
    165168    USE statistics,                                                            &
    … …  
    609612          CASE ( 'time_dvrp' )
    610613             READ ( 13 )  time_dvrp
     614          CASE ( 'time_radiation' )
     615             READ ( 13 )  time_radiation
    611616          CASE ( 'time_restart' )
    612617             READ ( 13 )  time_restart

  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r1556 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    112112
    113113    USE radiation_model_mod,                                                   &
    114         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av
     114        ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
     115               rad_sw_out_av, rad_lw_in, rad_lw_in_av, rad_lw_out,             &
     116               rad_lw_out_av
    115117
    116118    IMPLICIT NONE
    … …  
    326328                rad_net_av = 0.0_wp
    327329
    328              CASE ( 'rad_sw_in*' )
     330             CASE ( 'rad_lw_in' )
     331                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
     332                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     333                ENDIF
     334                rad_lw_in_av = 0.0_wp
     335
     336             CASE ( 'rad_lw_out' )
     337                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
     338                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     339                ENDIF
     340                rad_lw_out_av = 0.0_wp
     341
     342             CASE ( 'rad_sw_in' )
    329343                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    330                    ALLOCATE( rad_sw_in_av(nysg:nyng,nxlg:nxrg) )
     344                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    331345                ENDIF
    332346                rad_sw_in_av = 0.0_wp
     347
     348             CASE ( 'rad_sw_out' )
     349                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
     350                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     351                ENDIF
     352                rad_sw_out_av = 0.0_wp
    333353
    334354             CASE ( 'rho' )
    … …  
    732752             ENDDO
    733753
    734           CASE ( 'rad_sw_in*' )
    735              DO  i = nxlg, nxrg
    736                 DO  j = nysg, nyng
    737                    rad_sw_in_av(j,i) = rad_sw_in_av(j,i) + rad_sw_in(j,i)
     754          CASE ( 'rad_lw_in' )
     755             DO  i = nxlg, nxrg
     756                DO  j = nysg, nyng
     757                   DO  k = nzb, nzt+1
     758                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
     759                   ENDDO
     760                ENDDO
     761             ENDDO
     762
     763          CASE ( 'rad_lw_out' )
     764             DO  i = nxlg, nxrg
     765                DO  j = nysg, nyng
     766                   DO  k = nzb, nzt+1
     767                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
     768                   ENDDO
     769                ENDDO
     770             ENDDO
     771
     772
     773          CASE ( 'rad_sw_in' )
     774             DO  i = nxlg, nxrg
     775                DO  j = nysg, nyng
     776                   DO  k = nzb, nzt+1
     777                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
     778                   ENDDO
     779                ENDDO
     780             ENDDO
     781
     782          CASE ( 'rad_sw_out' )
     783             DO  i = nxlg, nxrg
     784                DO  j = nysg, nyng
     785                   DO  k = nzb, nzt+1
     786                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
     787                   ENDDO
    738788                ENDDO
    739789             ENDDO

  • palm/trunk/SOURCE/time_integration.f90

    r1552 r1585  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Moved call of radiation scheme. Added support for RRTM
    2323!
    2424! Former revisions:
    … …  
    226226
    227227    USE radiation_model_mod,                                                   &
    228         ONLY: dt_radiation, irad_scheme, radiation, radiation_clearsky,        &
    229               radiation_constant, radiation_rrtm, time_radiation
     228        ONLY: dt_radiation, radiation, radiation_clearsky,                     &
     229              radiation_rrtmg, radiation_scheme, time_radiation
    230230
    231231    USE statistics,                                                            &
    … …  
    362362          ENDIF
    363363
     364          IF ( radiation .AND. intermediate_timestep_count                     &
     365               == intermediate_timestep_count_max  )  THEN
     366
     367               time_radiation = time_radiation + dt_3d
     368
     369             IF ( time_radiation >= dt_radiation )  THEN
     370
     371                CALL cpu_log( log_point(50), 'radiation', 'start' )
     372
     373                time_radiation = time_radiation - dt_radiation
     374                IF ( radiation_scheme == 'clear-sky' )  THEN
     375                   CALL radiation_clearsky
     376                ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
     377                   CALL radiation_rrtmg
     378                ENDIF
     379
     380                CALL cpu_log( log_point(50), 'radiation', 'stop' )
     381             ENDIF
     382          ENDIF
    364383!
    365384!--       Solve the prognostic equations. A fast cache optimized version with
    … …  
    597616!--       2) run soil model
    598617          IF ( land_surface )  THEN
    599 
    600              time_radiation = time_radiation + dt_3d
    601 
    602              IF ( time_radiation >= dt_radiation )  THEN
    603 
    604                 CALL cpu_log( log_point(50), 'radiation', 'start' )
    605 
    606                 time_radiation = time_radiation - dt_radiation
    607                 IF ( irad_scheme == 0 )  THEN
    608                    CALL radiation_constant             
    609                 ELSEIF ( irad_scheme == 1 )  THEN
    610                    CALL radiation_clearsky
    611                 ELSEIF ( irad_scheme == 2 )  THEN
    612                    CALL radiation_rrtm
    613                 ENDIF
    614 
    615                 CALL cpu_log( log_point(50), 'radiation', 'stop' )
    616              ENDIF
    617618
    618619             CALL cpu_log( log_point(54), 'land_surface', 'start' )
    … …  
    693694             !$acc update device( vpt )
    694695          ENDIF
     696
    695697!
    696698!--       Compute the diffusion quantities

  • palm/trunk/SOURCE/user_init_land_surface.f90

    r1497 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Changed description text
    2323!
    2424! Former revisions:
    … …  
    5252
    5353!
    54 !-- Here the user-defined grid initializing actions follow:
     54!-- Here the user-defined land surface initialization actions follow:
    5555
    5656

  • palm/trunk/SOURCE/write_3d_binary.f90

    r1552 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    102102   
    103103    USE radiation_model_mod,                                                   &
    104         ONLY: radiation, rad_net_av, rad_sw_in_av
     104        ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
     105              rad_lw_out, rad_lw_out_av, rad_sw_in, rad_sw_in_av, rad_sw_out,  &
     106              rad_sw_out_av
    105107
    106108    USE random_function_mod,                                                   &
    … …  
    124126!
    125127!-- Write arrays.
    126     binary_version = '4.0'
     128    binary_version = '4.1'
    127129
    128130    WRITE ( 14 )  binary_version
    … …  
    262264       ENDIF
    263265    ENDIF
     266    IF ( ALLOCATED( rad_net ) )  THEN
     267       WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
     268    ENDIF
    264269    IF ( radiation )  THEN
    265270       IF ( ALLOCATED( rad_net_av ) )  THEN
    266271          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
    267272       ENDIF 
     273       IF ( ALLOCATED( rad_lw_in ) )  THEN
     274          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
     275       ENDIF
     276       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
     277          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
     278       ENDIF
     279       IF ( ALLOCATED( rad_lw_out ) )  THEN
     280          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out
     281       ENDIF
     282       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
     283          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
     284       ENDIF
     285       IF ( ALLOCATED( rad_sw_in ) )  THEN
     286          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
     287       ENDIF
    268288       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
    269           WRITE ( 14 )  'rad_sw_in_av          ';  WRITE ( 14 )  rad_sw_in_av 
     289          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
     290       ENDIF
     291       IF ( ALLOCATED( rad_sw_out ) )  THEN
     292          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
     293       ENDIF
     294       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
     295          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
    270296       ENDIF
    271297    ENDIF

  • palm/trunk/SOURCE/write_var_list.f90

    r1552 r1585  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adapted for RRTMG
    2323!
    2424! Former revisions:
    … …  
    147147   
    148148    USE pegrid
     149
     150    USE radiation_model_mod,                                                   &
     151        ONLY: time_radiation
    149152
    150153    USE statistics,                                                            &
    … …  
    522525    WRITE ( 14 )  'time_dvrp                     '
    523526    WRITE ( 14 )  time_dvrp
     527    WRITE ( 14 )  'time_radiation                '
     528    WRITE ( 14 )  time_radiation
    524529    WRITE ( 14 )  'time_restart                  '
    525530    WRITE ( 14 )  time_restart
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |