Home


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

Changeset 1976 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jul 27, 2016 1:28:04 PM (8 years ago)
Author:
maronga
Message:

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

Location:
palm/trunk/SOURCE
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/Makefile

    r1973 r1976  
    509509mod_particle_attributes.o: mod_particle_attributes.f90 mod_kinds.o
    510510netcdf_interface_mod.o: netcdf_interface_mod.f90 modules.o mod_kinds.o \
    511    land_surface_model_mod.o spectra_mod.o
     511   land_surface_model_mod.o radiation_model_mod.o spectra_mod.o
    512512nudging_mod.o: modules.o cpulog_mod.o mod_kinds.o
    513513package_parin.o: modules.o mod_kinds.o mod_particle_attributes.o

  • palm/trunk/SOURCE/average_3d_data.f90

    r1973 r1976  
    101101
    102102    USE radiation_model_mod,                                                   &
    103         ONLY:  rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,       &
    104                rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
    105                rad_lw_hr_av, rad_sw_in, rad_sw_in_av, rad_sw_out,              &
    106                rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
    107                rad_sw_hr_av
     103        ONLY:  radiation, radiation_3d_data_averaging
    108104
    109105
    … …  
    295291             ENDDO
    296292
    297          CASE ( 'rad_net*' )
    298              DO  i = nxlg, nxrg
    299                 DO  j = nysg, nyng
    300                    rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
    301                 ENDDO
    302              ENDDO
    303 
    304           CASE ( 'rad_lw_in' )
    305              DO  i = nxlg, nxrg
    306                 DO  j = nysg, nyng
    307                    DO  k = nzb, nzt+1
    308                       rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    309                    ENDDO
    310                 ENDDO
    311              ENDDO
    312 
    313           CASE ( 'rad_lw_out' )
    314              DO  i = nxlg, nxrg
    315                 DO  j = nysg, nyng
    316                    DO  k = nzb, nzt+1
    317                       rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    318                    ENDDO
    319                 ENDDO
    320              ENDDO
    321 
    322           CASE ( 'rad_lw_cs_hr' )
    323              DO  i = nxlg, nxrg
    324                 DO  j = nysg, nyng
    325                    DO  k = nzb, nzt+1
    326                       rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    327                    ENDDO
    328                 ENDDO
    329              ENDDO
    330 
    331           CASE ( 'rad_lw_hr' )
    332              DO  i = nxlg, nxrg
    333                 DO  j = nysg, nyng
    334                    DO  k = nzb, nzt+1
    335                       rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    336                    ENDDO
    337                 ENDDO
    338              ENDDO
    339 
    340           CASE ( 'rad_sw_in' )
    341              DO  i = nxlg, nxrg
    342                 DO  j = nysg, nyng
    343                    DO  k = nzb, nzt+1
    344                       rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    345                    ENDDO
    346                 ENDDO
    347              ENDDO
    348 
    349           CASE ( 'rad_sw_out' )
    350              DO  i = nxlg, nxrg
    351                 DO  j = nysg, nyng
    352                    DO  k = nzb, nzt+1
    353                       rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    354                    ENDDO
    355                 ENDDO
    356              ENDDO
    357 
    358           CASE ( 'rad_sw_cs_hr' )
    359              DO  i = nxlg, nxrg
    360                 DO  j = nysg, nyng
    361                    DO  k = nzb, nzt+1
    362                       rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    363                    ENDDO
    364                 ENDDO
    365              ENDDO
    366 
    367           CASE ( 'rad_sw_hr' )
    368              DO  i = nxlg, nxrg
    369                 DO  j = nysg, nyng
    370                    DO  k = nzb, nzt+1
    371                       rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
    372                    ENDDO
    373                 ENDDO
    374              ENDDO
    375 
    376293          CASE ( 'rho' )
    377294             DO  i = nxlg, nxrg
    … …  
    487404
    488405!
     406!--          Radiation quantity
     407             IF ( radiation )  THEN
     408                CALL radiation_3d_data_averaging( 'average', doav(ii) )
     409             ENDIF
     410
     411!
    489412!--          User-defined quantity
    490413             CALL user_3d_data_averaging( 'average', doav(ii) )

  • palm/trunk/SOURCE/data_output_2d.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Output of radiation quantities is now done directly in the respective module
    2222!
    2323! Former revisions:
    … …  
    2626!
    2727! 1972 2016-07-26 07:52:02Z maronga
    28 ! Output of land surface quantities is now done directly in the respective module
     28! Output of land surface quantities is now done directly in the respective
     29! module
    2930!
    3031! 1960 2016-07-12 16:34:24Z suehring
    … …  
    195196
    196197    USE radiation_model_mod,                                                   &
    197         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
    198                rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
    199                rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
    200                rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
    201                rad_lw_hr_av
     198        ONLY:  radiation, radiation_data_output_2d
    202199
    203200    IMPLICIT NONE
    … …  
    753750                IF ( mode == 'xy' )  level_z = zu
    754751
    755              CASE ( 'rad_net*_xy' )        ! 2d-array
    756                 IF ( av == 0 ) THEN
    757                    DO  i = nxlg, nxrg
    758                       DO  j = nysg, nyng
    759                          local_pf(i,j,nzb+1) =  rad_net(j,i)
    760                       ENDDO
    761                    ENDDO
    762                 ELSE
    763                    DO  i = nxlg, nxrg
    764                       DO  j = nysg, nyng
    765                          local_pf(i,j,nzb+1) =  rad_net_av(j,i)
    766                       ENDDO
    767                    ENDDO
    768                 ENDIF
    769                 resorted = .TRUE.
    770                 two_d = .TRUE.
    771                 level_z(nzb+1) = zu(nzb+1)
    772 
    773 
    774              CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
    775                 IF ( av == 0 )  THEN
    776                    to_be_resorted => rad_lw_in
    777                 ELSE
    778                    to_be_resorted => rad_lw_in_av
    779                 ENDIF
    780                 IF ( mode == 'xy' )  level_z = zu
    781 
    782              CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
    783                 IF ( av == 0 )  THEN
    784                    to_be_resorted => rad_lw_out
    785                 ELSE
    786                    to_be_resorted => rad_lw_out_av
    787                 ENDIF
    788                 IF ( mode == 'xy' )  level_z = zu
    789 
    790              CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
    791                 IF ( av == 0 )  THEN
    792                    to_be_resorted => rad_lw_cs_hr
    793                 ELSE
    794                    to_be_resorted => rad_lw_cs_hr_av
    795                 ENDIF
    796                 IF ( mode == 'xy' )  level_z = zw
    797 
    798              CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
    799                 IF ( av == 0 )  THEN
    800                    to_be_resorted => rad_lw_hr
    801                 ELSE
    802                    to_be_resorted => rad_lw_hr_av
    803                 ENDIF
    804                 IF ( mode == 'xy' )  level_z = zw
    805 
    806              CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
    807                 IF ( av == 0 )  THEN
    808                    to_be_resorted => rad_sw_in
    809                 ELSE
    810                    to_be_resorted => rad_sw_in_av
    811                 ENDIF
    812                 IF ( mode == 'xy' )  level_z = zu
    813 
    814              CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
    815                 IF ( av == 0 )  THEN
    816                    to_be_resorted => rad_sw_out
    817                 ELSE
    818                    to_be_resorted => rad_sw_out_av
    819                 ENDIF
    820                 IF ( mode == 'xy' )  level_z = zu
    821 
    822              CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
    823                 IF ( av == 0 )  THEN
    824                    to_be_resorted => rad_sw_cs_hr
    825                 ELSE
    826                    to_be_resorted => rad_sw_cs_hr_av
    827                 ENDIF
    828                 IF ( mode == 'xy' )  level_z = zw
    829 
    830              CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
    831                 IF ( av == 0 )  THEN
    832                    to_be_resorted => rad_sw_hr
    833                 ELSE
    834                    to_be_resorted => rad_sw_hr_av
    835                 ENDIF
    836                 IF ( mode == 'xy' )  level_z = zw
     752
    837753
    838754             CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
    … …  
    1034950                   CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
    1035951                                            local_pf, two_d, nzb_do, nzt_do )
     952                ENDIF
     953
     954!
     955!--             Radiation quantity
     956                IF ( .NOT. found  .AND.  radiation )  THEN
     957                   CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
     958                                                  mode, local_pf, two_d  )
    1036959                ENDIF
    1037960

  • palm/trunk/SOURCE/data_output_3d.f90

    r1973 r1976  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Output of radiation quantities is now done directly in the respective module
    2222!
    2323! Former revisions:
    … …  
    165165
    166166    USE radiation_model_mod,                                                   &
    167         ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
    168                rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,         &
    169                rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,             &
    170                rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
     167        ONLY:  radiation, radiation_data_output_3d
    171168
    172169
    … …  
    494491             ENDIF
    495492
    496           CASE ( 'rad_sw_in' )
    497              IF ( av == 0 )  THEN
    498                 to_be_resorted => rad_sw_in
    499              ELSE
    500                 to_be_resorted => rad_sw_in_av
    501              ENDIF
    502 
    503           CASE ( 'rad_sw_out' )
    504              IF ( av == 0 )  THEN
    505                 to_be_resorted => rad_sw_out
    506              ELSE
    507                 to_be_resorted => rad_sw_out_av
    508              ENDIF
    509 
    510           CASE ( 'rad_sw_cs_hr' )
    511              IF ( av == 0 )  THEN
    512                 to_be_resorted => rad_sw_cs_hr
    513              ELSE
    514                 to_be_resorted => rad_sw_cs_hr_av
    515              ENDIF
    516 
    517           CASE ( 'rad_sw_hr' )
    518              IF ( av == 0 )  THEN
    519                 to_be_resorted => rad_sw_hr
    520              ELSE
    521                 to_be_resorted => rad_sw_hr_av
    522              ENDIF
    523 
    524           CASE ( 'rad_lw_in' )
    525              IF ( av == 0 )  THEN
    526                 to_be_resorted => rad_lw_in
    527              ELSE
    528                 to_be_resorted => rad_lw_in_av
    529              ENDIF
    530 
    531           CASE ( 'rad_lw_out' )
    532              IF ( av == 0 )  THEN
    533                 to_be_resorted => rad_lw_out
    534              ELSE
    535                 to_be_resorted => rad_lw_out_av
    536              ENDIF
    537 
    538           CASE ( 'rad_lw_cs_hr' )
    539              IF ( av == 0 )  THEN
    540                 to_be_resorted => rad_lw_cs_hr
    541              ELSE
    542                 to_be_resorted => rad_lw_cs_hr_av
    543              ENDIF
    544 
    545           CASE ( 'rad_lw_hr' )
    546              IF ( av == 0 )  THEN
    547                 to_be_resorted => rad_lw_hr
    548              ELSE
    549                 to_be_resorted => rad_lw_hr_av
    550              ENDIF
    551 
    552493          CASE ( 'rho' )
    553494             IF ( av == 0 )  THEN
    … …  
    613554
    614555                CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
     556                resorted = .TRUE.
     557
     558!
     559!--             If no soil model variable was found, re-allocate local_pf
     560                IF ( .NOT. found )  THEN
     561                   nzb_do = nzb
     562                   nzt_do = nz_do3d
     563
     564                   DEALLOCATE ( local_pf )
     565                   ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb_do:nzt_do) )                 
     566                ENDIF
     567
     568             ENDIF
     569
     570!
     571!--          Radiation quantity
     572             IF ( .NOT. found  .AND.  radiation )  THEN
     573                CALL radiation_data_output_3d( av, do3d(av,if), found,         &
     574                                               local_pf )
    615575                resorted = .TRUE.
    616576             ENDIF

  • palm/trunk/SOURCE/data_output_mask.f90

    r1961 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Output of radiation quantities is now done directly in the respective module
    2222!
    2323! Former revisions:
    … …  
    127127
    128128    USE radiation_model_mod,                                                   &
    129         ONLY:  rad_lw_in, rad_lw_in_av, rad_lw_out, rad_lw_out_av,             &
    130                rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,         &
    131                rad_sw_in, rad_sw_in_av, rad_sw_out, rad_sw_out_av,             &
    132                rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
    133 
     129        ONLY:  radiation, radiation_data_output_mask
    134130
    135131    IMPLICIT NONE
    … …  
    425421             ENDIF
    426422
    427           CASE ( 'rad_lw_in' )
    428              IF ( av == 0 )  THEN
    429                 to_be_resorted => rad_lw_in
    430              ELSE
    431                 to_be_resorted => rad_lw_in_av
    432              ENDIF
    433 
    434           CASE ( 'rad_lw_out' )
    435              IF ( av == 0 )  THEN
    436                 to_be_resorted => rad_lw_out
    437              ELSE
    438                 to_be_resorted => rad_lw_out_av
    439              ENDIF
    440 
    441           CASE ( 'rad_lw_cs_hr' )
    442              IF ( av == 0 )  THEN
    443                 to_be_resorted => rad_lw_cs_hr
    444              ELSE
    445                 to_be_resorted => rad_lw_cs_hr_av
    446              ENDIF
    447 
    448           CASE ( 'rad_lw_hr' )
    449              IF ( av == 0 )  THEN
    450                 to_be_resorted => rad_lw_hr
    451              ELSE
    452                 to_be_resorted => rad_lw_hr_av
    453              ENDIF
    454 
    455           CASE ( 'rad_sw_in' )
    456              IF ( av == 0 )  THEN
    457                 to_be_resorted => rad_sw_in
    458              ELSE
    459                 to_be_resorted => rad_sw_in_av
    460              ENDIF
    461 
    462           CASE ( 'rad_sw_out' )
    463              IF ( av == 0 )  THEN
    464                 to_be_resorted => rad_sw_out
    465              ELSE
    466                 to_be_resorted => rad_sw_out_av
    467              ENDIF
    468 
    469           CASE ( 'rad_sw_cs_hr' )
    470              IF ( av == 0 )  THEN
    471                 to_be_resorted => rad_sw_cs_hr
    472              ELSE
    473                 to_be_resorted => rad_sw_cs_hr_av
    474              ENDIF
    475 
    476           CASE ( 'rad_sw_hr' )
    477              IF ( av == 0 )  THEN
    478                 to_be_resorted => rad_sw_hr
    479              ELSE
    480                 to_be_resorted => rad_sw_hr_av
    481              ENDIF
    482 
    483423          CASE ( 'rho' )
    484424             IF ( av == 0 )  THEN
    … …  
    531471
    532472          CASE DEFAULT
     473
     474!
     475!--          Radiation quantity
     476             IF ( radiation )  THEN
     477                CALL radiation_data_output_mask(av, domask(mid,av,if), found,  &
     478                                                local_pf )
     479             ENDIF
     480
    533481!
    534482!--          User defined quantity
    535              CALL user_data_output_mask(av, domask(mid,av,if), found, local_pf )
     483             IF ( .NOT. found )  THEN
     484                CALL user_data_output_mask(av, domask(mid,av,if), found,       &
     485                                           local_pf )
     486             ENDIF
     487
    536488             resorted = .TRUE.
    537489

  • palm/trunk/SOURCE/flow_statistics.f90

    r1961 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Removed some unneeded __rrtmg preprocessor directives
    2222!
    2323! Former revisions:
    … …  
    832832                ENDIF
    833833#endif
    834 
    835834             ENDIF
    836 
    837835!
    838836!--          Subgridscale fluxes at the top surface
    … …  
    14471445
    14481446          IF ( radiation_scheme == 'rrtmg' )  THEN
    1449 #if defined ( __rrtmg )
    14501447             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
    14511448             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
    … …  
    14571454             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
    14581455             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
    1459 #endif
    14601456          ENDIF
    1461 
    1462 
    14631457       ENDIF
    14641458
    … …  
    16431637          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
    16441638
    1645 #if defined ( __rrtmg )
    16461639          IF ( radiation_scheme == 'rrtmg' )  THEN
    16471640             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
    … …  
    16501643             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
    16511644          ENDIF
    1652 #endif
    16531645
    16541646       ENDIF
    … …  
    36193611          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
    36203612
    3621 #if defined ( __rrtmg )
    36223613          IF ( radiation_scheme == 'rrtmg' )  THEN
    36233614             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
    … …  
    36263617             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
    36273618          ENDIF
    3628 #endif
    36293619
    36303620       ENDIF

  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Parts of the code have been reformatted. Use of radiation model output is
     22! generalized and simplified. Added more output quantities due to modularization
    2223!
    2324! Former revisions:
    … …  
    178179
    179180    USE radiation_model_mod,                                                   &
    180         ONLY:  force_radiation_call, radiation_scheme, rad_net, rad_sw_in,     &
    181                rad_lw_out, rad_lw_out_change_0, sigma_sb,                      &
    182                unscheduled_radiation_calls
     181        ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
     182               rad_lw_out_change_0, unscheduled_radiation_calls
    183183       
    184 #if defined ( __rrtmg )
    185     USE radiation_model_mod,                                                   &
    186         ONLY:  rrtm_idrv
    187 #endif
    188 
    189184    USE statistics,                                                            &
    190185        ONLY:  hom, statistic_regions
    … …  
    636631    END INTERFACE lsm_read_restart_data
    637632
    638 
    639633    INTERFACE lsm_last_actions
    640634       MODULE PROCEDURE lsm_last_actions
    … …  
    648642!> Check data output for land surface model
    649643!------------------------------------------------------------------------------!
    650     SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
     644 SUBROUTINE lsm_check_data_output( var, unit, i, ilen, k )
    651645 
    652646 
    653        USE control_parameters,                                                 &
    654            ONLY: data_output, message_string
    655 
    656        IMPLICIT NONE
    657 
    658        CHARACTER (LEN=*) ::  unit     !<
    659        CHARACTER (LEN=*) ::  var !<
    660 
    661        INTEGER(iwp) :: i
    662        INTEGER(iwp) :: ilen   
    663        INTEGER(iwp) :: k
    664 
    665        SELECT CASE ( TRIM( var ) )
    666 
    667           CASE ( 'm_soil' )
    668              IF (  .NOT.  land_surface )  THEN
    669                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    670                          'res land_surface = .TRUE.'
    671                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    672              ENDIF
    673              unit = 'm3/m3'
     647    USE control_parameters,                                                 &
     648        ONLY:  data_output, message_string
     649
     650    IMPLICIT NONE
     651
     652    CHARACTER (LEN=*) ::  unit     !<
     653    CHARACTER (LEN=*) ::  var !<
     654
     655    INTEGER(iwp) :: i
     656    INTEGER(iwp) :: ilen   
     657    INTEGER(iwp) :: k
     658
     659    SELECT CASE ( TRIM( var ) )
     660
     661       CASE ( 'm_soil' )
     662          IF (  .NOT.  land_surface )  THEN
     663             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     664                      'res land_surface = .TRUE.'
     665             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     666          ENDIF
     667          unit = 'm3/m3'
     668           
     669       CASE ( 't_soil' )
     670          IF (  .NOT.  land_surface )  THEN
     671             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     672                      'res land_surface = .TRUE.'
     673             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     674          ENDIF
     675          unit = 'K'   
    674676             
    675           CASE ( 't_soil' )
    676              IF (  .NOT.  land_surface )  THEN
    677                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    678                          'res land_surface = .TRUE.'
    679                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    680              ENDIF
    681              unit = 'K'   
     677       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
     678              'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
     679              'r_a*', 'r_s*', 'shf_eb*' )
     680          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
     681             message_string = 'illegal value for data_output: "' //         &
     682                              TRIM( var ) // '" & only 2d-horizontal ' //   &
     683                              'cross sections are allowed for this value'
     684             CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
     685          ENDIF
     686          IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
     687             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     688                              'res land_surface = .TRUE.'
     689             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     690          ENDIF
     691          IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
     692             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     693                              'res land_surface = .TRUE.'
     694             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     695          ENDIF
     696          IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
     697             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     698                              'res land_surface = .TRUE.'
     699             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     700          ENDIF
     701          IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
     702             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     703                              'res land_surface = .TRUE.'
     704             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     705          ENDIF
     706          IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
     707             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     708                              'res land_surface = .TRUE.'
     709             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     710          ENDIF
     711          IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
     712             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     713                              'res land_surface = .TRUE.'
     714             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     715          ENDIF
     716          IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
     717             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     718                              'res land_surface = .TRUE.'
     719             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     720          ENDIF
     721          IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
     722          THEN
     723             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     724                              'res land_surface = .TRUE.'
     725             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     726          ENDIF
     727          IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
     728          THEN
     729             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     730                              'res land_surface = .TRUE.'
     731             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     732          ENDIF
     733          IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
     734          THEN
     735             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     736                              'res land_surface = .TRUE.'
     737             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     738          ENDIF
     739          IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
     740          THEN
     741             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     742                              'res land_surface = .TRUE.'
     743             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     744          ENDIF
     745          IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
     746          THEN
     747             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     748                              'res land_surface = .TRUE.'
     749             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
     750          ENDIF
     751
     752          IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
     753          IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
     754          IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
     755          IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
     756          IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
     757          IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
     758          IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
     759          IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
     760          IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
     761          IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
     762          IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
     763          IF ( TRIM( var ) == 'r_s*')     unit = 's/m'
     764          IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
    682765             
    683           CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
    684                  'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
    685                  'r_a*', 'r_s*', 'shf_eb*' )
    686              IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    687                 message_string = 'illegal value for data_output: "' //         &
    688                                  TRIM( var ) // '" & only 2d-horizontal ' //   &
    689                                  'cross sections are allowed for this value'
    690                 CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
    691              ENDIF
    692              IF ( TRIM( var ) == 'lai*'  .AND.  .NOT.  land_surface )  THEN
    693                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    694                                  'res land_surface = .TRUE.'
    695                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    696              ENDIF
    697              IF ( TRIM( var ) == 'c_liq*'  .AND.  .NOT.  land_surface )  THEN
    698                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    699                                  'res land_surface = .TRUE.'
    700                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    701              ENDIF
    702              IF ( TRIM( var ) == 'c_soil*'  .AND.  .NOT.  land_surface )  THEN
    703                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    704                                  'res land_surface = .TRUE.'
    705                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    706              ENDIF
    707              IF ( TRIM( var ) == 'c_veg*'  .AND.  .NOT. land_surface )  THEN
    708                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    709                                  'res land_surface = .TRUE.'
    710                 CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    711              ENDIF
    712              IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
    713                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    714                                  'res land_surface = .TRUE.'
    715                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    716              ENDIF
    717              IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
    718                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    719                                  'res land_surface = .TRUE.'
    720                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    721              ENDIF
    722              IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
    723                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    724                                  'res land_surface = .TRUE.'
    725                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    726              ENDIF
    727              IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
    728              THEN
    729                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    730                                  'res land_surface = .TRUE.'
    731                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    732              ENDIF
    733              IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
    734              THEN
    735                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    736                                  'res land_surface = .TRUE.'
    737                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    738              ENDIF
    739              IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
    740              THEN
    741                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    742                                  'res land_surface = .TRUE.'
    743                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    744              ENDIF
    745              IF ( TRIM( var ) == 'r_a*'  .AND.  .NOT.  land_surface ) &
    746              THEN
    747                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    748                                  'res land_surface = .TRUE.'
    749                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    750              ENDIF
    751              IF ( TRIM( var ) == 'r_s*'  .AND.  .NOT.  land_surface ) &
    752              THEN
    753                 message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    754                                  'res land_surface = .TRUE.'
    755                 CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    756              ENDIF
    757 
    758              IF ( TRIM( var ) == 'lai*'   )  unit = 'none'
    759              IF ( TRIM( var ) == 'c_liq*' )  unit = 'none'
    760              IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
    761              IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
    762              IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
    763              IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
    764              IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
    765              IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
    766              IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
    767              IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
    768              IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
    769              IF ( TRIM( var ) == 'r_s*')     unit = 's/m'
    770              IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
    771              
    772           CASE DEFAULT
    773              unit = 'illegal'
    774 
    775        END SELECT
    776 
    777 
    778     END SUBROUTINE lsm_check_data_output
    779 
     766       CASE DEFAULT
     767          unit = 'illegal'
     768
     769    END SELECT
     770
     771
     772 END SUBROUTINE lsm_check_data_output
     773
     774
     775!------------------------------------------------------------------------------!
     776! Description:
     777! ------------
     778!> Check data output of profiles for land surface model
     779!------------------------------------------------------------------------------!
    780780 SUBROUTINE lsm_check_data_output_pr( variable, var_count, unit, dopr_unit )
    781781 
    782        USE control_parameters,                                                 &
    783            ONLY: data_output_pr, message_string
    784 
    785        USE indices
    786 
    787        USE profil_parameter
    788 
    789        USE statistics
    790 
    791        IMPLICIT NONE
     782    USE control_parameters,                                                 &
     783        ONLY: data_output_pr, message_string
     784
     785    USE indices
     786
     787    USE profil_parameter
     788
     789    USE statistics
     790
     791    IMPLICIT NONE
    792792   
    793        CHARACTER (LEN=*) ::  unit      !<
    794        CHARACTER (LEN=*) ::  variable  !<
    795        CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
     793    CHARACTER (LEN=*) ::  unit      !<
     794    CHARACTER (LEN=*) ::  variable  !<
     795    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
    796796 
    797        INTEGER(iwp) ::  user_pr_index !<
    798        INTEGER(iwp) ::  var_count     !<
    799 
    800        SELECT CASE ( TRIM( variable ) )
     797    INTEGER(iwp) ::  user_pr_index !<
     798    INTEGER(iwp) ::  var_count     !<
     799
     800    SELECT CASE ( TRIM( variable ) )
    801801       
    802           CASE ( 't_soil', '#t_soil' )
    803              IF (  .NOT.  land_surface )  THEN
    804                 message_string = 'data_output_pr = ' //                        &
    805                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    806                                  'not implemented for land_surface = .FALSE.'
    807                 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
    808              ELSE
    809                 dopr_index(var_count) = 89
    810                 dopr_unit     = 'K'
    811                 hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
    812                 IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
    813                    dopr_initial_index(var_count) = 90
    814                    hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
    815                    data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    816                 ENDIF
    817                 unit = dopr_unit
     802       CASE ( 't_soil', '#t_soil' )
     803          IF (  .NOT.  land_surface )  THEN
     804             message_string = 'data_output_pr = ' //                        &
     805                              TRIM( data_output_pr(var_count) ) // ' is' // &
     806                              'not implemented for land_surface = .FALSE.'
     807             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
     808          ELSE
     809             dopr_index(var_count) = 89
     810             dopr_unit     = 'K'
     811             hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     812             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
     813                dopr_initial_index(var_count) = 90
     814                hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     815                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    818816             ENDIF
    819 
    820           CASE ( 'm_soil', '#m_soil' )
    821              IF (  .NOT.  land_surface )  THEN
    822                 message_string = 'data_output_pr = ' //                        &
    823                                  TRIM( data_output_pr(var_count) ) // ' is' // &
    824                                  ' not implemented for land_surface = .FALSE.'
    825                 CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
    826              ELSE
    827                 dopr_index(var_count) = 91
    828                 dopr_unit     = 'm3/m3'
    829                 hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
    830                 IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
    831                    dopr_initial_index(var_count) = 92
    832                    hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
    833                    data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    834                 ENDIF
    835                  unit = dopr_unit
     817             unit = dopr_unit
     818          ENDIF
     819
     820       CASE ( 'm_soil', '#m_soil' )
     821          IF (  .NOT.  land_surface )  THEN
     822             message_string = 'data_output_pr = ' //                        &
     823                              TRIM( data_output_pr(var_count) ) // ' is' // &
     824                              ' not implemented for land_surface = .FALSE.'
     825             CALL message( 'check_parameters', 'PA0402', 1, 2, 0, 6, 0 )
     826          ELSE
     827             dopr_index(var_count) = 91
     828             dopr_unit     = 'm3/m3'
     829             hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     830             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
     831                dopr_initial_index(var_count) = 92
     832                hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     833                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    836834             ENDIF
    837 
    838 
    839           CASE DEFAULT
    840              unit = 'illegal'
    841 
    842        END SELECT
    843 
    844 
    845     END SUBROUTINE lsm_check_data_output_pr
     835             unit = dopr_unit
     836          ENDIF
     837
     838
     839       CASE DEFAULT
     840          unit = 'illegal'
     841
     842    END SELECT
     843
     844
     845 END SUBROUTINE lsm_check_data_output_pr
    846846 
    847847 
    … …  
    851851!> Check parameters routine for land surface model
    852852!------------------------------------------------------------------------------!
    853     SUBROUTINE lsm_check_parameters
    854 
    855        USE control_parameters,                                                 &
    856            ONLY: bc_pt_b, bc_q_b, constant_flux_layer, message_string,         &
    857                  most_method, topography
     853 SUBROUTINE lsm_check_parameters
     854
     855    USE control_parameters,                                                    &
     856        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
     857               most_method, topography
    858858                 
    859        USE radiation_model_mod,                                                &
    860            ONLY: radiation
     859    USE radiation_model_mod,                                                   &
     860        ONLY: radiation
    861861   
    862862   
    863        IMPLICIT NONE
     863    IMPLICIT NONE
    864864
    865865 
    866866!
    867 !--    Dirichlet boundary conditions are required as the surface fluxes are
    868 !--    calculated from the temperature/humidity gradients in the land surface
    869 !--    model
    870        IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
    871           message_string = 'lsm requires setting of'//                         &
    872                            'bc_pt_b = "dirichlet" and '//                      &
    873                            'bc_q_b  = "dirichlet"'
    874           CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
     867!-- Dirichlet boundary conditions are required as the surface fluxes are
     868!-- calculated from the temperature/humidity gradients in the land surface
     869!-- model
     870    IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
     871       message_string = 'lsm requires setting of'//                         &
     872                        'bc_pt_b = "dirichlet" and '//                      &
     873                        'bc_q_b  = "dirichlet"'
     874       CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
     875    ENDIF
     876
     877    IF (  .NOT.  constant_flux_layer )  THEN
     878       message_string = 'lsm requires '//                                   &
     879                        'constant_flux_layer = .T.'
     880       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     881    ENDIF
     882
     883    IF ( topography /= 'flat' )  THEN
     884       message_string = 'lsm cannot be used ' //                            &
     885                        'in combination with  topography /= "flat"'
     886       CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
     887    ENDIF
     888
     889    IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
     890           most_method == 'lookup' )  THEN
     891        WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
     892                                   'allowed in combination with ',          &
     893                                   'most_method = ', most_method
     894       CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     895    ENDIF
     896
     897    IF ( veg_type == 0 )  THEN
     898       IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
     899          message_string = 'veg_type = 0 (user_defined)'//                  &
     900                           'requires setting of root_fraction(0:3)'//       &
     901                           '/= 9999999.9 and SUM(root_fraction) = 1'
     902          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    875903       ENDIF
    876 
    877        IF (  .NOT.  constant_flux_layer )  THEN
    878           message_string = 'lsm requires '//                                   &
    879                            'constant_flux_layer = .T.'
    880           CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     904 
     905       IF ( min_canopy_resistance == 9999999.9_wp )  THEN
     906          message_string = 'veg_type = 0 (user defined)'//                  &
     907                           'requires setting of min_canopy_resistance'//    &
     908                           '/= 9999999.9'
     909          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    881910       ENDIF
    882911
    883        IF ( topography /= 'flat' )  THEN
    884           message_string = 'lsm cannot be used ' //                            &
    885                            'in combination with  topography /= "flat"'
    886           CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
     912       IF ( leaf_area_index == 9999999.9_wp )  THEN
     913          message_string = 'veg_type = 0 (user_defined)'//                  &
     914                           'requires setting of leaf_area_index'//          &
     915                           '/= 9999999.9'
     916          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    887917       ENDIF
    888918
    889        IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
    890               most_method == 'lookup' )  THEN
    891            WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
    892                                       'allowed in combination with ',          &
    893                                       'most_method = ', most_method
    894           CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     919       IF ( vegetation_coverage == 9999999.9_wp )  THEN
     920          message_string = 'veg_type = 0 (user_defined)'//                  &
     921                           'requires setting of vegetation_coverage'//      &
     922                           '/= 9999999.9'
     923             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    895924       ENDIF
    896925
    897        IF ( veg_type == 0 )  THEN
    898           IF ( SUM( root_fraction ) /= 1.0_wp )  THEN
    899              message_string = 'veg_type = 0 (user_defined)'//                  &
    900                               'requires setting of root_fraction(0:3)'//       &
    901                               '/= 9999999.9 and SUM(root_fraction) = 1'
    902              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    903           ENDIF
    904  
    905           IF ( min_canopy_resistance == 9999999.9_wp )  THEN
    906              message_string = 'veg_type = 0 (user defined)'//                  &
    907                               'requires setting of min_canopy_resistance'//    &
    908                               '/= 9999999.9'
    909              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    910           ENDIF
    911 
    912           IF ( leaf_area_index == 9999999.9_wp )  THEN
    913              message_string = 'veg_type = 0 (user_defined)'//                  &
    914                               'requires setting of leaf_area_index'//          &
    915                               '/= 9999999.9'
    916              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    917           ENDIF
    918 
    919           IF ( vegetation_coverage == 9999999.9_wp )  THEN
    920              message_string = 'veg_type = 0 (user_defined)'//                  &
    921                               'requires setting of vegetation_coverage'//      &
    922                               '/= 9999999.9'
    923              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    924           ENDIF
    925 
    926           IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
    927              message_string = 'veg_type = 0 (user_defined)'//                  &
    928                               'requires setting of'//                          &
    929                               'canopy_resistance_coefficient /= 9999999.9'
    930              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    931           ENDIF
    932 
    933           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    934              message_string = 'veg_type = 0 (user_defined)'//                  &
    935                               'requires setting of lambda_surface_stable'//    &
    936                               '/= 9999999.9'
    937              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    938           ENDIF
    939 
    940           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    941              message_string = 'veg_type = 0 (user_defined)'//                  &
    942                               'requires setting of lambda_surface_unstable'//  &
    943                               '/= 9999999.9'
    944              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    945           ENDIF
    946 
    947           IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
    948              message_string = 'veg_type = 0 (user_defined)'//                  &
    949                               'requires setting of f_shortwave_incoming'//     &
    950                               '/= 9999999.9'
    951              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    952           ENDIF
    953 
    954           IF ( z0_eb == 9999999.9_wp )  THEN
    955              message_string = 'veg_type = 0 (user_defined)'//                  &
    956                               'requires setting of z0_eb'//                    &
    957                               '/= 9999999.9'
    958              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    959           ENDIF
    960 
    961           IF ( z0h_eb == 9999999.9_wp )  THEN
    962              message_string = 'veg_type = 0 (user_defined)'//                  &
    963                               'requires setting of z0h_eb'//                   &
    964                               '/= 9999999.9'
    965              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    966           ENDIF
    967 
    968 
     926       IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
     927          message_string = 'veg_type = 0 (user_defined)'//                  &
     928                           'requires setting of'//                          &
     929                           'canopy_resistance_coefficient /= 9999999.9'
     930          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    969931       ENDIF
    970932
    971        IF ( soil_type == 0 )  THEN
    972 
    973           IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
    974              message_string = 'soil_type = 0 (user_defined)'//                 &
    975                               'requires setting of alpha_vangenuchten'//       &
    976                               '/= 9999999.9'
    977              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    978           ENDIF
    979 
    980           IF ( l_vangenuchten == 9999999.9_wp )  THEN
    981              message_string = 'soil_type = 0 (user_defined)'//                 &
    982                               'requires setting of l_vangenuchten'//           &
    983                               '/= 9999999.9'
    984              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    985           ENDIF
    986 
    987           IF ( n_vangenuchten == 9999999.9_wp )  THEN
    988              message_string = 'soil_type = 0 (user_defined)'//                 &
    989                               'requires setting of n_vangenuchten'//           &
    990                               '/= 9999999.9'
    991              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    992           ENDIF
    993 
    994           IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
    995              message_string = 'soil_type = 0 (user_defined)'//                 &
    996                               'requires setting of hydraulic_conductivity'//   &
    997                               '/= 9999999.9'
    998              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    999           ENDIF
    1000 
    1001           IF ( saturation_moisture == 9999999.9_wp )  THEN
    1002              message_string = 'soil_type = 0 (user_defined)'//                 &
    1003                               'requires setting of saturation_moisture'//      &
    1004                               '/= 9999999.9'
    1005              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    1006           ENDIF
    1007 
    1008           IF ( field_capacity == 9999999.9_wp )  THEN
    1009              message_string = 'soil_type = 0 (user_defined)'//                 &
    1010                               'requires setting of field_capacity'//           &
    1011                               '/= 9999999.9'
    1012              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    1013           ENDIF
    1014 
    1015           IF ( wilting_point == 9999999.9_wp )  THEN
    1016              message_string = 'soil_type = 0 (user_defined)'//                 &
    1017                               'requires setting of wilting_point'//            &
    1018                               '/= 9999999.9'
    1019              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    1020           ENDIF
    1021 
    1022           IF ( residual_moisture == 9999999.9_wp )  THEN
    1023              message_string = 'soil_type = 0 (user_defined)'//                 &
    1024                               'requires setting of residual_moisture'//        &
    1025                               '/= 9999999.9'
    1026              CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    1027           ENDIF
    1028 
     933       IF ( lambda_surface_stable == 9999999.9_wp )  THEN
     934          message_string = 'veg_type = 0 (user_defined)'//                  &
     935                           'requires setting of lambda_surface_stable'//    &
     936                           '/= 9999999.9'
     937          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1029938       ENDIF
    1030939
    1031        IF (  .NOT.  radiation )  THEN
    1032           message_string = 'lsm requires '//                                   &
    1033                            'radiation = .T.'
    1034           CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     940       IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
     941          message_string = 'veg_type = 0 (user_defined)'//                  &
     942                           'requires setting of lambda_surface_unstable'//  &
     943                           '/= 9999999.9'
     944          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1035945       ENDIF
     946
     947       IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
     948          message_string = 'veg_type = 0 (user_defined)'//                  &
     949                           'requires setting of f_shortwave_incoming'//     &
     950                           '/= 9999999.9'
     951          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     952       ENDIF
     953
     954       IF ( z0_eb == 9999999.9_wp )  THEN
     955          message_string = 'veg_type = 0 (user_defined)'//                  &
     956                           'requires setting of z0_eb'//                    &
     957                           '/= 9999999.9'
     958          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     959       ENDIF
     960
     961       IF ( z0h_eb == 9999999.9_wp )  THEN
     962          message_string = 'veg_type = 0 (user_defined)'//                  &
     963                           'requires setting of z0h_eb'//                   &
     964                           '/= 9999999.9'
     965          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     966       ENDIF
     967
     968
     969    ENDIF
     970
     971    IF ( soil_type == 0 )  THEN
     972
     973       IF ( alpha_vangenuchten == 9999999.9_wp )  THEN
     974          message_string = 'soil_type = 0 (user_defined)'//                 &
     975                           'requires setting of alpha_vangenuchten'//       &
     976                           '/= 9999999.9'
     977          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     978       ENDIF
     979
     980       IF ( l_vangenuchten == 9999999.9_wp )  THEN
     981          message_string = 'soil_type = 0 (user_defined)'//                 &
     982                           'requires setting of l_vangenuchten'//           &
     983                           '/= 9999999.9'
     984          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     985       ENDIF
     986
     987       IF ( n_vangenuchten == 9999999.9_wp )  THEN
     988          message_string = 'soil_type = 0 (user_defined)'//                 &
     989                           'requires setting of n_vangenuchten'//           &
     990                           '/= 9999999.9'
     991          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     992       ENDIF
     993
     994       IF ( hydraulic_conductivity == 9999999.9_wp )  THEN
     995          message_string = 'soil_type = 0 (user_defined)'//                 &
     996                           'requires setting of hydraulic_conductivity'//   &
     997                           '/= 9999999.9'
     998          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     999       ENDIF
     1000
     1001       IF ( saturation_moisture == 9999999.9_wp )  THEN
     1002          message_string = 'soil_type = 0 (user_defined)'//                 &
     1003                           'requires setting of saturation_moisture'//      &
     1004                           '/= 9999999.9'
     1005          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1006       ENDIF
     1007
     1008       IF ( field_capacity == 9999999.9_wp )  THEN
     1009          message_string = 'soil_type = 0 (user_defined)'//                 &
     1010                           'requires setting of field_capacity'//           &
     1011                           '/= 9999999.9'
     1012          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1013       ENDIF
     1014
     1015       IF ( wilting_point == 9999999.9_wp )  THEN
     1016          message_string = 'soil_type = 0 (user_defined)'//                 &
     1017                           'requires setting of wilting_point'//            &
     1018                           '/= 9999999.9'
     1019          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1020       ENDIF
     1021
     1022       IF ( residual_moisture == 9999999.9_wp )  THEN
     1023          message_string = 'soil_type = 0 (user_defined)'//                 &
     1024                           'requires setting of residual_moisture'//        &
     1025                           '/= 9999999.9'
     1026          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1027       ENDIF
     1028
     1029    ENDIF
     1030
     1031    IF (  .NOT.  radiation )  THEN
     1032       message_string = 'lsm requires '//                                   &
     1033                        'radiation = .T.'
     1034       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
     1035    ENDIF
    10361036       
    10371037       
    1038     END SUBROUTINE lsm_check_parameters
     1038 END SUBROUTINE lsm_check_parameters
    10391039 
    10401040!------------------------------------------------------------------------------!
    … …  
    10431043!> Solver for the energy balance at the surface.
    10441044!------------------------------------------------------------------------------!
    1045     SUBROUTINE lsm_energy_balance
    1046 
    1047 
    1048        IMPLICIT NONE
    1049 
    1050        INTEGER(iwp) ::  i         !< running index
    1051        INTEGER(iwp) ::  j         !< running index
    1052        INTEGER(iwp) ::  k, ks     !< running index
    1053 
    1054        REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
    1055                    f1,          & !< resistance correction term 1
    1056                    f2,          & !< resistance correction term 2
    1057                    f3,          & !< resistance correction term 3
    1058                    m_min,       & !< minimum soil moisture
    1059                    e,           & !< water vapour pressure
    1060                    e_s,         & !< water vapour saturation pressure
    1061                    e_s_dt,      & !< derivate of e_s with respect to T
    1062                    tend,        & !< tendency
    1063                    dq_s_dt,     & !< derivate of q_s with respect to T
    1064                    coef_1,      & !< coef. for prognostic equation
    1065                    coef_2,      & !< coef. for prognostic equation
    1066                    f_qsws,      & !< factor for qsws_eb
    1067                    f_qsws_veg,  & !< factor for qsws_veg_eb
    1068                    f_qsws_soil, & !< factor for qsws_soil_eb
    1069                    f_qsws_liq,  & !< factor for qsws_liq_eb
    1070                    f_shf,       & !< factor for shf_eb
    1071                    lambda_surface, & !< Current value of lambda_surface
    1072                    m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
    1073                    pt1,         & !< potential temperature at first grid level
    1074                    qv1            !< specific humidity at first grid level
    1075 
    1076 !
    1077 !--    Calculate the exner function for the current time step
    1078        exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    1079 
    1080        DO  i = nxlg, nxrg
    1081           DO  j = nysg, nyng
    1082              k = nzb_s_inner(j,i)
    1083 
    1084 !
    1085 !--          Set lambda_surface according to stratification between skin layer and soil
    1086              IF (  .NOT.  pave_surface(j,i) )  THEN
    1087 
    1088                 c_surface_tmp = c_surface
    1089 
    1090                 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
    1091                    lambda_surface = lambda_surface_s(j,i)
    1092                 ELSE
    1093                    lambda_surface = lambda_surface_u(j,i)
    1094                 ENDIF
     1045 SUBROUTINE lsm_energy_balance
     1046
     1047
     1048    IMPLICIT NONE
     1049
     1050    INTEGER(iwp) ::  i         !< running index
     1051    INTEGER(iwp) ::  j         !< running index
     1052    INTEGER(iwp) ::  k, ks     !< running index
     1053
     1054    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
     1055                f1,          & !< resistance correction term 1
     1056                f2,          & !< resistance correction term 2
     1057                f3,          & !< resistance correction term 3
     1058                m_min,       & !< minimum soil moisture
     1059                e,           & !< water vapour pressure
     1060                e_s,         & !< water vapour saturation pressure
     1061                e_s_dt,      & !< derivate of e_s with respect to T
     1062                tend,        & !< tendency
     1063                dq_s_dt,     & !< derivate of q_s with respect to T
     1064                coef_1,      & !< coef. for prognostic equation
     1065                coef_2,      & !< coef. for prognostic equation
     1066                f_qsws,      & !< factor for qsws_eb
     1067                f_qsws_veg,  & !< factor for qsws_veg_eb
     1068                f_qsws_soil, & !< factor for qsws_soil_eb
     1069                f_qsws_liq,  & !< factor for qsws_liq_eb
     1070                f_shf,       & !< factor for shf_eb
     1071                lambda_surface, & !< Current value of lambda_surface
     1072                m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
     1073                pt1,         & !< potential temperature at first grid level
     1074                qv1            !< specific humidity at first grid level
     1075
     1076!
     1077!-- Calculate the exner function for the current time step
     1078    exn = ( surface_pressure / 1000.0_wp )**0.286_wp
     1079
     1080    DO  i = nxlg, nxrg
     1081       DO  j = nysg, nyng
     1082          k = nzb_s_inner(j,i)
     1083
     1084!
     1085!--       Set lambda_surface according to stratification between skin layer and soil
     1086          IF (  .NOT.  pave_surface(j,i) )  THEN
     1087
     1088             c_surface_tmp = c_surface
     1089
     1090             IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
     1091                lambda_surface = lambda_surface_s(j,i)
    10951092             ELSE
    1096 
    1097                 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
    1098                 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
    1099 
     1093                lambda_surface = lambda_surface_u(j,i)
    11001094             ENDIF
    1101 
    1102 !
    1103 !--          First step: calculate aerodyamic resistance. As pt, us, ts
    1104 !--          are not available for the prognostic time step, data from the last
    1105 !--          time step is used here. Note that this formulation is the
    1106 !--          equivalent to the ECMWF formulation using drag coefficients
    1107              IF ( cloud_physics )  THEN
    1108                 pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
    1109                 qv1 = q(k+1,j,i) - ql(k+1,j,i)
    1110              ELSE
    1111                 pt1 = pt(k+1,j,i)
    1112                 qv1 = q(k+1,j,i)
    1113              ENDIF
    1114 
    1115              r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
    1116 
    1117 !
    1118 !--          Make sure that the resistance does not drop to zero for neutral
    1119 !--          stratification
    1120              IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
    1121 
    1122 !
    1123 !--          Second step: calculate canopy resistance r_canopy
    1124 !--          f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
     1095          ELSE
     1096
     1097             c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
     1098             lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
     1099
     1100          ENDIF
     1101
     1102!
     1103!--       First step: calculate aerodyamic resistance. As pt, us, ts
     1104!--       are not available for the prognostic time step, data from the last
     1105!--       time step is used here. Note that this formulation is the
     1106!--       equivalent to the ECMWF formulation using drag coefficients
     1107          IF ( cloud_physics )  THEN
     1108             pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
     1109             qv1 = q(k+1,j,i) - ql(k+1,j,i)
     1110          ELSE
     1111             pt1 = pt(k+1,j,i)
     1112             qv1 = q(k+1,j,i)
     1113          ENDIF
     1114
     1115          r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
     1116
     1117!
     1118!--       Make sure that the resistance does not drop to zero for neutral
     1119!--       stratification
     1120          IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
     1121
     1122!
     1123!--       Second step: calculate canopy resistance r_canopy
     1124!--       f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
    11251125 
    1126 !--          f1: correction for incoming shortwave radiation (stomata close at
    1127 !--          night)
    1128              IF ( radiation_scheme /= 'constant' )  THEN
    1129                 f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /  &
    1130                               (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)          &
    1131                                + 1.0_wp)) )
    1132              ELSE
    1133                 f1 = 1.0_wp
    1134              ENDIF
    1135 
    1136 
    1137 !
    1138 !--          f2: correction for soil moisture availability to plants (the
    1139 !--          integrated soil moisture must thus be considered here)
    1140 !--          f2 = 0 for very dry soils
    1141              m_total = 0.0_wp
    1142              DO  ks = nzb_soil, nzt_soil
    1143                  m_total = m_total + root_fr(ks,j,i)                           &
    1144                            * MAX(m_soil(ks,j,i),m_wilt(j,i))
    1145              ENDDO
    1146 
    1147              IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
    1148                 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
    1149              ELSEIF ( m_total >= m_fc(j,i) )  THEN
    1150                 f2 = 1.0_wp
    1151              ELSE
    1152                 f2 = 1.0E-20_wp
    1153              ENDIF
    1154 
    1155 !
    1156 !--          Calculate water vapour pressure at saturation
    1157              e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)     &
    1158                            - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
    1159 
    1160 !
    1161 !--          f3: correction for vapour pressure deficit
    1162              IF ( g_d(j,i) /= 0.0_wp )  THEN
    1163 !
    1164 !--             Calculate vapour pressure
    1165                 e  = qv1 * surface_pressure / 0.622_wp
    1166                 f3 = EXP ( -g_d(j,i) * (e_s - e) )
    1167              ELSE
    1168                 f3 = 1.0_wp
    1169              ENDIF
    1170 
    1171 !
    1172 !--          Calculate canopy resistance. In case that c_veg is 0 (bare soils),
    1173 !--          this calculation is obsolete, as r_canopy is not used below.
    1174 !--          To do: check for very dry soil -> r_canopy goes to infinity
    1175              r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3      &
    1176                                              + 1.0E-20_wp)
    1177 
    1178 !
    1179 !--          Third step: calculate bare soil resistance r_soil. The Clapp &
    1180 !--          Hornberger parametrization does not consider c_veg.
    1181              IF ( soil_type_2d(j,i) /= 7 )  THEN
    1182                 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *     &
    1183                         m_res(j,i)
    1184              ELSE
    1185                 m_min = m_wilt(j,i)
    1186              ENDIF
    1187 
    1188              f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
    1189              f2 = MAX(f2,1.0E-20_wp)
    1190              f2 = MIN(f2,1.0_wp)
    1191 
    1192              r_soil(j,i) = r_soil_min(j,i) / f2
    1193 
    1194 !
    1195 !--          Calculate the maximum possible liquid water amount on plants and
    1196 !--          bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
    1197 !--          assumed, while paved surfaces might hold up 1 mm of water. The
    1198 !--          liquid water fraction for paved surfaces is calculated after
    1199 !--          Noilhan & Planton (1989), while the ECMWF formulation is used for
    1200 !--          vegetated surfaces and bare soils.
    1201              IF ( pave_surface(j,i) )  THEN
    1202                 m_liq_eb_max = m_max_depth * 5.0_wp
    1203                 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
    1204              ELSE
    1205                 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)           &
    1206                                + (1.0_wp - c_veg(j,i)) )
    1207                 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
    1208              ENDIF
    1209 
    1210 !
    1211 !--          Calculate saturation specific humidity
    1212              q_s = 0.622_wp * e_s / surface_pressure
    1213 
    1214 !
    1215 !--          In case of dewfall, set evapotranspiration to zero
    1216 !--          All super-saturated water is then removed from the air
    1217              IF ( humidity  .AND.  q_s <= qv1 )  THEN
    1218                 r_canopy(j,i) = 0.0_wp
    1219                 r_soil(j,i)   = 0.0_wp
    1220              ENDIF
    1221 
    1222 !
    1223 !--          Calculate coefficients for the total evapotranspiration
    1224 !--          In case of water surface, set vegetation and soil fluxes to zero.
    1225 !--          For pavements, only evaporation of liquid water is possible.
    1226              IF ( water_surface(j,i) )  THEN
    1227                 f_qsws_veg  = 0.0_wp
    1228                 f_qsws_soil = 0.0_wp
    1229                 f_qsws_liq  = rho_lv / r_a(j,i)
    1230              ELSEIF ( pave_surface (j,i) )  THEN
    1231                 f_qsws_veg  = 0.0_wp
    1232                 f_qsws_soil = 0.0_wp
    1233                 f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
    1234              ELSE
    1235                 f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/     &
    1236                               (r_a(j,i) + r_canopy(j,i))
    1237                 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +     &
    1238                                                           r_soil(j,i))
    1239                 f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
    1240              ENDIF
    1241 !
    1242 !--          If soil moisture is below wilting point, plants do no longer
    1243 !--          transpirate.
    1244 !              IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
    1245 !                 f_qsws_veg = 0.0_wp
    1246 !              ENDIF
    1247 
    1248              f_shf  = rho_cp / r_a(j,i)
    1249              f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
    1250 
    1251 !
    1252 !--          Calculate derivative of q_s for Taylor series expansion
    1253              e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -        &
    1254                               17.269_wp*(t_surface(j,i) - 273.16_wp)           &
    1255                               / (t_surface(j,i) - 35.86_wp)**2 )
    1256 
    1257              dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
    1258 
    1259 !
    1260 !--          Add LW up so that it can be removed in prognostic equation
    1261              rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
    1262 
    1263 !
    1264 !--          Calculate new skin temperature
    1265              IF ( humidity )  THEN
    1266 #if defined ( __rrtmg )
    1267 !
    1268 !--             Numerator of the prognostic equation
    1269                 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
    1270                          * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
    1271                          + f_shf * pt1 + f_qsws * ( qv1 - q_s                  &
    1272                          + dq_s_dt * t_surface(j,i) ) + lambda_surface         &
    1273                          * t_soil(nzb_soil,j,i)
    1274 
    1275 !
    1276 !--             Denominator of the prognostic equation
    1277                 coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt           &
    1278                          + lambda_surface + f_shf / exn
    1279 #else
    1280 
    1281 !
    1282 !--             Numerator of the prognostic equation
    1283                 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
    1284                          * t_surface(j,i) ** 4                                 &
    1285                          + f_shf * pt1 + f_qsws * ( qv1                        &
    1286                          - q_s + dq_s_dt * t_surface(j,i) )                    &
    1287                          + lambda_surface * t_soil(nzb_soil,j,i)
    1288 
    1289 !
    1290 !--             Denominator of the prognostic equation
    1291                 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws      &
    1292                          * dq_s_dt + lambda_surface + f_shf / exn
    1293  
    1294 #endif
    1295              ELSE
    1296 
    1297 #if defined ( __rrtmg )
    1298 !
    1299 !--             Numerator of the prognostic equation
    1300                 coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)             &
    1301                          * t_surface(j,i) - rad_lw_out(nzb,j,i)                &
    1302                          + f_shf * pt1  + lambda_surface                       &
    1303                          * t_soil(nzb_soil,j,i)
    1304 
    1305 !
    1306 !--             Denominator of the prognostic equation
    1307                 coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
    1308 #else
    1309 
    1310 !
    1311 !--             Numerator of the prognostic equation
    1312                 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb                    &
    1313                           * t_surface(j,i) ** 4 + f_shf * pt1                  &
    1314                           + lambda_surface * t_soil(nzb_soil,j,i)
    1315 
    1316 !
    1317 !--             Denominator of the prognostic equation
    1318                 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3               &
    1319                          + lambda_surface + f_shf / exn
    1320 #endif
    1321              ENDIF
    1322 
    1323              tend = 0.0_wp
    1324 
    1325 !
    1326 !--          Implicit solution when the surface layer has no heat capacity,
    1327 !--          otherwise use RK3 scheme.
    1328              t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *    &
    1329                                 t_surface(j,i) ) / ( c_surface_tmp + coef_2    &
     1126!--       f1: correction for incoming shortwave radiation (stomata close at
     1127!--       night)
     1128          f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /        &
     1129                           (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)             &
     1130                            + 1.0_wp)) )
     1131
     1132
     1133
     1134!
     1135!--       f2: correction for soil moisture availability to plants (the
     1136!--       integrated soil moisture must thus be considered here)
     1137!--       f2 = 0 for very dry soils
     1138          m_total = 0.0_wp
     1139          DO  ks = nzb_soil, nzt_soil
     1140              m_total = m_total + root_fr(ks,j,i)                              &
     1141                        * MAX(m_soil(ks,j,i),m_wilt(j,i))
     1142          ENDDO
     1143
     1144          IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
     1145             f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
     1146          ELSEIF ( m_total >= m_fc(j,i) )  THEN
     1147             f2 = 1.0_wp
     1148          ELSE
     1149             f2 = 1.0E-20_wp
     1150          ENDIF
     1151
     1152!
     1153!--       Calculate water vapour pressure at saturation
     1154          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)        &
     1155                        - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
     1156
     1157!
     1158!--       f3: correction for vapour pressure deficit
     1159          IF ( g_d(j,i) /= 0.0_wp )  THEN
     1160!
     1161!--          Calculate vapour pressure
     1162             e  = qv1 * surface_pressure / 0.622_wp
     1163             f3 = EXP ( -g_d(j,i) * (e_s - e) )
     1164          ELSE
     1165             f3 = 1.0_wp
     1166          ENDIF
     1167
     1168!
     1169!--       Calculate canopy resistance. In case that c_veg is 0 (bare soils),
     1170!--       this calculation is obsolete, as r_canopy is not used below.
     1171!--       To do: check for very dry soil -> r_canopy goes to infinity
     1172          r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3         &
     1173                                          + 1.0E-20_wp)
     1174
     1175!
     1176!--       Third step: calculate bare soil resistance r_soil. The Clapp &
     1177!--       Hornberger parametrization does not consider c_veg.
     1178          IF ( soil_type_2d(j,i) /= 7 )  THEN
     1179             m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
     1180                     m_res(j,i)
     1181          ELSE
     1182             m_min = m_wilt(j,i)
     1183          ENDIF
     1184
     1185          f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
     1186          f2 = MAX(f2,1.0E-20_wp)
     1187          f2 = MIN(f2,1.0_wp)
     1188
     1189          r_soil(j,i) = r_soil_min(j,i) / f2
     1190
     1191!
     1192!--       Calculate the maximum possible liquid water amount on plants and
     1193!--       bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
     1194!--       assumed, while paved surfaces might hold up 1 mm of water. The
     1195!--       liquid water fraction for paved surfaces is calculated after
     1196!--       Noilhan & Planton (1989), while the ECMWF formulation is used for
     1197!--       vegetated surfaces and bare soils.
     1198          IF ( pave_surface(j,i) )  THEN
     1199             m_liq_eb_max = m_max_depth * 5.0_wp
     1200             c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
     1201          ELSE
     1202             m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)              &
     1203                            + (1.0_wp - c_veg(j,i)) )
     1204             c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
     1205          ENDIF
     1206
     1207!
     1208!--       Calculate saturation specific humidity
     1209          q_s = 0.622_wp * e_s / surface_pressure
     1210
     1211!
     1212!--       In case of dewfall, set evapotranspiration to zero
     1213!--       All super-saturated water is then removed from the air
     1214          IF ( humidity  .AND.  q_s <= qv1 )  THEN
     1215             r_canopy(j,i) = 0.0_wp
     1216             r_soil(j,i)   = 0.0_wp
     1217          ENDIF
     1218
     1219!
     1220!--       Calculate coefficients for the total evapotranspiration
     1221!--       In case of water surface, set vegetation and soil fluxes to zero.
     1222!--       For pavements, only evaporation of liquid water is possible.
     1223          IF ( water_surface(j,i) )  THEN
     1224             f_qsws_veg  = 0.0_wp
     1225             f_qsws_soil = 0.0_wp
     1226             f_qsws_liq  = rho_lv / r_a(j,i)
     1227          ELSEIF ( pave_surface (j,i) )  THEN
     1228             f_qsws_veg  = 0.0_wp
     1229             f_qsws_soil = 0.0_wp
     1230             f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
     1231          ELSE
     1232             f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
     1233                           (r_a(j,i) + r_canopy(j,i))
     1234             f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
     1235                                                             r_soil(j,i))
     1236             f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
     1237          ENDIF
     1238!
     1239!--       If soil moisture is below wilting point, plants do no longer
     1240!--       transpirate.
     1241!           IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
     1242!              f_qsws_veg = 0.0_wp
     1243!           ENDIF
     1244
     1245          f_shf  = rho_cp / r_a(j,i)
     1246          f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
     1247
     1248!
     1249!--       Calculate derivative of q_s for Taylor series expansion
     1250          e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -           &
     1251                           17.269_wp*(t_surface(j,i) - 273.16_wp)              &
     1252                           / (t_surface(j,i) - 35.86_wp)**2 )
     1253
     1254          dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
     1255
     1256!
     1257!--       Add LW up so that it can be removed in prognostic equation
     1258          rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
     1259
     1260!
     1261!--       Calculate new skin temperature
     1262          IF ( humidity )  THEN
     1263
     1264!
     1265!--          Numerator of the prognostic equation
     1266             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
     1267                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
     1268                      + f_shf * pt1 + f_qsws * ( qv1 - q_s                     &
     1269                      + dq_s_dt * t_surface(j,i) ) + lambda_surface            &
     1270                      * t_soil(nzb_soil,j,i)
     1271
     1272!
     1273!--          Denominator of the prognostic equation
     1274             coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt              &
     1275                      + lambda_surface + f_shf / exn
     1276          ELSE
     1277
     1278!
     1279!--          Numerator of the prognostic equation
     1280             coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
     1281                      * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
     1282                      + f_shf * pt1  + lambda_surface                          &
     1283                      * t_soil(nzb_soil,j,i)
     1284
     1285!
     1286!--          Denominator of the prognostic equation
     1287             coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
     1288
     1289          ENDIF
     1290
     1291          tend = 0.0_wp
     1292
     1293!
     1294!--       Implicit solution when the surface layer has no heat capacity,
     1295!--       otherwise use RK3 scheme.
     1296          t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *       &
     1297                             t_surface(j,i) ) / ( c_surface_tmp + coef_2       &
    13301298                                * dt_3d * tsc(2) )
    13311299
    13321300!
    1333 !--          Add RK3 term
    1334              IF ( c_surface_tmp /= 0.0_wp )  THEN
    1335 
    1336                 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)           &
    1337                                    * tt_surface_m(j,i)
    1338 
    1339 !
    1340 !--             Calculate true tendency
    1341                 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)     &
    1342                        * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
    1343 !
    1344 !--             Calculate t_surface tendencies for the next Runge-Kutta step
    1345                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1346                    IF ( intermediate_timestep_count == 1 )  THEN
    1347                       tt_surface_m(j,i) = tend
    1348                    ELSEIF ( intermediate_timestep_count <                      &
    1349                             intermediate_timestep_count_max )  THEN
    1350                       tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp        &
    1351                                           * tt_surface_m(j,i)
    1352                    ENDIF
     1301!--       Add RK3 term
     1302          IF ( c_surface_tmp /= 0.0_wp )  THEN
     1303
     1304             t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
     1305                                * tt_surface_m(j,i)
     1306
     1307!
     1308!--          Calculate true tendency
     1309             tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
     1310                    * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
     1311!
     1312!--          Calculate t_surface tendencies for the next Runge-Kutta step
     1313             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1314                IF ( intermediate_timestep_count == 1 )  THEN
     1315                   tt_surface_m(j,i) = tend
     1316                ELSEIF ( intermediate_timestep_count <                         &
     1317                         intermediate_timestep_count_max )  THEN
     1318                   tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
     1319                                       * tt_surface_m(j,i)
    13531320                ENDIF
    13541321             ENDIF
    1355 
    1356 !
    1357 !--          In case of fast changes in the skin temperature, it is possible to
    1358 !--          update the radiative fluxes independently from the prescribed
    1359 !--          radiation call frequency. This effectively prevents oscillations,
    1360 !--          especially when setting skip_time_do_radiation /= 0. The threshold
    1361 !--          value of 0.2 used here is just a first guess. This method should be
    1362 !--          revised in the future as tests have shown that the threshold is
    1363 !--          often reached, when no oscillations would occur (causes immense
    1364 !--          computing time for the radiation code).
    1365              IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.     &
    1366                   unscheduled_radiation_calls )  THEN
    1367                 force_radiation_call_l = .TRUE.
     1322          ENDIF
     1323
     1324!
     1325!--       In case of fast changes in the skin temperature, it is possible to
     1326!--       update the radiative fluxes independently from the prescribed
     1327!--       radiation call frequency. This effectively prevents oscillations,
     1328!--       especially when setting skip_time_do_radiation /= 0. The threshold
     1329!--       value of 0.2 used here is just a first guess. This method should be
     1330!--       revised in the future as tests have shown that the threshold is
     1331!--       often reached, when no oscillations would occur (causes immense
     1332!--       computing time for the radiation code).
     1333          IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.        &
     1334               unscheduled_radiation_calls )  THEN
     1335             force_radiation_call_l = .TRUE.
     1336          ENDIF
     1337
     1338          pt(k,j,i) = t_surface_p(j,i) / exn
     1339
     1340!
     1341!--       Calculate fluxes
     1342          rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)         &
     1343                             * t_surface(j,i) - rad_lw_out(nzb,j,i)            &
     1344                             - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
     1345
     1346          rad_net(j,i) = rad_net_l(j,i)
     1347          rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) &
     1348                                * ( t_surface_p(j,i) - t_surface(j,i) )
     1349
     1350          ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)                  &
     1351                           - t_soil(nzb_soil,j,i))
     1352
     1353          shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
     1354
     1355          shf(j,i) = shf_eb(j,i) / rho_cp
     1356
     1357          IF ( humidity )  THEN
     1358             qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt               &
     1359                             * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
     1360
     1361             qsws(j,i) = qsws_eb(j,i) / rho_lv
     1362
     1363             qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                   &
     1364                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
     1365                                 * t_surface_p(j,i) )
     1366
     1367             qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                   &
     1368                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
     1369                                 * t_surface_p(j,i) )
     1370
     1371             qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                   &
     1372                                 + dq_s_dt * t_surface(j,i) - dq_s_dt          &
     1373                                 * t_surface_p(j,i) )
     1374          ENDIF
     1375
     1376!
     1377!--       Calculate the true surface resistance
     1378          IF ( qsws_eb(j,i) == 0.0_wp )  THEN
     1379             r_s(j,i) = 1.0E10_wp
     1380          ELSE
     1381             r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
     1382                        * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )        &
     1383                        / qsws_eb(j,i) - r_a(j,i)
     1384          ENDIF
     1385
     1386!
     1387!--       Calculate change in liquid water reservoir due to dew fall or
     1388!--       evaporation of liquid water
     1389          IF ( humidity )  THEN
     1390!
     1391!--          If precipitation is activated, add rain water to qsws_liq_eb
     1392!--          and qsws_soil_eb according the the vegetation coverage.
     1393!--          precipitation_rate is given in mm.
     1394             IF ( precipitation )  THEN
     1395
     1396!
     1397!--             Add precipitation to liquid water reservoir, if possible.
     1398!--             Otherwise, add the water to soil. In case of
     1399!--             pavements, the exceeding water amount is implicitely removed
     1400!--             as runoff as qsws_soil_eb is then not used in the soil model
     1401                IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
     1402                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
     1403                                    + c_veg(j,i) * prr(k,j,i) * hyrho(k)       &
     1404                                    * 0.001_wp * rho_l * l_v
     1405                ELSE
     1406                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
     1407                                     + c_veg(j,i) * prr(k,j,i) * hyrho(k)      &
     1408                                     * 0.001_wp * rho_l * l_v
     1409                ENDIF
     1410
     1411!--             Add precipitation to bare soil according to the bare soil
     1412!--             coverage.
     1413                qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp                &
     1414                                    - c_veg(j,i)) * prr(k,j,i) * hyrho(k)      &
     1415                                    * 0.001_wp * rho_l * l_v
    13681416             ENDIF
    13691417
    1370              pt(k,j,i) = t_surface_p(j,i) / exn
    1371 
    1372 !
    1373 !--          Calculate fluxes
    1374 #if defined ( __rrtmg )
    1375              rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)      &
    1376                                 * t_surface(j,i) - rad_lw_out(nzb,j,i)         &
    1377                                 - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
    1378 
    1379              IF ( rrtm_idrv == 1 )  THEN
    1380                 rad_net(j,i) = rad_net_l(j,i)
    1381                 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i)                      &
    1382                                       + rad_lw_out_change_0(j,i)               &
    1383                                       * ( t_surface_p(j,i) - t_surface(j,i) )
     1418!
     1419!--          If the air is saturated, check the reservoir water level
     1420             IF ( qsws_eb(j,i) < 0.0_wp )  THEN
     1421
     1422!
     1423!--             Check if reservoir is full (avoid values > m_liq_eb_max)
     1424!--             In that case, qsws_liq_eb goes to qsws_soil_eb. In this
     1425!--             case qsws_veg_eb is zero anyway (because c_liq = 1),       
     1426!--             so that tend is zero and no further check is needed
     1427                IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
     1428                   qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
     1429                                        + qsws_liq_eb(j,i)
     1430
     1431                   qsws_liq_eb(j,i)  = 0.0_wp
     1432                ENDIF
     1433
     1434!
     1435!--             In case qsws_veg_eb becomes negative (unphysical behavior),
     1436!--             let the water enter the liquid water reservoir as dew on the
     1437!--             plant
     1438                IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
     1439                   qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
     1440                   qsws_veg_eb(j,i) = 0.0_wp
     1441                ENDIF
     1442             ENDIF                   
     1443 
     1444             tend = - qsws_liq_eb(j,i) * drho_l_lv
     1445             m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend         &
     1446                                                + tsc(3) * tm_liq_eb_m(j,i) )
     1447
     1448!
     1449!--          Check if reservoir is overfull -> reduce to maximum
     1450!--          (conservation of water is violated here)
     1451             m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
     1452
     1453!
     1454!--          Check if reservoir is empty (avoid values < 0.0)
     1455!--          (conservation of water is violated here)
     1456             m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
     1457
     1458
     1459!
     1460!--          Calculate m_liq_eb tendencies for the next Runge-Kutta step
     1461             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1462                IF ( intermediate_timestep_count == 1 )  THEN
     1463                   tm_liq_eb_m(j,i) = tend
     1464                ELSEIF ( intermediate_timestep_count <                         &
     1465                         intermediate_timestep_count_max )  THEN
     1466                   tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
     1467                                    * tm_liq_eb_m(j,i)
     1468                ENDIF
    13841469             ENDIF
     1470
     1471          ENDIF
     1472
     1473       ENDDO
     1474    ENDDO
     1475
     1476!
     1477!-- Make a logical OR for all processes. Force radiation call if at
     1478!-- least one processor reached the threshold change in skin temperature
     1479    IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count       &
     1480         == intermediate_timestep_count_max-1 )  THEN
     1481#if defined( __parallel )
     1482       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     1483       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
     1484                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
    13851485#else
    1386              rad_net_l(j,i)   = rad_net_l(j,i) + 3.0_wp * sigma_sb             &
    1387                                 * t_surface(j,i)**4 - 4.0_wp * sigma_sb        &
    1388                                 * t_surface(j,i)**3 * t_surface_p(j,i)
     1486       force_radiation_call = force_radiation_call_l
    13891487#endif
    1390 
    1391              ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)               &
    1392                               - t_soil(nzb_soil,j,i))
    1393 
    1394              shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
    1395 
    1396              shf(j,i) = shf_eb(j,i) / rho_cp
    1397 
    1398              IF ( humidity )  THEN
    1399                 qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt            &
    1400                                 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
    1401 
    1402                 qsws(j,i) = qsws_eb(j,i) / rho_lv
    1403 
    1404                 qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                &
    1405                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    1406                                     * t_surface_p(j,i) )
    1407 
    1408                 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                &
    1409                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    1410                                     * t_surface_p(j,i) )
    1411 
    1412                 qsws_liq_eb(j,i)  = - f_qsws_liq  * ( qv1 - q_s                &
    1413                                     + dq_s_dt * t_surface(j,i) - dq_s_dt       &
    1414                                     * t_surface_p(j,i) )
    1415              ENDIF
    1416 
    1417 !
    1418 !--          Calculate the true surface resistance
    1419              IF ( qsws_eb(j,i) == 0.0_wp )  THEN
    1420                 r_s(j,i) = 1.0E10_wp
    1421              ELSE
    1422                 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                    &
    1423                            * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )     &
    1424                            / qsws_eb(j,i) - r_a(j,i)
    1425              ENDIF
    1426 
    1427 !
    1428 !--          Calculate change in liquid water reservoir due to dew fall or
    1429 !--          evaporation of liquid water
    1430              IF ( humidity )  THEN
    1431 !
    1432 !--             If precipitation is activated, add rain water to qsws_liq_eb
    1433 !--             and qsws_soil_eb according the the vegetation coverage.
    1434 !--             precipitation_rate is given in mm.
    1435                 IF ( precipitation )  THEN
    1436 
    1437 !
    1438 !--                Add precipitation to liquid water reservoir, if possible.
    1439 !--                Otherwise, add the water to soil. In case of
    1440 !--                pavements, the exceeding water amount is implicitely removed
    1441 !--                as runoff as qsws_soil_eb is then not used in the soil model
    1442                    IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
    1443                       qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                      &
    1444                                        + c_veg(j,i) * prr(k,j,i) * hyrho(k)    &
    1445                                        * 0.001_wp * rho_l * l_v
    1446                    ELSE
    1447                       qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
    1448                                         + c_veg(j,i) * prr(k,j,i) * hyrho(k)   &
    1449                                         * 0.001_wp * rho_l * l_v
    1450                    ENDIF
    1451 
    1452 !--                Add precipitation to bare soil according to the bare soil
    1453 !--                coverage.
    1454                    qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp             &
    1455                                        - c_veg(j,i)) * prr(k,j,i) * hyrho(k)   &
    1456                                        * 0.001_wp * rho_l * l_v
    1457                 ENDIF
    1458 
    1459 !
    1460 !--             If the air is saturated, check the reservoir water level
    1461                 IF ( qsws_eb(j,i) < 0.0_wp )  THEN
    1462 
    1463 !
    1464 !--                Check if reservoir is full (avoid values > m_liq_eb_max)
    1465 !--                In that case, qsws_liq_eb goes to qsws_soil_eb. In this
    1466 !--                case qsws_veg_eb is zero anyway (because c_liq = 1),       
    1467 !--                so that tend is zero and no further check is needed
    1468                    IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
    1469                       qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                    &
    1470                                            + qsws_liq_eb(j,i)
    1471 
    1472                       qsws_liq_eb(j,i)  = 0.0_wp
    1473                    ENDIF
    1474 
    1475 !
    1476 !--                In case qsws_veg_eb becomes negative (unphysical behavior),
    1477 !--                let the water enter the liquid water reservoir as dew on the
    1478 !--                plant
    1479                    IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
    1480                       qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
    1481                       qsws_veg_eb(j,i) = 0.0_wp
    1482                    ENDIF
    1483                 ENDIF                   
    1484  
    1485                 tend = - qsws_liq_eb(j,i) * drho_l_lv
    1486 
    1487                 m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend      &
    1488                                                    + tsc(3) * tm_liq_eb_m(j,i) )
    1489 
    1490 !
    1491 !--             Check if reservoir is overfull -> reduce to maximum
    1492 !--             (conservation of water is violated here)
    1493                 m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
    1494 
    1495 !
    1496 !--             Check if reservoir is empty (avoid values < 0.0)
    1497 !--             (conservation of water is violated here)
    1498                 m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
    1499 
    1500 
    1501 !
    1502 !--             Calculate m_liq_eb tendencies for the next Runge-Kutta step
    1503                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1504                    IF ( intermediate_timestep_count == 1 )  THEN
    1505                       tm_liq_eb_m(j,i) = tend
    1506                    ELSEIF ( intermediate_timestep_count <                      &
    1507                             intermediate_timestep_count_max )  THEN
    1508                       tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp         &
    1509                                       * tm_liq_eb_m(j,i)
    1510                    ENDIF
    1511                 ENDIF
    1512 
    1513              ENDIF
    1514 
    1515           ENDDO
    1516        ENDDO
    1517 
    1518 !
    1519 !--    Make a logical OR for all processes. Force radiation call if at
    1520 !--    least one processor reached the threshold change in skin temperature
    1521        IF ( unscheduled_radiation_calls  .AND.  intermediate_timestep_count    &
    1522             == intermediate_timestep_count_max-1 )  THEN
    1523 #if defined( __parallel )
    1524           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1525           CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
    1526                               1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
    1527 #else
    1528           force_radiation_call = force_radiation_call_l
    1529 #endif
    1530           force_radiation_call_l = .FALSE.
    1531        ENDIF
    1532 
    1533 !
    1534 !--    Calculate surface specific humidity
    1535        IF ( humidity )  THEN
    1536           CALL calc_q_surface
    1537        ENDIF
    1538 
    1539 !
    1540 !--    Calculate new roughness lengths (for water surfaces only)
    1541        CALL calc_z0_water_surface
    1542 
    1543 
    1544     END SUBROUTINE lsm_energy_balance
     1488       force_radiation_call_l = .FALSE.
     1489    ENDIF
     1490
     1491!
     1492!-- Calculate surface specific humidity
     1493    IF ( humidity )  THEN
     1494       CALL calc_q_surface
     1495    ENDIF
     1496
     1497!
     1498!-- Calculate new roughness lengths (for water surfaces only)
     1499    CALL calc_z0_water_surface
     1500
     1501
     1502 END SUBROUTINE lsm_energy_balance
     1503
    15451504
    15461505!------------------------------------------------------------------------------!
    … …  
    24512410! Description:
    24522411! ------------
    2453 !> Soubroutine for averaging 3D data
     2412!> Subroutine for averaging 3D data
    24542413!------------------------------------------------------------------------------!
    24552414SUBROUTINE lsm_3d_data_averaging( mode, variable )
    … …  
    28052764! Description:
    28062765! ------------
    2807 !> Soubroutine defines appropriate grid for netcdf variables.
     2766!> Subroutine defining appropriate grid for netcdf variables.
    28082767!> It is called out from subroutine netcdf.
    28092768!------------------------------------------------------------------------------!
    2810     SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     2769 SUBROUTINE lsm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    28112770   
    2812         IMPLICIT NONE
    2813 
    2814         CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
    2815         LOGICAL, INTENT(OUT)           ::  found       !<
    2816         CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
    2817         CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
    2818         CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
    2819 
    2820 
    2821 
    2822 !
    2823 !--     Check for the grid
    2824         SELECT CASE ( TRIM( var ) )
    2825 
    2826            CASE ( 'm_soil', 't_soil' )
    2827               found  = .TRUE.
    2828               grid_x = 'x'
    2829               grid_y = 'y'
    2830               grid_z = 'zs'
    2831 
    2832            CASE DEFAULT
    2833               found  = .FALSE.
    2834               grid_x = 'none'
    2835               grid_y = 'none'
    2836               grid_z = 'none'
    2837         END SELECT
    2838 
    2839     END SUBROUTINE lsm_define_netcdf_grid
     2771     IMPLICIT NONE
     2772
     2773     CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
     2774     LOGICAL, INTENT(OUT)           ::  found       !<
     2775     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     2776     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
     2777     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
     2778
     2779     found  = .TRUE.
     2780
     2781!
     2782!--  Check for the grid
     2783     SELECT CASE ( TRIM( var ) )
     2784
     2785        CASE ( 'm_soil', 't_soil', 'm_soil_xy', 't_soil_xy', 'm_soil_xz',      &
     2786               't_soil_xz', 'm_soil_yz', 't_soil_yz' )
     2787           grid_x = 'x'
     2788           grid_y = 'y'
     2789           grid_z = 'zs'
     2790
     2791        CASE DEFAULT
     2792           found  = .FALSE.
     2793           grid_x = 'none'
     2794           grid_y = 'none'
     2795           grid_z = 'none'
     2796     END SELECT
     2797
     2798 END SUBROUTINE lsm_define_netcdf_grid
    28402799
    28412800!------------------------------------------------------------------------------!
    … …  
    28432802! Description:
    28442803! ------------
    2845 !> Soubroutine defines 3D output variables
     2804!> Subroutine defining 3D output variables
    28462805!------------------------------------------------------------------------------!
    28472806 SUBROUTINE lsm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
    28482807                                two_d, nzb_do, nzt_do )
    28492808 
    2850 
    28512809    USE indices
    28522810
    28532811    USE kinds
    28542812
    2855     USE user
    28562813
    28572814    IMPLICIT NONE
    … …  
    31743131! Description:
    31753132! ------------
    3176 !> Soubroutine defines 3D output variables
     3133!> Subroutine defining 3D output variables
    31773134!------------------------------------------------------------------------------!
    31783135 SUBROUTINE lsm_data_output_3d( av, variable, found, local_pf )

  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r1973 r1976  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Removed remaining 2D land surface quantities. Definition of radiation
     22! quantities is now done directly in the respective module
    2223!
    2324! Former revisions:
    … …  
    148149!> @todo calculation of output time levels for parallel NetCDF still does not
    149150!>       cover every exception (change of dt_do, end_time in restart)
     151!> @todo timeseries and profile output still needs to be rewritten to allow
     152!>       modularization
    150153!------------------------------------------------------------------------------!
    151154 MODULE netcdf_interface
    … …  
    414417    USE profil_parameter,                                                      &
    415418        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
     419
     420    USE radiation_model_mod,                                                   &
     421        ONLY: radiation, radiation_define_netcdf_grid
    416422
    417423    USE spectra_mod,                                                           &
    … …  
    735741                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
    736742                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
    737                        'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr',            &
    738                        'rad_sw_hr', 'rho', 's', 'sa', 'vpt'  )
     743                       'rho', 's', 'sa', 'vpt'  )
    739744
    740745                   grid_x = 'x'
    … …  
    757762!
    758763!--             w grid
    759                 CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
    760                        'w' )
     764                CASE ( 'w' )
    761765
    762766                   grid_x = 'x'
    … …  
    772776                      CALL lsm_define_netcdf_grid( domask(mid,av,i), found,    &
    773777                                                   grid_x, grid_y, grid_z )
     778                   ENDIF
     779
     780!
     781!--                Check for radiation quantities
     782                   IF ( .NOT. found  .AND.  radiation )  THEN
     783                      CALL radiation_define_netcdf_grid( domask(mid,av,i),     &
     784                                                         found, grid_x, grid_y,&
     785                                                         grid_z )
    774786                   ENDIF
    775787
    … …  
    12391251                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
    12401252                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho', &
    1241                        's', 'sa', 'vpt' , 'rad_lw_cs_hr', 'rad_lw_hr',         &
    1242                        'rad_sw_cs_hr', 'rad_sw_hr' )
     1253                       's', 'sa', 'vpt' )
    12431254
    12441255                   grid_x = 'x'
    … …  
    12611272!
    12621273!--             w grid
    1263                 CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',   &
    1264                        'w' )
     1274                CASE ( 'w' )
    12651275
    12661276                   grid_x = 'x'
    … …  
    12751285                   IF ( land_surface )  THEN
    12761286                      CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
    1277                                                     grid_y, grid_z )
     1287                                                   grid_y, grid_z )
     1288                   ENDIF
     1289
     1290!
     1291!--                Check for radiation quantities
     1292                   IF ( .NOT. found  .AND.  radiation )  THEN
     1293                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
     1294                                                         grid_x, grid_y,       &
     1295                                                         grid_z )
    12781296                   ENDIF
    12791297                   
    … …  
    18371855                             'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
    18381856                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
    1839                              'qr_xy', 'qv_xy', 'rad_lw_cs_hr_xy',              &
    1840                              'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', 'rad_sw_hr_xy',&
    1841                              'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
     1857                             'qr_xy', 'qv_xy', 'rho_xy', 's_xy', 'sa_xy',      &
     1858                             'vpt_xy' )
    18421859
    18431860                         grid_x = 'x'
    … …  
    18601877!
    18611878!--                   w grid
    1862                       CASE ( 'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy',  &
    1863                              'rad_sw_out_xy' , 'w_xy' )
     1879                      CASE ( 'w_xy' )
    18641880
    18651881                         grid_x = 'x'
    18661882                         grid_y = 'y'
    18671883                         grid_z = 'zw'
    1868 !
    1869 !--                   soil grid
    1870                       CASE ( 'm_soil_xy', 't_soil_xy' )
    1871                          grid_x = 'x'
    1872                          grid_y = 'y'
    1873                          grid_z = 'zs'
     1884
    18741885
    18751886                      CASE DEFAULT
    18761887!
     1888!--                      Check for land surface quantities
     1889                         IF ( land_surface )  THEN
     1890                            CALL lsm_define_netcdf_grid( do2d(av,i), found,    &
     1891                                                   grid_x, grid_y, grid_z )
     1892                         ENDIF
     1893
     1894!
     1895!--                      Check for radiation quantities
     1896                         IF ( .NOT. found  .AND.  radiation )  THEN
     1897                            CALL radiation_define_netcdf_grid( do2d(av,i),     &
     1898                                                         found, grid_x, grid_y,&
     1899                                                         grid_z )
     1900                         ENDIF
     1901
     1902!
    18771903!--                      Check for user-defined quantities
    1878                          CALL user_define_netcdf_grid( do2d(av,i), found, &
    1879                                                        grid_x, grid_y, grid_z )
     1904                         IF ( .NOT. found )  THEN
     1905                            CALL user_define_netcdf_grid( do2d(av,i), found,   &
     1906                                                          grid_x, grid_y,      &
     1907                                                          grid_z )
     1908                         ENDIF
    18801909
    18811910                         IF ( .NOT. found )  THEN
    … …  
    24912520                          'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
    24922521                          'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
    2493                           'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',                   &
    2494                           'rad_sw_cs_hr_xz', 'rad_sw_hr_xz''rho_xz', 's_xz',   &
    2495                           'sa_xz', 'vpt_xz' )
     2522                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
    24962523
    24972524                      grid_x = 'x'
    … …  
    25142541!
    25152542!--                w grid
    2516                    CASE ( 'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz',     &
    2517                           'rad_sw_out_xz', 'w_xz' )
     2543                   CASE ( 'w_xz' )
    25182544
    25192545                      grid_x = 'x'
    … …  
    25212547                      grid_z = 'zw'
    25222548
    2523 !
    2524 !--                soil grid
    2525                    CASE ( 'm_soil_xz', 't_soil_xz' )
    2526 
    2527                       grid_x = 'x'
    2528                       grid_y = 'y'
    2529                       grid_z = 'zs'
    2530 
    25312549                   CASE DEFAULT
     2550
     2551!
     2552!--                   Check for land surface quantities
     2553                      IF ( land_surface )  THEN
     2554                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
     2555                                                      grid_x, grid_y, grid_z )
     2556                      ENDIF
     2557
     2558!
     2559!--                   Check for radiation quantities
     2560                      IF ( .NOT. found  .AND.  radiation )  THEN
     2561                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
     2562                                                            grid_x, grid_y,    &
     2563                                                            grid_z )
     2564                      ENDIF
     2565
    25322566!
    25332567!--                   Check for user-defined quantities
    2534                       CALL user_define_netcdf_grid( do2d(av,i), found, &
    2535                                                     grid_x, grid_y, grid_z )
     2568                      IF ( .NOT. found )  THEN
     2569                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
     2570                                                       grid_x, grid_y, grid_z )
     2571                      ENDIF
     2572
    25362573                      IF ( .NOT. found )  THEN
    25372574                         WRITE ( message_string, * ) 'no grid defined for', &
    … …  
    31373174                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',         &
    31383175                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
    3139                           'rad_lw_cs_hr_yz', 'rad_lw_hr_yz',                   &
    3140                           'rad_sw_cs_hr_yz', 'rad_sw_hr_yz''rho_yz', 's_yz',   &
    3141                           'sa_yz', 'vpt_yz' )
     3176                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
    31423177
    31433178                      grid_x = 'x'
    … …  
    31603195!
    31613196!--                w grid
    3162                    CASE ( 'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz',     &
    3163                           'rad_sw_out_yz', 'w_yz' )
     3197                   CASE ( 'w_yz' )
    31643198
    31653199                      grid_x = 'x'
    31663200                      grid_y = 'y'
    31673201                      grid_z = 'zw'
    3168 !
    3169 !--                soil grid
    3170                    CASE ( 'm_soil_yz', 't_soil_yz' )
    3171 
    3172                       grid_x = 'x'
    3173                       grid_y = 'y'
    3174                       grid_z = 'zs'
     3202
    31753203
    31763204                   CASE DEFAULT
    31773205!
     3206!--                   Check for land surface quantities
     3207                      IF ( land_surface )  THEN
     3208                         CALL lsm_define_netcdf_grid( do2d(av,i), found,       &
     3209                                                      grid_x, grid_y, grid_z )
     3210                      ENDIF
     3211
     3212!
     3213!--                   Check for radiation quantities
     3214                      IF ( .NOT. found  .AND.  radiation )  THEN
     3215                         CALL radiation_define_netcdf_grid( do2d(av,i), found, &
     3216                                                            grid_x, grid_y,    &
     3217                                                            grid_z )
     3218                      ENDIF
     3219
     3220!
    31783221!--                   Check for user-defined quantities
    3179                       CALL user_define_netcdf_grid( do2d(av,i), found, &
    3180                                                     grid_x, grid_y, grid_z )
     3222                      IF ( .NOT. found )  THEN
     3223                         CALL user_define_netcdf_grid( do2d(av,i), found,      &
     3224                                                       grid_x, grid_y, grid_z )
     3225                      ENDIF
    31813226
    31823227                      IF ( .NOT. found )  THEN
    … …  
    43134358!--                Check for user-defined quantities (found, grid_x and grid_y
    43144359!--                are dummies)
    4315                    CALL user_define_netcdf_grid( data_output_sp(i), found,  &
     4360                   CALL user_define_netcdf_grid( data_output_sp(i), found,     &
    43164361                                                 grid_x, grid_y, grid_z )
    43174362

  • palm/trunk/SOURCE/palm.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Added call to radiation_last_actions for binary output of land surface model
     22! data
    2223!
    2324! Former revisions:
    … …  
    134135!> machines are controlled via cpp-directives.
    135136!> Model runs are only feasible using the ksh-script mrun.
     137!>
     138!> @todo create routine last_actions instead of calling lsm_last_actions etc.
    136139!------------------------------------------------------------------------------!
    137140 PROGRAM palm
    … …  
    182185        ONLY:  cpl_id, nested_run, pmci_child_initialize, pmci_init,           &
    183186               pmci_modelconfiguration, pmci_parent_initialize
     187
     188    USE radiation_model_mod,                                                   &
     189        ONLY:  radiation, radiation_last_actions
    184190
    185191    USE statistics,                                                            &
    … …  
    376382    ENDIF
    377383
    378 
    379384!
    380385!-- Output of program header
    … …  
    394399       CALL data_output_2d( 'yz', 0 )
    395400    ENDIF
     401
    396402    IF ( do3d_at_begin )  THEN
    397403       CALL data_output_3d( 0 )
    … …  
    458464             CALL lsm_last_actions
    459465          ENDIF
     466          IF ( radiation )  THEN
     467             CALL radiation_last_actions
     468          ENDIF
    460469          CALL user_last_actions
    461470          IF ( write_binary(1:4) == 'true' )  CALL close_file( 14 )

  • palm/trunk/SOURCE/prognostic_equations.f90

    r1961 r1976  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Simplied calls to radiation model
    2222!
    2323! Former revisions:
    … …  
    305305
    306306    USE radiation_model_mod,                                                   &
    307         ONLY:  radiation, radiation_scheme, radiation_tendency,                &
     307        ONLY:  radiation, radiation_tendency,                                  &
    308308               skip_time_do_radiation
    309309
    … …  
    638638!
    639639!--          If required, add tendency due to radiative heating/cooling
    640              IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.            &
     640             IF ( radiation  .AND.                                             &
    641641                  simulated_time > skip_time_do_radiation )  THEN
    642642                CALL radiation_tendency ( i, j, tend )
    … …  
    13711371!
    13721372!--    If required, add tendency due to radiative heating/cooling
    1373        IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
     1373       IF ( radiation  .AND.                                                   &
    13741374            simulated_time > skip_time_do_radiation )  THEN
    13751375            CALL radiation_tendency ( tend )
    … …  
    22932293       ENDIF
    22942294
    2295        IF ( radiation .AND. radiation_scheme == 'rrtmg' .AND.                  &
     2295       IF ( radiation .AND.                                                    &
    22962296            simulated_time > skip_time_do_radiation )  THEN
    22972297            CALL radiation_tendency ( tend )

  • palm/trunk/SOURCE/radiation_model_mod.f90

    r1857 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Output of 2D/3D/masked data is now directly done within this module. The
     22! radiation schemes have been simplified for better usability so that
     23! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
     24! the radiation code used.
    2225!
    2326! Former revisions:
    … …  
    130133
    131134#if defined ( __rrtmg )
    132 
    133135    USE parrrsw,                                                               &
    134136        ONLY:  naerec, nbndsw
    … …  
    276278!
    277279!-- Flag parameters for RRTMGS (should not be changed)
    278     INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
     280    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
     281                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
    279282                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
    280283                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
    … …  
    289292    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
    290293                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
    291                     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)
    292                     rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
     294                    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)
    293295
    294296    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
    … …  
    385387    END INTERFACE radiation_constant
    386388 
     389    INTERFACE radiation_control
     390       MODULE PROCEDURE radiation_control
     391    END INTERFACE radiation_control
     392
     393    INTERFACE radiation_3d_data_averaging
     394       MODULE PROCEDURE radiation_3d_data_averaging
     395    END INTERFACE radiation_3d_data_averaging
     396
     397    INTERFACE radiation_data_output_2d
     398       MODULE PROCEDURE radiation_data_output_2d
     399    END INTERFACE radiation_data_output_2d
     400
     401    INTERFACE radiation_data_output_3d
     402       MODULE PROCEDURE radiation_data_output_3d
     403    END INTERFACE radiation_data_output_3d
     404
     405    INTERFACE radiation_data_output_mask
     406       MODULE PROCEDURE radiation_data_output_mask
     407    END INTERFACE radiation_data_output_mask
     408
     409    INTERFACE radiation_define_netcdf_grid
     410       MODULE PROCEDURE radiation_define_netcdf_grid
     411    END INTERFACE radiation_define_netcdf_grid
     412
    387413    INTERFACE radiation_header
    388414       MODULE PROCEDURE radiation_header
    … …  
    406432    END INTERFACE radiation_tendency
    407433
     434    INTERFACE radiation_read_restart_data
     435       MODULE PROCEDURE radiation_read_restart_data
     436    END INTERFACE radiation_read_restart_data
     437
     438    INTERFACE radiation_last_actions
     439       MODULE PROCEDURE radiation_last_actions
     440    END INTERFACE radiation_last_actions
     441
    408442    SAVE
    409443
    … …  
    411445
    412446!
    413 !-- Public functions
     447!-- Public functions / NEEDS SORTING
    414448    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
    415            radiation_check_parameters, radiation_clearsky, radiation_constant, &
    416            radiation_header, radiation_init, radiation_parin, radiation_rrtmg, &
    417            radiation_tendency
     449           radiation_check_parameters, radiation_control,                      &
     450           radiation_header, radiation_init, radiation_parin,                  &
     451           radiation_3d_data_averaging, radiation_tendency,                    &
     452           radiation_data_output_2d, radiation_data_output_3d,                 &
     453           radiation_define_netcdf_grid, radiation_last_actions,               &
     454           radiation_read_restart_data, radiation_data_output_mask
    418455   
    419456!
    420 !-- Public variables and constants
     457!-- Public variables and constants / NEEDS SORTING
    421458    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
    422459           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
    … …  
    424461           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
    425462           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
    426            rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
     463           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av,                           &
    427464           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
    428465
    429466
    430467#if defined ( __rrtmg )
    431     PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
     468    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
    432469#endif
    433470
    434471 CONTAINS
     472
     473
     474!------------------------------------------------------------------------------!
     475! Description:
     476! ------------
     477!> This subroutine controls the calls of the radiation schemes
     478!------------------------------------------------------------------------------!
     479    SUBROUTINE radiation_control
     480 
     481 
     482       IMPLICIT NONE
     483
     484
     485       SELECT CASE ( TRIM( radiation_scheme ) )
     486
     487          CASE ( 'constant' )
     488             CALL radiation_constant
     489         
     490          CASE ( 'clear-sky' )
     491             CALL radiation_clearsky
     492       
     493          CASE ( 'rrtmg' )
     494             CALL radiation_rrtmg
     495
     496          CASE DEFAULT
     497
     498       END SELECT
     499
     500
     501    END SUBROUTINE radiation_control
    435502
    436503!------------------------------------------------------------------------------!
    … …  
    456523       SELECT CASE ( TRIM( var ) )
    457524
    458          CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
    459                  'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
     525         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
    460526             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    461527                message_string = '"output of "' // TRIM( var ) // '" requi' // &
    … …  
    779845       ENDIF
    780846
    781 
    782        IF ( radiation_scheme == 'constant' )  THEN
    783 
    784           IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
    785              ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
    786           ENDIF
    787 
    788        ENDIF
    789 
    790        IF ( radiation_scheme == 'clear-sky' )  THEN
     847       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
     848            radiation_scheme == 'constant')  THEN
    791849
    792850          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
    … …  
    10751133                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
    10761134
     1135
     1136             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
     1137                                        * (pt(k,j,i) * exn) ** 3
     1138
    10771139          ENDDO
    10781140       ENDDO
    … …  
    10931155       INTEGER(iwp) :: i, j, k   !< loop indices
    10941156       REAL(wp)     :: exn,   &  !< Exner functions at surface
     1157                       exn1,  &  !< Exner functions at first grid level
    10951158                       pt1       !< potential temperature at first grid level
    10961159
    … …  
    10991162       exn = (surface_pressure / 1000.0_wp )**0.286_wp
    11001163!
    1101 !--    Prescribe net radiation and estimate a longwave outgoing radiative
    1102 !--    flux (needed in land surface model)
     1164!--    Prescribe net radiation and estimate the remaining radiative fluxes
    11031165       DO i = nxlg, nxrg
    11041166          DO j = nysg, nyng
    … …  
    11061168
    11071169             rad_net(j,i)      = net_radiation
     1170
     1171             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
     1172
     1173             IF ( cloud_physics )  THEN
     1174                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
     1175                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
     1176             ELSE
     1177                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
     1178             ENDIF
     1179
    11081180             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
     1181
     1182             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
     1183                                  + rad_lw_out(0,j,i) )                        &
     1184                                  / ( 1.0_wp - alpha(j,i) )
    11091185
    11101186          ENDDO
    … …  
    21472223!> Cache-optimized version.
    21482224!------------------------------------------------------------------------------!
    2149     SUBROUTINE radiation_tendency_ij ( i, j, tend )
    2150 
    2151        USE cloud_parameters,                                                   &
    2152            ONLY:  pt_d_t
    2153 
    2154        IMPLICIT NONE
    2155 
    2156        INTEGER(iwp) :: i, j, k !< loop indices
    2157 
    2158        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
    2159 
    2160 #if defined ( __rrtmg )
     2225 SUBROUTINE radiation_tendency_ij ( i, j, tend )
     2226
     2227    USE cloud_parameters,                                                      &
     2228        ONLY:  pt_d_t
     2229
     2230    IMPLICIT NONE
     2231
     2232    INTEGER(iwp) :: i, j, k !< loop indices
     2233
     2234    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
     2235
     2236    IF ( radiation_scheme == 'rrtmg' )  THEN
     2237#if defined  ( __rrtmg )
    21612238!
    21622239!--    Calculate tendency based on heating rate
    21632240       DO k = nzb+1, nzt+1
    21642241          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
    2165                                       * pt_d_t(k) * d_seconds_hour
     2242                                         * pt_d_t(k) * d_seconds_hour
    21662243       ENDDO
    2167 
    21682244#endif
     2245    ENDIF
    21692246
    21702247    END SUBROUTINE radiation_tendency_ij
    … …  
    21772254!> Vector-optimized version
    21782255!------------------------------------------------------------------------------!
    2179     SUBROUTINE radiation_tendency ( tend )
    2180 
    2181        USE cloud_parameters,                                                   &
    2182            ONLY:  pt_d_t
    2183 
    2184        USE indices,                                                            &
    2185            ONLY:  nxl, nxr, nyn, nys
    2186 
    2187        IMPLICIT NONE
    2188 
    2189        INTEGER(iwp) :: i, j, k !< loop indices
    2190 
    2191        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
    2192 
    2193 #if defined ( __rrtmg )
     2256 SUBROUTINE radiation_tendency ( tend )
     2257
     2258    USE cloud_parameters,                                                      &
     2259        ONLY:  pt_d_t
     2260
     2261    USE indices,                                                               &
     2262        ONLY:  nxl, nxr, nyn, nys
     2263
     2264    IMPLICIT NONE
     2265
     2266    INTEGER(iwp) :: i, j, k !< loop indices
     2267
     2268    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
     2269
     2270    IF ( radiation_scheme == 'rrtmg' )  THEN
     2271#if defined  ( __rrtmg )
    21942272!
    21952273!--    Calculate tendency based on heating rate
    … …  
    21982276             DO k = nzb+1, nzt+1
    21992277                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
    2200                                             +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
    2201                                             * d_seconds_hour
    2202              ENDDO
    2203          ENDDO
     2278                                          +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
     2279                                          * d_seconds_hour
     2280             ENDDO
     2281          ENDDO
    22042282       ENDDO
    22052283#endif
    2206 
    2207     END SUBROUTINE radiation_tendency
     2284    ENDIF
     2285
     2286
     2287 END SUBROUTINE radiation_tendency
     2288
     2289!------------------------------------------------------------------------------!
     2290!
     2291! Description:
     2292! ------------
     2293!> Subroutine for averaging 3D data
     2294!------------------------------------------------------------------------------!
     2295SUBROUTINE radiation_3d_data_averaging( mode, variable )
     2296 
     2297
     2298    USE control_parameters
     2299
     2300    USE indices
     2301
     2302    USE kinds
     2303
     2304    IMPLICIT NONE
     2305
     2306    CHARACTER (LEN=*) ::  mode    !<
     2307    CHARACTER (LEN=*) :: variable !<
     2308
     2309    INTEGER(iwp) ::  i !<
     2310    INTEGER(iwp) ::  j !<
     2311    INTEGER(iwp) ::  k !<
     2312
     2313    IF ( mode == 'allocate' )  THEN
     2314
     2315       SELECT CASE ( TRIM( variable ) )
     2316
     2317             CASE ( 'rad_net*' )
     2318                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
     2319                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
     2320                ENDIF
     2321                rad_net_av = 0.0_wp
     2322
     2323             CASE ( 'rad_lw_in' )
     2324                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
     2325                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2326                ENDIF
     2327                rad_lw_in_av = 0.0_wp
     2328
     2329             CASE ( 'rad_lw_out' )
     2330                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
     2331                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2332                ENDIF
     2333                rad_lw_out_av = 0.0_wp
     2334
     2335             CASE ( 'rad_lw_cs_hr' )
     2336                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
     2337                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     2338                ENDIF
     2339                rad_lw_cs_hr_av = 0.0_wp
     2340
     2341             CASE ( 'rad_lw_hr' )
     2342                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
     2343                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     2344                ENDIF
     2345                rad_lw_hr_av = 0.0_wp
     2346
     2347             CASE ( 'rad_sw_in' )
     2348                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
     2349                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2350                ENDIF
     2351                rad_sw_in_av = 0.0_wp
     2352
     2353             CASE ( 'rad_sw_out' )
     2354                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
     2355                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     2356                ENDIF
     2357                rad_sw_out_av = 0.0_wp
     2358
     2359             CASE ( 'rad_sw_cs_hr' )
     2360                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
     2361                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     2362                ENDIF
     2363                rad_sw_cs_hr_av = 0.0_wp
     2364
     2365             CASE ( 'rad_sw_hr' )
     2366                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
     2367                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
     2368                ENDIF
     2369                rad_sw_hr_av = 0.0_wp
     2370
     2371          CASE DEFAULT
     2372             CONTINUE
     2373
     2374       END SELECT
     2375
     2376    ELSEIF ( mode == 'sum' )  THEN
     2377
     2378       SELECT CASE ( TRIM( variable ) )
     2379
     2380          CASE ( 'rad_net*' )
     2381             DO  i = nxlg, nxrg
     2382                DO  j = nysg, nyng
     2383                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
     2384                ENDDO
     2385             ENDDO
     2386
     2387          CASE ( 'rad_lw_in' )
     2388             DO  i = nxlg, nxrg
     2389                DO  j = nysg, nyng
     2390                   DO  k = nzb, nzt+1
     2391                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
     2392                   ENDDO
     2393                ENDDO
     2394             ENDDO
     2395
     2396          CASE ( 'rad_lw_out' )
     2397             DO  i = nxlg, nxrg
     2398                DO  j = nysg, nyng
     2399                   DO  k = nzb, nzt+1
     2400                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
     2401                   ENDDO
     2402                ENDDO
     2403             ENDDO
     2404
     2405          CASE ( 'rad_lw_cs_hr' )
     2406             DO  i = nxlg, nxrg
     2407                DO  j = nysg, nyng
     2408                   DO  k = nzb, nzt+1
     2409                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
     2410                   ENDDO
     2411                ENDDO
     2412             ENDDO
     2413
     2414          CASE ( 'rad_lw_hr' )
     2415             DO  i = nxlg, nxrg
     2416                DO  j = nysg, nyng
     2417                   DO  k = nzb, nzt+1
     2418                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
     2419                   ENDDO
     2420                ENDDO
     2421             ENDDO
     2422
     2423          CASE ( 'rad_sw_in' )
     2424             DO  i = nxlg, nxrg
     2425                DO  j = nysg, nyng
     2426                   DO  k = nzb, nzt+1
     2427                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
     2428                   ENDDO
     2429                ENDDO
     2430             ENDDO
     2431
     2432          CASE ( 'rad_sw_out' )
     2433             DO  i = nxlg, nxrg
     2434                DO  j = nysg, nyng
     2435                   DO  k = nzb, nzt+1
     2436                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
     2437                   ENDDO
     2438                ENDDO
     2439             ENDDO
     2440
     2441          CASE ( 'rad_sw_cs_hr' )
     2442             DO  i = nxlg, nxrg
     2443                DO  j = nysg, nyng
     2444                   DO  k = nzb, nzt+1
     2445                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
     2446                   ENDDO
     2447                ENDDO
     2448             ENDDO
     2449
     2450          CASE ( 'rad_sw_hr' )
     2451             DO  i = nxlg, nxrg
     2452                DO  j = nysg, nyng
     2453                   DO  k = nzb, nzt+1
     2454                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
     2455                   ENDDO
     2456                ENDDO
     2457             ENDDO
     2458
     2459          CASE DEFAULT
     2460             CONTINUE
     2461
     2462       END SELECT
     2463
     2464    ELSEIF ( mode == 'average' )  THEN
     2465
     2466       SELECT CASE ( TRIM( variable ) )
     2467
     2468         CASE ( 'rad_net*' )
     2469             DO  i = nxlg, nxrg
     2470                DO  j = nysg, nyng
     2471                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
     2472                ENDDO
     2473             ENDDO
     2474
     2475          CASE ( 'rad_lw_in' )
     2476             DO  i = nxlg, nxrg
     2477                DO  j = nysg, nyng
     2478                   DO  k = nzb, nzt+1
     2479                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2480                   ENDDO
     2481                ENDDO
     2482             ENDDO
     2483
     2484          CASE ( 'rad_lw_out' )
     2485             DO  i = nxlg, nxrg
     2486                DO  j = nysg, nyng
     2487                   DO  k = nzb, nzt+1
     2488                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2489                   ENDDO
     2490                ENDDO
     2491             ENDDO
     2492
     2493          CASE ( 'rad_lw_cs_hr' )
     2494             DO  i = nxlg, nxrg
     2495                DO  j = nysg, nyng
     2496                   DO  k = nzb, nzt+1
     2497                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2498                   ENDDO
     2499                ENDDO
     2500             ENDDO
     2501
     2502          CASE ( 'rad_lw_hr' )
     2503             DO  i = nxlg, nxrg
     2504                DO  j = nysg, nyng
     2505                   DO  k = nzb, nzt+1
     2506                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2507                   ENDDO
     2508                ENDDO
     2509             ENDDO
     2510
     2511          CASE ( 'rad_sw_in' )
     2512             DO  i = nxlg, nxrg
     2513                DO  j = nysg, nyng
     2514                   DO  k = nzb, nzt+1
     2515                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2516                   ENDDO
     2517                ENDDO
     2518             ENDDO
     2519
     2520          CASE ( 'rad_sw_out' )
     2521             DO  i = nxlg, nxrg
     2522                DO  j = nysg, nyng
     2523                   DO  k = nzb, nzt+1
     2524                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2525                   ENDDO
     2526                ENDDO
     2527             ENDDO
     2528
     2529          CASE ( 'rad_sw_cs_hr' )
     2530             DO  i = nxlg, nxrg
     2531                DO  j = nysg, nyng
     2532                   DO  k = nzb, nzt+1
     2533                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2534                   ENDDO
     2535                ENDDO
     2536             ENDDO
     2537
     2538          CASE ( 'rad_sw_hr' )
     2539             DO  i = nxlg, nxrg
     2540                DO  j = nysg, nyng
     2541                   DO  k = nzb, nzt+1
     2542                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2543                   ENDDO
     2544                ENDDO
     2545             ENDDO
     2546
     2547       END SELECT
     2548
     2549    ENDIF
     2550
     2551END SUBROUTINE radiation_3d_data_averaging
     2552
     2553
     2554!------------------------------------------------------------------------------!
     2555!
     2556! Description:
     2557! ------------
     2558!> Subroutine defining appropriate grid for netcdf variables.
     2559!> It is called out from subroutine netcdf.
     2560!------------------------------------------------------------------------------!
     2561SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
     2562   
     2563    IMPLICIT NONE
     2564
     2565    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
     2566    LOGICAL, INTENT(OUT)           ::  found       !<
     2567    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
     2568    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
     2569    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
     2570
     2571    found  = .TRUE.
     2572
     2573
     2574!
     2575!-- Check for the grid
     2576    SELECT CASE ( TRIM( var ) )
     2577
     2578       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
     2579              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
     2580              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
     2581              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
     2582              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
     2583          grid_x = 'x'
     2584          grid_y = 'y'
     2585          grid_z = 'zu'
     2586
     2587       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
     2588              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
     2589              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
     2590              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
     2591          grid_x = 'x'
     2592          grid_y = 'y'
     2593          grid_z = 'zw'
     2594
     2595
     2596       CASE DEFAULT
     2597          found  = .FALSE.
     2598          grid_x = 'none'
     2599          grid_y = 'none'
     2600          grid_z = 'none'
     2601
     2602        END SELECT
     2603
     2604    END SUBROUTINE radiation_define_netcdf_grid
     2605
     2606!------------------------------------------------------------------------------!
     2607!
     2608! Description:
     2609! ------------
     2610!> Subroutine defining 3D output variables
     2611!------------------------------------------------------------------------------!
     2612 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
     2613                                      local_pf, two_d )
     2614 
     2615    USE indices
     2616
     2617    USE kinds
     2618
     2619
     2620    IMPLICIT NONE
     2621
     2622    CHARACTER (LEN=*) ::  grid     !<
     2623    CHARACTER (LEN=*) ::  mode     !<
     2624    CHARACTER (LEN=*) ::  variable !<
     2625
     2626    INTEGER(iwp) ::  av !<
     2627    INTEGER(iwp) ::  i  !<
     2628    INTEGER(iwp) ::  j  !<
     2629    INTEGER(iwp) ::  k  !<
     2630
     2631    LOGICAL      ::  found !<
     2632    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
     2633
     2634    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
     2635
     2636    found = .TRUE.
     2637
     2638    SELECT CASE ( TRIM( variable ) )
     2639
     2640       CASE ( 'rad_net*_xy' )        ! 2d-array
     2641          IF ( av == 0 ) THEN
     2642             DO  i = nxlg, nxrg
     2643                DO  j = nysg, nyng
     2644                   local_pf(i,j,nzb+1) = rad_net(j,i)
     2645                ENDDO
     2646             ENDDO
     2647          ELSE
     2648             DO  i = nxlg, nxrg
     2649                DO  j = nysg, nyng
     2650                   local_pf(i,j,nzb+1) = rad_net_av(j,i)
     2651                ENDDO
     2652             ENDDO
     2653          ENDIF
     2654          two_d = .TRUE.
     2655          grid = 'zu1'
     2656
     2657 
     2658       CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
     2659          IF ( av == 0 ) THEN
     2660             DO  i = nxlg, nxrg
     2661                DO  j = nysg, nyng
     2662                   DO  k = nzb, nzt+1
     2663                      local_pf(i,j,k) = rad_lw_in(k,j,i)
     2664                   ENDDO
     2665                ENDDO
     2666             ENDDO
     2667          ELSE
     2668             DO  i = nxlg, nxrg
     2669                DO  j = nysg, nyng
     2670                   DO  k = nzb, nzt+1
     2671                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
     2672                   ENDDO
     2673                ENDDO
     2674             ENDDO
     2675          ENDIF
     2676          IF ( mode == 'xy' )  grid = 'zu'
     2677
     2678       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
     2679          IF ( av == 0 ) THEN
     2680             DO  i = nxlg, nxrg
     2681                DO  j = nysg, nyng
     2682                   DO  k = nzb, nzt+1
     2683                      local_pf(i,j,k) = rad_lw_out(k,j,i)
     2684                   ENDDO
     2685                ENDDO
     2686             ENDDO
     2687          ELSE
     2688             DO  i = nxlg, nxrg
     2689                DO  j = nysg, nyng
     2690                   DO  k = nzb, nzt+1
     2691                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
     2692                   ENDDO
     2693                ENDDO
     2694             ENDDO
     2695          ENDIF   
     2696          IF ( mode == 'xy' )  grid = 'zu'
     2697
     2698       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
     2699          IF ( av == 0 ) THEN
     2700             DO  i = nxlg, nxrg
     2701                DO  j = nysg, nyng
     2702                   DO  k = nzb, nzt+1
     2703                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
     2704                   ENDDO
     2705                ENDDO
     2706             ENDDO
     2707          ELSE
     2708             DO  i = nxlg, nxrg
     2709                DO  j = nysg, nyng
     2710                   DO  k = nzb, nzt+1
     2711                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
     2712                   ENDDO
     2713                ENDDO
     2714             ENDDO
     2715          ENDIF
     2716          IF ( mode == 'xy' )  grid = 'zw'
     2717
     2718       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
     2719          IF ( av == 0 ) THEN
     2720             DO  i = nxlg, nxrg
     2721                DO  j = nysg, nyng
     2722                   DO  k = nzb, nzt+1
     2723                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
     2724                   ENDDO
     2725                ENDDO
     2726             ENDDO
     2727          ELSE
     2728             DO  i = nxlg, nxrg
     2729                DO  j = nysg, nyng
     2730                   DO  k = nzb, nzt+1
     2731                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
     2732                   ENDDO
     2733                ENDDO
     2734             ENDDO
     2735          ENDIF
     2736          IF ( mode == 'xy' )  grid = 'zw'
     2737
     2738       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
     2739          IF ( av == 0 ) THEN
     2740             DO  i = nxlg, nxrg
     2741                DO  j = nysg, nyng
     2742                   DO  k = nzb, nzt+1
     2743                      local_pf(i,j,k) = rad_sw_in(k,j,i)
     2744                   ENDDO
     2745                ENDDO
     2746             ENDDO
     2747          ELSE
     2748             DO  i = nxlg, nxrg
     2749                DO  j = nysg, nyng
     2750                   DO  k = nzb, nzt+1
     2751                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
     2752                   ENDDO
     2753                ENDDO
     2754             ENDDO
     2755          ENDIF
     2756          IF ( mode == 'xy' )  grid = 'zu'
     2757
     2758       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
     2759          IF ( av == 0 ) THEN
     2760             DO  i = nxlg, nxrg
     2761                DO  j = nysg, nyng
     2762                   DO  k = nzb, nzt+1
     2763                      local_pf(i,j,k) = rad_sw_out(k,j,i)
     2764                   ENDDO
     2765                ENDDO
     2766             ENDDO
     2767          ELSE
     2768             DO  i = nxlg, nxrg
     2769                DO  j = nysg, nyng
     2770                   DO  k = nzb, nzt+1
     2771                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
     2772                   ENDDO
     2773                ENDDO
     2774             ENDDO
     2775          ENDIF
     2776          IF ( mode == 'xy' )  grid = 'zu'
     2777
     2778       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
     2779          IF ( av == 0 ) THEN
     2780             DO  i = nxlg, nxrg
     2781                DO  j = nysg, nyng
     2782                   DO  k = nzb, nzt+1
     2783                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
     2784                   ENDDO
     2785                ENDDO
     2786             ENDDO
     2787          ELSE
     2788             DO  i = nxlg, nxrg
     2789                DO  j = nysg, nyng
     2790                   DO  k = nzb, nzt+1
     2791                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
     2792                   ENDDO
     2793                ENDDO
     2794             ENDDO
     2795          ENDIF
     2796          IF ( mode == 'xy' )  grid = 'zw'
     2797
     2798       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
     2799          IF ( av == 0 ) THEN
     2800             DO  i = nxlg, nxrg
     2801                DO  j = nysg, nyng
     2802                   DO  k = nzb, nzt+1
     2803                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
     2804                   ENDDO
     2805                ENDDO
     2806             ENDDO
     2807          ELSE
     2808             DO  i = nxlg, nxrg
     2809                DO  j = nysg, nyng
     2810                   DO  k = nzb, nzt+1
     2811                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
     2812                   ENDDO
     2813                ENDDO
     2814             ENDDO
     2815          ENDIF
     2816          IF ( mode == 'xy' )  grid = 'zw'
     2817
     2818       CASE DEFAULT
     2819          found = .FALSE.
     2820          grid  = 'none'
     2821
     2822    END SELECT
     2823 
     2824 END SUBROUTINE radiation_data_output_2d
     2825
     2826
     2827!------------------------------------------------------------------------------!
     2828!
     2829! Description:
     2830! ------------
     2831!> Subroutine defining 3D output variables
     2832!------------------------------------------------------------------------------!
     2833 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
     2834 
     2835
     2836    USE indices
     2837
     2838    USE kinds
     2839
     2840
     2841    IMPLICIT NONE
     2842
     2843    CHARACTER (LEN=*) ::  variable !<
     2844
     2845    INTEGER(iwp) ::  av    !<
     2846    INTEGER(iwp) ::  i     !<
     2847    INTEGER(iwp) ::  j     !<
     2848    INTEGER(iwp) ::  k     !<
     2849
     2850    LOGICAL      ::  found !<
     2851
     2852    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
     2853
     2854
     2855    found = .TRUE.
     2856
     2857
     2858    SELECT CASE ( TRIM( variable ) )
     2859
     2860      CASE ( 'rad_sw_in' )
     2861         IF ( av == 0 )  THEN
     2862            DO  i = nxlg, nxrg
     2863               DO  j = nysg, nyng
     2864                  DO  k = nzb, nzt+1
     2865                     local_pf(i,j,k) = rad_sw_in(k,j,i)
     2866                  ENDDO
     2867               ENDDO
     2868            ENDDO
     2869         ELSE
     2870            DO  i = nxlg, nxrg
     2871               DO  j = nysg, nyng
     2872                  DO  k = nzb, nzt+1
     2873                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
     2874                  ENDDO
     2875               ENDDO
     2876            ENDDO
     2877         ENDIF
     2878
     2879      CASE ( 'rad_sw_out' )
     2880         IF ( av == 0 )  THEN
     2881            DO  i = nxlg, nxrg
     2882               DO  j = nysg, nyng
     2883                  DO  k = nzb, nzt+1
     2884                     local_pf(i,j,k) = rad_sw_out(k,j,i)
     2885                  ENDDO
     2886               ENDDO
     2887            ENDDO
     2888         ELSE
     2889            DO  i = nxlg, nxrg
     2890               DO  j = nysg, nyng
     2891                  DO  k = nzb, nzt+1
     2892                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
     2893                  ENDDO
     2894               ENDDO
     2895            ENDDO
     2896         ENDIF
     2897
     2898      CASE ( 'rad_sw_cs_hr' )
     2899         IF ( av == 0 )  THEN
     2900            DO  i = nxlg, nxrg
     2901               DO  j = nysg, nyng
     2902                  DO  k = nzb, nzt+1
     2903                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
     2904                  ENDDO
     2905               ENDDO
     2906            ENDDO
     2907         ELSE
     2908            DO  i = nxlg, nxrg
     2909               DO  j = nysg, nyng
     2910                  DO  k = nzb, nzt+1
     2911                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
     2912                  ENDDO
     2913               ENDDO
     2914            ENDDO
     2915         ENDIF
     2916
     2917      CASE ( 'rad_sw_hr' )
     2918         IF ( av == 0 )  THEN
     2919            DO  i = nxlg, nxrg
     2920               DO  j = nysg, nyng
     2921                  DO  k = nzb, nzt+1
     2922                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
     2923                  ENDDO
     2924               ENDDO
     2925            ENDDO
     2926         ELSE
     2927            DO  i = nxlg, nxrg
     2928               DO  j = nysg, nyng
     2929                  DO  k = nzb, nzt+1
     2930                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
     2931                  ENDDO
     2932               ENDDO
     2933            ENDDO
     2934         ENDIF
     2935
     2936      CASE ( 'rad_lw_in' )
     2937         IF ( av == 0 )  THEN
     2938            DO  i = nxlg, nxrg
     2939               DO  j = nysg, nyng
     2940                  DO  k = nzb, nzt+1
     2941                     local_pf(i,j,k) = rad_lw_in(k,j,i)
     2942                  ENDDO
     2943               ENDDO
     2944            ENDDO
     2945         ELSE
     2946            DO  i = nxlg, nxrg
     2947               DO  j = nysg, nyng
     2948                  DO  k = nzb, nzt+1
     2949                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
     2950                  ENDDO
     2951               ENDDO
     2952            ENDDO
     2953         ENDIF
     2954
     2955      CASE ( 'rad_lw_out' )
     2956         IF ( av == 0 )  THEN
     2957            DO  i = nxlg, nxrg
     2958               DO  j = nysg, nyng
     2959                  DO  k = nzb, nzt+1
     2960                     local_pf(i,j,k) = rad_lw_out(k,j,i)
     2961                  ENDDO
     2962               ENDDO
     2963            ENDDO
     2964         ELSE
     2965            DO  i = nxlg, nxrg
     2966               DO  j = nysg, nyng
     2967                  DO  k = nzb, nzt+1
     2968                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
     2969                  ENDDO
     2970               ENDDO
     2971            ENDDO
     2972         ENDIF
     2973
     2974      CASE ( 'rad_lw_cs_hr' )
     2975         IF ( av == 0 )  THEN
     2976            DO  i = nxlg, nxrg
     2977               DO  j = nysg, nyng
     2978                  DO  k = nzb, nzt+1
     2979                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
     2980                  ENDDO
     2981               ENDDO
     2982            ENDDO
     2983         ELSE
     2984            DO  i = nxlg, nxrg
     2985               DO  j = nysg, nyng
     2986                  DO  k = nzb, nzt+1
     2987                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
     2988                  ENDDO
     2989               ENDDO
     2990            ENDDO
     2991         ENDIF
     2992
     2993      CASE ( 'rad_lw_hr' )
     2994         IF ( av == 0 )  THEN
     2995            DO  i = nxlg, nxrg
     2996               DO  j = nysg, nyng
     2997                  DO  k = nzb, nzt+1
     2998                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
     2999                  ENDDO
     3000               ENDDO
     3001            ENDDO
     3002         ELSE
     3003            DO  i = nxlg, nxrg
     3004               DO  j = nysg, nyng
     3005                  DO  k = nzb, nzt+1
     3006                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
     3007                  ENDDO
     3008               ENDDO
     3009            ENDDO
     3010         ENDIF
     3011
     3012       CASE DEFAULT
     3013          found = .FALSE.
     3014
     3015    END SELECT
     3016
     3017
     3018 END SUBROUTINE radiation_data_output_3d
     3019
     3020!------------------------------------------------------------------------------!
     3021!
     3022! Description:
     3023! ------------
     3024!> Subroutine defining masked data output
     3025!------------------------------------------------------------------------------!
     3026 SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
     3027 
     3028    USE control_parameters
     3029       
     3030    USE indices
     3031   
     3032    USE kinds
     3033   
     3034
     3035    IMPLICIT NONE
     3036
     3037    CHARACTER (LEN=*) ::  variable   !<
     3038
     3039    INTEGER(iwp) ::  av   !<
     3040    INTEGER(iwp) ::  i    !<
     3041    INTEGER(iwp) ::  j    !<
     3042    INTEGER(iwp) ::  k    !<
     3043
     3044    LOGICAL ::  found     !<
     3045
     3046    REAL(wp),                                                                  &
     3047       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
     3048          local_pf   !<
     3049
     3050
     3051    found = .TRUE.
     3052
     3053    SELECT CASE ( TRIM( variable ) )
     3054
     3055
     3056       CASE ( 'rad_lw_in' )
     3057          IF ( av == 0 )  THEN
     3058             DO  i = 1, mask_size_l(mid,1)
     3059                DO  j = 1, mask_size_l(mid,2)
     3060                   DO  k = 1, mask_size_l(mid,3)
     3061                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
     3062                                            mask_j(mid,j),mask_i(mid,i))
     3063                    ENDDO
     3064                 ENDDO
     3065              ENDDO
     3066          ELSE
     3067             DO  i = 1, mask_size_l(mid,1)
     3068                DO  j = 1, mask_size_l(mid,2)
     3069                   DO  k = 1, mask_size_l(mid,3)
     3070                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
     3071                                               mask_j(mid,j),mask_i(mid,i))
     3072                   ENDDO
     3073                ENDDO
     3074             ENDDO
     3075          ENDIF
     3076
     3077       CASE ( 'rad_lw_out' )
     3078          IF ( av == 0 )  THEN
     3079             DO  i = 1, mask_size_l(mid,1)
     3080                DO  j = 1, mask_size_l(mid,2)
     3081                   DO  k = 1, mask_size_l(mid,3)
     3082                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
     3083                                            mask_j(mid,j),mask_i(mid,i))
     3084                    ENDDO
     3085                 ENDDO
     3086              ENDDO
     3087          ELSE
     3088             DO  i = 1, mask_size_l(mid,1)
     3089                DO  j = 1, mask_size_l(mid,2)
     3090                   DO  k = 1, mask_size_l(mid,3)
     3091                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
     3092                                               mask_j(mid,j),mask_i(mid,i))
     3093                   ENDDO
     3094                ENDDO
     3095             ENDDO
     3096          ENDIF
     3097
     3098       CASE ( 'rad_lw_cs_hr' )
     3099          IF ( av == 0 )  THEN
     3100             DO  i = 1, mask_size_l(mid,1)
     3101                DO  j = 1, mask_size_l(mid,2)
     3102                   DO  k = 1, mask_size_l(mid,3)
     3103                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
     3104                                            mask_j(mid,j),mask_i(mid,i))
     3105                    ENDDO
     3106                 ENDDO
     3107              ENDDO
     3108          ELSE
     3109             DO  i = 1, mask_size_l(mid,1)
     3110                DO  j = 1, mask_size_l(mid,2)
     3111                   DO  k = 1, mask_size_l(mid,3)
     3112                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
     3113                                               mask_j(mid,j),mask_i(mid,i))
     3114                   ENDDO
     3115                ENDDO
     3116             ENDDO
     3117          ENDIF
     3118
     3119       CASE ( 'rad_lw_hr' )
     3120          IF ( av == 0 )  THEN
     3121             DO  i = 1, mask_size_l(mid,1)
     3122                DO  j = 1, mask_size_l(mid,2)
     3123                   DO  k = 1, mask_size_l(mid,3)
     3124                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
     3125                                            mask_j(mid,j),mask_i(mid,i))
     3126                    ENDDO
     3127                 ENDDO
     3128              ENDDO
     3129          ELSE
     3130             DO  i = 1, mask_size_l(mid,1)
     3131                DO  j = 1, mask_size_l(mid,2)
     3132                   DO  k = 1, mask_size_l(mid,3)
     3133                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
     3134                                               mask_j(mid,j),mask_i(mid,i))
     3135                   ENDDO
     3136                ENDDO
     3137             ENDDO
     3138          ENDIF
     3139
     3140       CASE ( 'rad_sw_in' )
     3141          IF ( av == 0 )  THEN
     3142             DO  i = 1, mask_size_l(mid,1)
     3143                DO  j = 1, mask_size_l(mid,2)
     3144                   DO  k = 1, mask_size_l(mid,3)
     3145                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
     3146                                            mask_j(mid,j),mask_i(mid,i))
     3147                    ENDDO
     3148                 ENDDO
     3149              ENDDO
     3150          ELSE
     3151             DO  i = 1, mask_size_l(mid,1)
     3152                DO  j = 1, mask_size_l(mid,2)
     3153                   DO  k = 1, mask_size_l(mid,3)
     3154                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
     3155                                               mask_j(mid,j),mask_i(mid,i))
     3156                   ENDDO
     3157                ENDDO
     3158             ENDDO
     3159          ENDIF
     3160
     3161       CASE ( 'rad_sw_out' )
     3162          IF ( av == 0 )  THEN
     3163             DO  i = 1, mask_size_l(mid,1)
     3164                DO  j = 1, mask_size_l(mid,2)
     3165                   DO  k = 1, mask_size_l(mid,3)
     3166                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
     3167                                            mask_j(mid,j),mask_i(mid,i))
     3168                    ENDDO
     3169                 ENDDO
     3170              ENDDO
     3171          ELSE
     3172             DO  i = 1, mask_size_l(mid,1)
     3173                DO  j = 1, mask_size_l(mid,2)
     3174                   DO  k = 1, mask_size_l(mid,3)
     3175                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
     3176                                               mask_j(mid,j),mask_i(mid,i))
     3177                   ENDDO
     3178                ENDDO
     3179             ENDDO
     3180          ENDIF
     3181
     3182       CASE ( 'rad_sw_cs_hr' )
     3183          IF ( av == 0 )  THEN
     3184             DO  i = 1, mask_size_l(mid,1)
     3185                DO  j = 1, mask_size_l(mid,2)
     3186                   DO  k = 1, mask_size_l(mid,3)
     3187                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
     3188                                            mask_j(mid,j),mask_i(mid,i))
     3189                    ENDDO
     3190                 ENDDO
     3191              ENDDO
     3192          ELSE
     3193             DO  i = 1, mask_size_l(mid,1)
     3194                DO  j = 1, mask_size_l(mid,2)
     3195                   DO  k = 1, mask_size_l(mid,3)
     3196                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
     3197                                               mask_j(mid,j),mask_i(mid,i))
     3198                   ENDDO
     3199                ENDDO
     3200             ENDDO
     3201          ENDIF
     3202
     3203       CASE ( 'rad_sw_hr' )
     3204          IF ( av == 0 )  THEN
     3205             DO  i = 1, mask_size_l(mid,1)
     3206                DO  j = 1, mask_size_l(mid,2)
     3207                   DO  k = 1, mask_size_l(mid,3)
     3208                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
     3209                                            mask_j(mid,j),mask_i(mid,i))
     3210                    ENDDO
     3211                 ENDDO
     3212              ENDDO
     3213          ELSE
     3214             DO  i = 1, mask_size_l(mid,1)
     3215                DO  j = 1, mask_size_l(mid,2)
     3216                   DO  k = 1, mask_size_l(mid,3)
     3217                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
     3218                                               mask_j(mid,j),mask_i(mid,i))
     3219                   ENDDO
     3220                ENDDO
     3221             ENDDO
     3222          ENDIF
     3223
     3224       CASE DEFAULT
     3225          found = .FALSE.
     3226
     3227    END SELECT
     3228
     3229
     3230 END SUBROUTINE radiation_data_output_mask
     3231
     3232
     3233!------------------------------------------------------------------------------!
     3234!
     3235! Description:
     3236! ------------
     3237!> Subroutine defines masked output variables
     3238!------------------------------------------------------------------------------!
     3239 SUBROUTINE radiation_last_actions
     3240 
     3241
     3242    USE control_parameters
     3243       
     3244    USE kinds
     3245
     3246    IMPLICIT NONE
     3247
     3248    IF ( write_binary(1:4) == 'true' )  THEN
     3249       IF ( ALLOCATED( rad_net ) )  THEN
     3250          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
     3251       ENDIF
     3252       IF ( ALLOCATED( rad_net_av ) )  THEN
     3253          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
     3254       ENDIF 
     3255       IF ( ALLOCATED( rad_lw_in ) )  THEN
     3256          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
     3257       ENDIF
     3258       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
     3259          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
     3260       ENDIF
     3261       IF ( ALLOCATED( rad_lw_out ) )  THEN
     3262          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out
     3263       ENDIF
     3264       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
     3265          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
     3266       ENDIF
     3267       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
     3268          WRITE ( 14 )  'rad_lw_out_change_0 '
     3269          WRITE ( 14 )  rad_lw_out_change_0
     3270       ENDIF
     3271       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
     3272          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
     3273       ENDIF
     3274       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
     3275          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
     3276       ENDIF
     3277       IF ( ALLOCATED( rad_lw_hr ) )  THEN
     3278          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
     3279       ENDIF
     3280       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
     3281          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
     3282       ENDIF
     3283       IF ( ALLOCATED( rad_sw_in ) )  THEN
     3284          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
     3285       ENDIF
     3286       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
     3287          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
     3288       ENDIF
     3289       IF ( ALLOCATED( rad_sw_out ) )  THEN
     3290          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
     3291       ENDIF
     3292       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
     3293          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
     3294       ENDIF
     3295       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
     3296          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
     3297       ENDIF
     3298       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
     3299          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
     3300       ENDIF
     3301       IF ( ALLOCATED( rad_sw_hr ) )  THEN
     3302          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
     3303       ENDIF
     3304       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
     3305          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
     3306       ENDIF
     3307
     3308       WRITE ( 14 )  '*** end rad ***     '
     3309
     3310    ENDIF
     3311
     3312 END SUBROUTINE radiation_last_actions
     3313
     3314
     3315SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
     3316                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
     3317                                     offset_xa, offset_ya, overlap_count,      &
     3318                                     tmp_2d, tmp_3d )
     3319 
     3320
     3321    USE control_parameters
     3322       
     3323    USE indices
     3324   
     3325    USE kinds
     3326   
     3327    USE pegrid
     3328
     3329    IMPLICIT NONE
     3330
     3331    CHARACTER (LEN=20) :: field_char   !<
     3332
     3333    INTEGER(iwp) ::  i               !<
     3334    INTEGER(iwp) ::  k               !<
     3335    INTEGER(iwp) ::  nxlc            !<
     3336    INTEGER(iwp) ::  nxlf            !<
     3337    INTEGER(iwp) ::  nxl_on_file     !<
     3338    INTEGER(iwp) ::  nxrc            !<
     3339    INTEGER(iwp) ::  nxrf            !<
     3340    INTEGER(iwp) ::  nxr_on_file     !<
     3341    INTEGER(iwp) ::  nync            !<
     3342    INTEGER(iwp) ::  nynf            !<
     3343    INTEGER(iwp) ::  nyn_on_file     !<
     3344    INTEGER(iwp) ::  nysc            !<
     3345    INTEGER(iwp) ::  nysf            !<
     3346    INTEGER(iwp) ::  nys_on_file     !<
     3347    INTEGER(iwp) ::  overlap_count   !<
     3348
     3349    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
     3350    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
     3351    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
     3352    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
     3353    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
     3354    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
     3355
     3356    REAL(wp),                                                                  &
     3357       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
     3358          tmp_2d   !<
     3359
     3360    REAL(wp),                                                                  &
     3361       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
     3362          tmp_3d   !<
     3363
     3364
     3365
     3366   IF ( initializing_actions == 'read_restart_data' )  THEN
     3367      READ ( 13 )  field_char
     3368
     3369      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
     3370
     3371         DO  k = 1, overlap_count
     3372
     3373            nxlf = nxlfa(i,k)
     3374            nxlc = nxlfa(i,k) + offset_xa(i,k)
     3375            nxrf = nxrfa(i,k)
     3376            nxrc = nxrfa(i,k) + offset_xa(i,k)
     3377            nysf = nysfa(i,k)
     3378            nysc = nysfa(i,k) + offset_ya(i,k)
     3379            nynf = nynfa(i,k)
     3380            nync = nynfa(i,k) + offset_ya(i,k)
     3381
     3382
     3383            SELECT CASE ( TRIM( field_char ) )
     3384
     3385                CASE ( 'rad_net' )
     3386                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
     3387                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
     3388                   ENDIF 
     3389                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     3390                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     3391                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3392
     3393                CASE ( 'rad_net_av' )
     3394                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
     3395                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
     3396                   ENDIF 
     3397                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     3398                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     3399                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3400                CASE ( 'rad_lw_in' )
     3401                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
     3402                      ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3403                   ENDIF 
     3404                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3405                   rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3406                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3407
     3408                CASE ( 'rad_lw_in_av' )
     3409                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
     3410                      ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3411                   ENDIF 
     3412                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3413                   rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3414                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3415
     3416                CASE ( 'rad_lw_out' )
     3417                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
     3418                      ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3419                   ENDIF 
     3420                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3421                   rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3422                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3423
     3424                CASE ( 'rad_lw_out_av' )
     3425                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
     3426                      ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3427                   ENDIF 
     3428                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3429                   rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3430                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3431
     3432                CASE ( 'rad_lw_out_change_0' )
     3433                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
     3434                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
     3435                   ENDIF
     3436                   IF ( k == 1 )  READ ( 13 )  tmp_2d
     3437                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
     3438                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3439
     3440                CASE ( 'rad_lw_cs_hr' )
     3441                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
     3442                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3443                   ENDIF
     3444                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3445                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3446                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3447
     3448                CASE ( 'rad_lw_cs_hr_av' )
     3449                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
     3450                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3451                   ENDIF
     3452                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3453                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3454                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3455
     3456                CASE ( 'rad_lw_hr' )
     3457                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
     3458                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3459                   ENDIF
     3460                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3461                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3462                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3463
     3464                CASE ( 'rad_lw_hr_av' )
     3465                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
     3466                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3467                   ENDIF
     3468                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3469                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3470                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3471
     3472                CASE ( 'rad_sw_in' )
     3473                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
     3474                      ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3475                   ENDIF 
     3476                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3477                   rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3478                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3479
     3480                CASE ( 'rad_sw_in_av' )
     3481                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
     3482                      ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3483                   ENDIF 
     3484                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3485                   rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3486                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3487
     3488                CASE ( 'rad_sw_out' )
     3489                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
     3490                      ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3491                   ENDIF 
     3492                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3493                   rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3494                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3495
     3496                CASE ( 'rad_sw_out_av' )
     3497                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
     3498                      ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3499                   ENDIF 
     3500                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3501                   rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3502                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3503
     3504                CASE ( 'rad_sw_cs_hr' )
     3505                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
     3506                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3507                   ENDIF
     3508                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3509                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3510                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3511
     3512                CASE ( 'rad_sw_cs_hr_av' )
     3513                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
     3514                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3515                   ENDIF
     3516                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3517                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3518                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3519
     3520                CASE ( 'rad_sw_hr' )
     3521                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
     3522                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3523                   ENDIF
     3524                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3525                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3526                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3527
     3528                CASE ( 'rad_sw_hr_av' )
     3529                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
     3530                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     3531                   ENDIF
     3532                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     3533                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3534                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3535
     3536               CASE DEFAULT
     3537                  WRITE( message_string, * ) 'unknown variable named "',       &
     3538                                        TRIM( field_char ), '" found in',      &
     3539                                        '&data from prior run on PE ', myid
     3540                  CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
     3541                                 0 )
     3542
     3543            END SELECT
     3544
     3545         ENDDO
     3546
     3547         READ ( 13 )  field_char
     3548
     3549      ENDDO
     3550   ENDIF
     3551
     3552 END SUBROUTINE radiation_read_restart_data
     3553
    22083554
    22093555 END MODULE radiation_model_mod

  • palm/trunk/SOURCE/read_3d_binary.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Bugfix: read of land surface data only when module is switched on.
     22! Radiation parts are now done in the respective module.
     23! Binary version increased to 4.5.
    2224!
    2325! Former revisions:
    … …  
    138140
    139141    USE radiation_model_mod,                                                   &
    140         ONLY: rad_net, rad_net_av, rad_lw_in, rad_lw_in_av, rad_lw_out,        &
    141               rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av,          &
    142               rad_lw_out_av, rad_lw_out_change_0, rad_sw_in, rad_sw_in_av,     &
    143               rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av,        &
    144               rad_sw_hr, rad_sw_hr_av
     142        ONLY: radiation, radiation_read_restart_data
    145143
    146144    USE random_function_mod,                                                   &
    … …  
    328326!--    First compare the version numbers
    329327       READ ( 13 )  version_on_file
    330        binary_version = '4.4'
     328       binary_version = '4.5'
    331329       IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
    332330          WRITE( message_string, * ) 'version mismatch concerning data ',      &
    … …  
    698696                   qv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    699697                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    700 
    701                 CASE ( 'rad_net' )
    702                    IF ( .NOT. ALLOCATED( rad_net ) )  THEN
    703                       ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
    704                    ENDIF 
    705                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    706                    rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    707                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    708 
    709                 CASE ( 'rad_net_av' )
    710                    IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
    711                       ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
    712                    ENDIF 
    713                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    714                    rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    715                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    716                 CASE ( 'rad_lw_in' )
    717                    IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
    718                       ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    719                    ENDIF 
    720                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    721                    rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    722                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    723 
    724                 CASE ( 'rad_lw_in_av' )
    725                    IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
    726                       ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    727                    ENDIF 
    728                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    729                    rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    730                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    731 
    732                 CASE ( 'rad_lw_out' )
    733                    IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
    734                       ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    735                    ENDIF 
    736                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    737                    rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    738                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    739 
    740                 CASE ( 'rad_lw_out_av' )
    741                    IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
    742                       ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    743                    ENDIF 
    744                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    745                    rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    746                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    747 
    748                 CASE ( 'rad_lw_out_change_0' )
    749                    IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
    750                       ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
    751                    ENDIF
    752                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    753                    rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
    754                               = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    755 
    756                 CASE ( 'rad_lw_cs_hr' )
    757                    IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
    758                       ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    759                    ENDIF
    760                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    761                    rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    762                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    763 
    764                 CASE ( 'rad_lw_cs_hr_av' )
    765                    IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
    766                       ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    767                    ENDIF
    768                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    769                    rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    770                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    771 
    772                 CASE ( 'rad_lw_hr' )
    773                    IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
    774                       ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    775                    ENDIF
    776                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    777                    rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    778                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    779 
    780                 CASE ( 'rad_lw_hr_av' )
    781                    IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
    782                       ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    783                    ENDIF
    784                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    785                    rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    786                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    787 
    788                 CASE ( 'rad_sw_in' )
    789                    IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
    790                       ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    791                    ENDIF 
    792                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    793                    rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    794                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    795 
    796                 CASE ( 'rad_sw_in_av' )
    797                    IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    798                       ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    799                    ENDIF 
    800                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    801                    rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    802                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    803 
    804                 CASE ( 'rad_sw_out' )
    805                    IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
    806                       ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    807                    ENDIF 
    808                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    809                    rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    810                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    811 
    812                 CASE ( 'rad_sw_out_av' )
    813                    IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
    814                       ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    815                    ENDIF 
    816                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    817                    rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    818                              tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    819 
    820                 CASE ( 'rad_sw_cs_hr' )
    821                    IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
    822                       ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    823                    ENDIF
    824                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    825                    rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    826                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    827 
    828                 CASE ( 'rad_sw_cs_hr_av' )
    829                    IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
    830                       ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    831                    ENDIF
    832                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    833                    rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    834                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    835 
    836                 CASE ( 'rad_sw_hr' )
    837                    IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
    838                       ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    839                    ENDIF
    840                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    841                    rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    842                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    843 
    844                 CASE ( 'rad_sw_hr_av' )
    845                    IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
    846                       ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    847                    ENDIF
    848                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    849                    rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    850                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    851698
    852699                CASE ( 'random_iv' )  ! still unresolved issue
    … …  
    12741121!
    12751122!--    Read land surface restart data
    1276        CALL lsm_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,  &
    1277                                     nynfa, nyn_on_file, nysfa, nys_on_file,    &
    1278                                     offset_xa, offset_ya, overlap_count(i),    &
    1279                                     tmp_2d )
     1123       IF ( land_surface )  THEN
     1124          CALL lsm_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,            &
     1125                                      nxr_on_file, nynfa, nyn_on_file, nysfa,  &
     1126                                      nys_on_file, offset_xa, offset_ya,       &
     1127                                      overlap_count(i), tmp_2d )
     1128       ENDIF
     1129
     1130!
     1131!--    Read radiation restart data
     1132       IF ( radiation )  THEN
     1133          CALL radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,      &
     1134                                            nxr_on_file, nynfa, nyn_on_file,   &
     1135                                            nysfa, nys_on_file, offset_xa,     &
     1136                                            offset_ya, overlap_count(i),       &
     1137                                            tmp_2d, tmp_3d )
     1138       ENDIF
    12801139
    12811140!

  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Radiation actions are now done directly in the respective module
    2222!
    2323! Former revisions:
    … …  
    138138
    139139    USE radiation_model_mod,                                                   &
    140         ONLY:  rad_net, rad_net_av, rad_sw_in, rad_sw_in_av, rad_sw_out,       &
    141                rad_sw_out_av, rad_sw_cs_hr, rad_sw_cs_hr_av, rad_sw_hr,        &
    142                rad_sw_hr_av, rad_lw_in, rad_lw_in_av, rad_lw_out,              &
    143                rad_lw_out_av, rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr,        &
    144                rad_lw_hr_av
     140        ONLY:  radiation, radiation_3d_data_averaging
    145141
    146142
    … …  
    291287                qv_av = 0.0_wp
    292288
    293              CASE ( 'rad_net*' )
    294                 IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
    295                    ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
    296                 ENDIF
    297                 rad_net_av = 0.0_wp
    298 
    299              CASE ( 'rad_lw_in' )
    300                 IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
    301                    ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    302                 ENDIF
    303                 rad_lw_in_av = 0.0_wp
    304 
    305              CASE ( 'rad_lw_out' )
    306                 IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
    307                    ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    308                 ENDIF
    309                 rad_lw_out_av = 0.0_wp
    310 
    311              CASE ( 'rad_lw_cs_hr' )
    312                 IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
    313                    ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
    314                 ENDIF
    315                 rad_lw_cs_hr_av = 0.0_wp
    316 
    317              CASE ( 'rad_lw_hr' )
    318                 IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
    319                    ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
    320                 ENDIF
    321                 rad_lw_hr_av = 0.0_wp
    322 
    323              CASE ( 'rad_sw_in' )
    324                 IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    325                    ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    326                 ENDIF
    327                 rad_sw_in_av = 0.0_wp
    328 
    329              CASE ( 'rad_sw_out' )
    330                 IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
    331                    ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    332                 ENDIF
    333                 rad_sw_out_av = 0.0_wp
    334 
    335              CASE ( 'rad_sw_cs_hr' )
    336                 IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
    337                    ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
    338                 ENDIF
    339                 rad_sw_cs_hr_av = 0.0_wp
    340 
    341              CASE ( 'rad_sw_hr' )
    342                 IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
    343                    ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
    344                 ENDIF
    345                 rad_sw_hr_av = 0.0_wp
    346 
    347289             CASE ( 'rho' )
    348290                IF ( .NOT. ALLOCATED( rho_av ) )  THEN
    … …  
    429371                IF ( land_surface )  THEN
    430372                   CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
     373                ENDIF
     374
     375!
     376!--             Radiation quantity
     377                IF ( radiation )  THEN
     378                   CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
    431379                ENDIF
    432380
    … …  
    655603             ENDDO
    656604
    657           CASE ( 'rad_net*' )
    658              DO  i = nxlg, nxrg
    659                 DO  j = nysg, nyng
    660                    rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
    661                 ENDDO
    662              ENDDO
    663 
    664           CASE ( 'rad_lw_in' )
    665              DO  i = nxlg, nxrg
    666                 DO  j = nysg, nyng
    667                    DO  k = nzb, nzt+1
    668                       rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
    669                    ENDDO
    670                 ENDDO
    671              ENDDO
    672 
    673           CASE ( 'rad_lw_out' )
    674              DO  i = nxlg, nxrg
    675                 DO  j = nysg, nyng
    676                    DO  k = nzb, nzt+1
    677                       rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
    678                    ENDDO
    679                 ENDDO
    680              ENDDO
    681 
    682           CASE ( 'rad_lw_cs_hr' )
    683              DO  i = nxlg, nxrg
    684                 DO  j = nysg, nyng
    685                    DO  k = nzb, nzt+1
    686                       rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
    687                    ENDDO
    688                 ENDDO
    689              ENDDO
    690 
    691           CASE ( 'rad_lw_hr' )
    692              DO  i = nxlg, nxrg
    693                 DO  j = nysg, nyng
    694                    DO  k = nzb, nzt+1
    695                       rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
    696                    ENDDO
    697                 ENDDO
    698              ENDDO
    699 
    700           CASE ( 'rad_sw_in' )
    701              DO  i = nxlg, nxrg
    702                 DO  j = nysg, nyng
    703                    DO  k = nzb, nzt+1
    704                       rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
    705                    ENDDO
    706                 ENDDO
    707              ENDDO
    708 
    709           CASE ( 'rad_sw_out' )
    710              DO  i = nxlg, nxrg
    711                 DO  j = nysg, nyng
    712                    DO  k = nzb, nzt+1
    713                       rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
    714                    ENDDO
    715                 ENDDO
    716              ENDDO
    717 
    718           CASE ( 'rad_sw_cs_hr' )
    719              DO  i = nxlg, nxrg
    720                 DO  j = nysg, nyng
    721                    DO  k = nzb, nzt+1
    722                       rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
    723                    ENDDO
    724                 ENDDO
    725              ENDDO
    726 
    727           CASE ( 'rad_sw_hr' )
    728              DO  i = nxlg, nxrg
    729                 DO  j = nysg, nyng
    730                    DO  k = nzb, nzt+1
    731                       rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
    732                    ENDDO
    733                 ENDDO
    734              ENDDO
    735 
    736605          CASE ( 'rho' )
    737606             DO  i = nxlg, nxrg
    … …  
    854723
    855724!
     725!--          Radiation quantity
     726             IF ( radiation )  THEN
     727                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
     728             ENDIF
     729
     730!
    856731!--          User-defined quantity
    857732             CALL user_3d_data_averaging( 'sum', doav(ii) )

  • palm/trunk/SOURCE/time_integration.f90

    r1961 r1976  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Simplified calls to radiation model
    2222!
    2323! Former revisions:
    … …  
    321321
    322322    USE radiation_model_mod,                                                   &
    323         ONLY: dt_radiation, force_radiation_call, radiation,                   &
    324               radiation_clearsky, radiation_constant, radiation_rrtmg,         &
    325               radiation_scheme, skip_time_do_radiation, time_radiation
     323        ONLY: dt_radiation, force_radiation_call, radiation, radiation_control,&
     324              skip_time_do_radiation, time_radiation
    326325
    327326    USE spectra_mod,                                                           &
    … …  
    900899                ENDIF
    901900
    902                 IF ( radiation_scheme == 'clear-sky' )  THEN
    903                    CALL radiation_clearsky
    904                 ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
    905                    CALL radiation_rrtmg
    906                 ELSE
    907                    CALL radiation_constant
    908                 ENDIF
     901                CALL radiation_control
    909902
    910903                CALL cpu_log( log_point(50), 'radiation', 'stop' )

  • palm/trunk/SOURCE/write_3d_binary.f90

    r1973 r1976  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Radiation actions are now done directly in the respective module.
     22! Binary version increased to 4.5.
    2223!
    2324! Former revisions:
    … …  
    122123    USE pegrid
    123124   
    124     USE radiation_model_mod,                                                   &
    125         ONLY: radiation, rad_net, rad_net_av, rad_lw_in, rad_lw_in_av,         &
    126               rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, rad_lw_cs_hr,   &
    127               rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,             &
    128               rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,           &
    129               rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av
    130 
    131125    USE random_function_mod,                                                   &
    132126        ONLY:  random_iv, random_iy
    … …  
    149143!
    150144!-- Write arrays.
    151     binary_version = '4.4'
     145    binary_version = '4.5'
    152146
    153147    WRITE ( 14 )  binary_version
    … …  
    254248       WRITE ( 14 )  'sswst               ';  WRITE ( 14 ) sswst
    255249    ENDIF   
    256     IF ( ALLOCATED( rad_net ) )  THEN
    257        WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
    258     ENDIF
    259     IF ( radiation )  THEN
    260        IF ( ALLOCATED( rad_net_av ) )  THEN
    261           WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
    262        ENDIF 
    263        IF ( ALLOCATED( rad_lw_in ) )  THEN
    264           WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
    265        ENDIF
    266        IF ( ALLOCATED( rad_lw_in_av ) )  THEN
    267           WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
    268        ENDIF
    269        IF ( ALLOCATED( rad_lw_out ) )  THEN
    270           WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out
    271        ENDIF
    272        IF ( ALLOCATED( rad_lw_out_av ) )  THEN
    273           WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
    274        ENDIF
    275        IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
    276           WRITE ( 14 )  'rad_lw_out_change_0 '
    277           WRITE ( 14 )  rad_lw_out_change_0
    278        ENDIF
    279        IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
    280           WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
    281        ENDIF
    282        IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
    283           WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
    284        ENDIF
    285        IF ( ALLOCATED( rad_lw_hr ) )  THEN
    286           WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
    287        ENDIF
    288        IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
    289           WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
    290        ENDIF
    291        IF ( ALLOCATED( rad_sw_in ) )  THEN
    292           WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
    293        ENDIF
    294        IF ( ALLOCATED( rad_sw_in_av ) )  THEN
    295           WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
    296        ENDIF
    297        IF ( ALLOCATED( rad_sw_out ) )  THEN
    298           WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
    299        ENDIF
    300        IF ( ALLOCATED( rad_sw_out_av ) )  THEN
    301           WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
    302        ENDIF
    303        IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
    304           WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
    305        ENDIF
    306        IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
    307           WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
    308        ENDIF
    309        IF ( ALLOCATED( rad_sw_hr ) )  THEN
    310           WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
    311        ENDIF
    312        IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
    313           WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
    314        ENDIF
    315     ENDIF
    316250    IF ( ocean )  THEN
    317251       IF ( ALLOCATED( rho_av ) )  THEN
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |