Changeset 2299


Ignore:
Timestamp:
Jun 29, 2017 10:14:38 AM (4 years ago)
Author:
maronga
Message:

improvements for spinup mechanism

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

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

    r2298 r2299  
    2525! -----------------
    2626! $Id$
     27! Modified output for spinups
     28!
     29! 2298 2017-06-29 09:28:18Z raasch
    2730! MPI2 related parts removed
    2831!
     
    507510    ENDIF
    508511#if defined( __parallel )
    509     IF ( coupling_start_time /= 0.0_wp )  THEN
     512    IF ( coupling_start_time /= 0.0_wp  .AND. .NOT. spinup )  THEN
    510513       IF ( coupling_start_time > simulated_time_at_begin )  THEN
    511514          WRITE ( io, 109 )
  • palm/trunk/SOURCE/init_1d_model.f90

    r2101 r2299  
    2121! -----------------
    2222!
    23 !
     23! 
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
     27! Removed german text
     28!
     29! 2101 2017-01-05 16:42:31Z suehring
    2730!
    2831! 2059 2016-11-10 14:20:40Z maronga
     
    851854!
    852855!-- formats
    853 100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ &
     856100 FORMAT (///'1D run control output:'/ &
    854857              &'------------------------------'// &
    855858           &'ITER.  HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2298 r2299  
    2525! -----------------
    2626! $Id$
     27! Removed pt_p from USE statement. Adjusted call to lsm_soil_model to allow
     28! spinups without soil moisture prediction
     29!
     30! 2298 2017-06-29 09:28:18Z raasch
    2731! type of write_binary changed from CHARACTER to LOGICAL
    2832!
     
    228232 
    229233    USE arrays_3d,                                                             &
    230         ONLY:  hyp, pt, pt_p, prr, q, q_p, ql, vpt, u, v, w
     234        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
    231235
    232236    USE cloud_parameters,                                                      &
     
    28842888!> temperature and water content.
    28852889!------------------------------------------------------------------------------!
    2886     SUBROUTINE lsm_soil_model( horizontal, l )
     2890    SUBROUTINE lsm_soil_model( horizontal, l, calc_soil_moisture )
    28872891
    28882892
     
    28922896       INTEGER(iwp) ::  l       !< surface-data type index indication facing
    28932897       INTEGER(iwp) ::  m       !< running index
     2898
     2899       LOGICAL, INTENT(IN) ::  calc_soil_moisture !< flag indicating whether soil moisture shall be calculated or not.
    28942900
    28952901       LOGICAL      ::  horizontal !< flag indication horizontal wall, required to set pointer accordingly
     
    30593065
    30603066             ENDDO
     3067
     3068          ENDIF
     3069
     3070       ENDDO
     3071
     3072
     3073       DO  m = 1, surf%ns
     3074
     3075          IF (  .NOT.  surf%water_surface(m)  .AND.  calc_soil_moisture )  THEN
     3076
    30613077
    30623078!
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2298 r2299  
    2525! -----------------
    2626! $Id$
     27! Improved syntax layout
     28!
     29! 2298 2017-06-29 09:28:18Z raasch
    2730! type of write_binary changed from CHARACTER to LOGICAL
    2831!
     
    505508!
    506509!-- Public variables and constants / NEEDS SORTING
    507     PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
    508            rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
     510    PUBLIC decl_1, decl_2, decl_3, dots_rad, dt_radiation, force_radiation_call,&
     511           lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
    509512           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
    510513           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
     
    632635!> Check data output of profiles for radiation model
    633636!------------------------------------------------------------------------------! 
    634     SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
     637    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
     638               dopr_unit )
    635639 
    636640       USE arrays_3d,                                                          &
     
    12371241          DO j = nysg, nyng
    12381242!
    1239 !--          Obtain vertical index of topography top. So far it is identical to nzb.
     1243!--          Obtain vertical index of topography top. So far it is identical to
     1244!--          nzb.
    12401245             k = MAXLOC(                                                       &
    12411246                        MERGE( 1, 0,                                           &
     
    13141319
    13151320 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
    1316  2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
     1321 2 FORMAT ('    Day of the year at model start :   day_init = ',I3             &
    13171322            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
    13181323 3 FORMAT (//' Radiation model information:'/                                  &
    13191324              ' ----------------------------'/)
    1320  4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
     1325 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
    13211326           // 'W/m**2')
    1322  5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
     1327 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
    13231328                   ' default)')
    13241329 6 FORMAT ('    --> RRTMG scheme is used')
     
    14931498!--                   Calculate cloud droplet effective radius
    14941499                      IF ( cloud_physics )  THEN
    1495                          rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
    1496                                            * rho_surface                          &
    1497                                            / ( 4.0_wp * pi * nc_const * rho_l )   &
    1498                                            )**0.33333333333333_wp                 &
     1500                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
     1501                                           * rho_surface                       &
     1502                                           / ( 4.0_wp * pi * nc_const * rho_l )&
     1503                                           )**0.33333333333333_wp              &
    14991504                                           * EXP( LOG( sigma_gc )**2 )
    15001505
     
    16361641!--    Calculate current day and time based on the initial values and simulation
    16371642!--    time
    1638        day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
     1643       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)&
    16391644                               / 86400.0_wp ), KIND=iwp)
    16401645       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
     
    16481653!
    16491654!--    Calculate cosine of solar zenith angle
    1650        zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
     1655       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
    16511656                                            * COS(hour_angle)
    16521657       zenith(0) = MAX(0.0_wp,zenith(0))
     
    16551660!--    Calculate solar directional vector
    16561661       IF ( sun_direction )  THEN
     1662
     1663!
    16571664!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
    16581665          sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
     1666
     1667!
    16591668!--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
    16601669          sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
     
    19421951          rrtm_h2ovmr(0,k) = q_snd(k)
    19431952       ENDDO
    1944        rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
     1953       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
    19451954                                - rrtm_tlay(0,nzt_rad-1)
    19461955       DO k = nzt+9, nzt_rad+1
     
    20582067
    20592068
    2060        REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
    2061                                                  rrtm_play_tmp, & !< temporary array for pressure zu-levels
    2062                                                  rrtm_plev_tmp, & !< temporary array for pressure zw-levels
    2063                                                  trace_path_tmp   !< temporary array for storing trace gas path data
     2069       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,          & !< pressure levels for the absorbers
     2070                                                 rrtm_play_tmp,  & !< temporary array for pressure zu-levels
     2071                                                 rrtm_plev_tmp,  & !< temporary array for pressure zw-levels
     2072                                                 trace_path_tmp    !< temporary array for storing trace gas path data
    20642073
    20652074       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
     
    24972506                DO  j = nysg, nyng
    24982507                   DO  k = nzb, nzt+1
    2499                       rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
     2508                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
     2509                                             + rad_lw_out(k,j,i)
    25002510                   ENDDO
    25012511                ENDDO
     
    25062516                DO  j = nysg, nyng
    25072517                   DO  k = nzb, nzt+1
    2508                       rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
     2518                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
     2519                                               + rad_lw_cs_hr(k,j,i)
    25092520                   ENDDO
    25102521                ENDDO
     
    25152526                DO  j = nysg, nyng
    25162527                   DO  k = nzb, nzt+1
    2517                       rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
     2528                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
     2529                                            + rad_lw_hr(k,j,i)
    25182530                   ENDDO
    25192531                ENDDO
     
    25242536                DO  j = nysg, nyng
    25252537                   DO  k = nzb, nzt+1
    2526                       rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
     2538                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
     2539                                            + rad_sw_in(k,j,i)
    25272540                   ENDDO
    25282541                ENDDO
     
    25332546                DO  j = nysg, nyng
    25342547                   DO  k = nzb, nzt+1
    2535                       rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
     2548                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
     2549                                             + rad_sw_out(k,j,i)
    25362550                   ENDDO
    25372551                ENDDO
     
    25422556                DO  j = nysg, nyng
    25432557                   DO  k = nzb, nzt+1
    2544                       rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
     2558                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
     2559                                               + rad_sw_cs_hr(k,j,i)
    25452560                   ENDDO
    25462561                ENDDO
     
    25512566                DO  j = nysg, nyng
    25522567                   DO  k = nzb, nzt+1
    2553                       rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
     2568                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
     2569                                            + rad_sw_hr(k,j,i)
    25542570                   ENDDO
    25552571                ENDDO
     
    25682584             DO  i = nxlg, nxrg
    25692585                DO  j = nysg, nyng
    2570                    rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
     2586                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, &
     2587                                     KIND=wp )
    25712588                ENDDO
    25722589             ENDDO
     
    25762593                DO  j = nysg, nyng
    25772594                   DO  k = nzb, nzt+1
    2578                       rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2595                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i)                &
     2596                                            / REAL( average_count_3d, KIND=wp )
    25792597                   ENDDO
    25802598                ENDDO
     
    25852603                DO  j = nysg, nyng
    25862604                   DO  k = nzb, nzt+1
    2587                       rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2605                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
     2606                                             / REAL( average_count_3d, KIND=wp )
    25882607                   ENDDO
    25892608                ENDDO
     
    25942613                DO  j = nysg, nyng
    25952614                   DO  k = nzb, nzt+1
    2596                       rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2615                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
     2616                                             / REAL( average_count_3d, KIND=wp )
    25972617                   ENDDO
    25982618                ENDDO
     
    26032623                DO  j = nysg, nyng
    26042624                   DO  k = nzb, nzt+1
    2605                       rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2625                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
     2626                                            / REAL( average_count_3d, KIND=wp )
    26062627                   ENDDO
    26072628                ENDDO
     
    26122633                DO  j = nysg, nyng
    26132634                   DO  k = nzb, nzt+1
    2614                       rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2635                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
     2636                                            / REAL( average_count_3d, KIND=wp )
    26152637                   ENDDO
    26162638                ENDDO
     
    26212643                DO  j = nysg, nyng
    26222644                   DO  k = nzb, nzt+1
    2623                       rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2645                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
     2646                                             / REAL( average_count_3d, KIND=wp )
    26242647                   ENDDO
    26252648                ENDDO
     
    26302653                DO  j = nysg, nyng
    26312654                   DO  k = nzb, nzt+1
    2632                       rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2655                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
     2656                                             / REAL( average_count_3d, KIND=wp )
    26332657                   ENDDO
    26342658                ENDDO
     
    26392663                DO  j = nysg, nyng
    26402664                   DO  k = nzb, nzt+1
    2641                       rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     2665                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
     2666                                            / REAL( average_count_3d, KIND=wp )
    26422667                   ENDDO
    26432668                ENDDO
     
    34123437
    34133438
    3414 SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
    3415                                      nynfa, nyn_on_file, nysfa, nys_on_file,   &
    3416                                      offset_xa, offset_ya, overlap_count,      &
    3417                                      tmp_2d, tmp_3d )
     3439SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,          &
     3440                                        nxr_on_file, nynfa, nyn_on_file, nysfa,&
     3441                                        nys_on_file, offset_xa, offset_ya,     &
     3442                                        overlap_count, tmp_2d, tmp_3d )
    34183443 
    34193444
     
    34913516                   ENDIF 
    34923517                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3493                    rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    3494                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3518                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =         &
     3519                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34953520
    34963521                CASE ( 'rad_net_av' )
     
    34993524                   ENDIF 
    35003525                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3501                    rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    3502                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3526                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
     3527                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35033528                CASE ( 'rad_lw_in' )
    35043529                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
     
    35153540                         READ ( 13 )  tmp_3d2
    35163541                         rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
    3517                              tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3542                            tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    35183543                      ELSE
    35193544                         READ ( 13 )  tmp_3d
     
    36023627                   ENDIF
    36033628                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3604                    rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3629                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    36053630                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    36063631
     
    36183643                   ENDIF
    36193644                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3620                    rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3645                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
    36213646                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    36223647
     
    36263651                   ENDIF
    36273652                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3628                    rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3653                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    36293654                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    36303655
     
    37223747                   ENDIF
    37233748                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3724                    rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3749                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    37253750                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    37263751
     
    37383763                   ENDIF
    37393764                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3740                    rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3765                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
    37413766                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    37423767
     
    37463771                   ENDIF
    37473772                   IF ( k == 1 )  READ ( 13 )  tmp_3d
    3748                    rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3773                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    37493774                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    37503775
     
    37533778                                        TRIM( field_char ), '" found in',      &
    37543779                                        '&data from prior run on PE ', myid
    3755                   CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, 0, 6, &
    3756                                 0 )
     3780                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, &
     3781                                0, 6, 0 )
    37573782
    37583783            END SELECT
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r2292 r2299  
    2525! -----------------
    2626! $Id$
     27! Adjusted for allow separate spinups of LSM and atmosphere code
     28!
     29! 2292 2017-06-20 09:51:42Z schwenkel
    2730! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
    2831! includes two more prognostic equations for cloud drop concentration (nc) 
     
    203206               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
    204207               run_coupled, surface_pressure, simulated_time, terminate_run,   &
    205                urban_surface, zeta_max, zeta_min
     208               time_since_reference_point, urban_surface, zeta_max, zeta_min
    206209
    207210    USE grid_variables,                                                        &
     
    19211924!
    19221925!--       Compute the vertical kinematic heat flux
    1923           IF (  .NOT.  constant_heatflux  .AND.  ( simulated_time <=           &
    1924                skip_time_do_lsm  .OR.  .NOT.  land_surface )  .AND.            &
    1925                .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
     1926          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
     1927               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
     1928               .NOT.  land_surface )  .AND.  .NOT.  urban_surface  .AND.       &
     1929               .NOT. downward )  THEN
    19261930             !$OMP PARALLEL DO PRIVATE( i, j, k )
    19271931             DO  m = 1, surf%ns
     
    19351939!--       Compute the vertical water flux
    19361940          IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
    1937                 ( simulated_time <= skip_time_do_lsm                           &
    1938                .OR.  .NOT.  land_surface )  .AND.  .NOT. downward )  THEN
     1941               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
     1942               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface )  .AND.    &
     1943               .NOT. downward )  THEN
    19391944             !$OMP PARALLEL DO PRIVATE( i, j, k )
    19401945             DO  m = 1, surf%ns
  • palm/trunk/SOURCE/time_integration.f90

    r2292 r2299  
    2525! -----------------
    2626! $Id$
     27! Call of soil model adjusted
     28!
     29! 2292 2017-06-20 09:51:42Z schwenkel
    2730! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
    2831! includes two more prognostic equations for cloud drop concentration (nc) 
     
    770773!--          If required, solve the energy balance for the surface and run soil
    771774!--          model. Call for horizontal as well as vertical surfaces
    772              IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
     775             IF ( land_surface .AND. time_since_reference_point > skip_time_do_lsm)  THEN
    773776
    774777                CALL cpu_log( log_point(54), 'land_surface', 'start' )
     
    776779!--             Call for horizontal upward-facing surfaces
    777780                CALL lsm_energy_balance( .TRUE., -1 )
    778                 CALL lsm_soil_model( .TRUE., -1 )
     781                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
    779782!
    780783!--             Call for northward-facing surfaces
    781784                CALL lsm_energy_balance( .FALSE., 0 )
    782                 CALL lsm_soil_model( .FALSE., 0 )
     785                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
    783786!
    784787!--             Call for southward-facing surfaces
    785788                CALL lsm_energy_balance( .FALSE., 1 )
    786                 CALL lsm_soil_model( .FALSE., 1 )
     789                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
    787790!
    788791!--             Call for eastward-facing surfaces
    789792                CALL lsm_energy_balance( .FALSE., 2 )
    790                 CALL lsm_soil_model( .FALSE., 2 )
     793                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
    791794!
    792795!--             Call for westward-facing surfaces
    793796                CALL lsm_energy_balance( .FALSE., 3 )
    794                 CALL lsm_soil_model( .FALSE., 3 )
     797                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
    795798
    796799                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
     
    829832!--       If required, calculate radiative fluxes and heating rates
    830833          IF ( radiation .AND. intermediate_timestep_count                     &
    831                == intermediate_timestep_count_max .AND. simulated_time >    &
     834               == intermediate_timestep_count_max .AND. time_since_reference_point >    &
    832835               skip_time_do_radiation )  THEN
    833836
  • palm/trunk/SOURCE/time_integration_spinup.f90

    r2297 r2299  
    2525! -----------------
    2626! $Id$
     27! Call of soil model adjusted to avoid prognostic equation for soil moisture
     28! during spinup.
     29! Better representation of diurnal cycle of near-surface temperature.
     30! Excluded prognostic equation for soil moisture during spinup.
     31! Added output of run control data for spinup.
     32!
     33! 2297 2017-06-28 14:35:57Z scharf
    2734! bugfixes
    2835!
    2936! 2296 2017-06-28 07:53:56Z maronga
    3037! Initial revision
    31 !
    32 !
    3338!
    3439!
     
    5055               dt_spinup, humidity, intermediate_timestep_count,               &
    5156               intermediate_timestep_count_max, land_surface,                  &
    52                nr_timesteps_this_run, simulated_time, simulated_time_chr,      &
     57               simulated_time, simulated_time_chr,      &
    5358               skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
    5459               spinup_time, timestep_count, timestep_scheme, time_dopr,        &
     
    6772
    6873    USE land_surface_model_mod,                                                &
    69         ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
    70                skip_time_do_lsm
    71 
    72     USE pegrid
     74        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
     75
     76    USE pegrid,                                                                &
     77        ONLY:  myid
    7378
    7479    USE kinds
     
    7681    USE radiation_model_mod,                                                   &
    7782        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
    78                radiation_control, skip_time_do_radiation, time_radiation
     83               radiation_control, rad_sw_in, time_radiation, time_utc_init
    7984
    8085    USE statistics,                                                            &
     
    99104    CHARACTER (LEN=9) ::  time_to_string          !<
    100105 
    101     INTEGER(iwp) :: i
    102     INTEGER(iwp) :: j
    103     INTEGER(iwp) :: k
    104     INTEGER(iwp) :: l
    105     INTEGER(iwp) :: m
     106    INTEGER(iwp) ::  i !< running index
     107    INTEGER(iwp) ::  j !< running index
     108    INTEGER(iwp) ::  k !< running index
     109    INTEGER(iwp) ::  l !< running index
     110    INTEGER(iwp) ::  m !< running index
     111
     112    INTEGER(iwp) :: current_timestep_number_spinup = 0  !< number if timestep during spinup
    106113 
     114    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
     115
    107116    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
    108117                 
    109118    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
    110119
    111     REAL(wp) ::  day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
    112 
    113120    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    114121
    115     CALL exchange_horiz( pt_p, nbgp )   
    116     pt_save = pt_p
     122    CALL exchange_horiz( pt, nbgp )   
     123    pt_save = pt
    117124
    118125    CALL location_message( 'starting spinup-sequence', .TRUE. )
     
    137144
    138145
    139 !!!!      Set new values of pt_p here
    140           pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
    141 
     146!
     147!--       Estimate a near-surface air temperature based on the position of the
     148!--       sun and user input about mean temperature and amplitude. The time is
     149!--       shifted by one hour to simulate a lag between air temperature and
     150!--       incoming radiation
     151          pt_spinup = spinup_pt_mean + spinup_pt_amplitude                     &
     152             * solar_angle (time_utc_init + time_since_reference_point - 3600.0)
     153
     154!
     155!--       Map air temperature to all grid points in the vicinity of a surface
     156!--       element
    142157          IF ( land_surface )  THEN
    143158             DO  m = 1, surf_lsm_h%ns
     
    145160                j   = surf_lsm_h%j(m)
    146161                k   = surf_lsm_h%k(m)
    147                 pt_p(k,j,i) = pt_spinup
     162                pt(k,j,i) = pt_spinup
    148163             ENDDO
    149164
     
    153168                   j   = surf_lsm_v(l)%j(m)
    154169                   k   = surf_lsm_v(l)%k(m)
    155                    pt_p(k,j,i) = pt_spinup
     170                   pt(k,j,i) = pt_spinup
    156171                ENDDO
    157172             ENDDO
     
    163178                j   = surf_usm_h%j(m)
    164179                k   = surf_usm_h%k(m)
    165                 pt_p(k,j,i) = pt_spinup
     180                pt(k,j,i) = pt_spinup
    166181             ENDDO
    167182
     
    171186                   j   = surf_usm_v(l)%j(m)
    172187                   k   = surf_usm_v(l)%k(m)
    173                    pt_p(k,j,i) = pt_spinup
     188                   pt(k,j,i) = pt_spinup
    174189                ENDDO
    175190             ENDDO
     
    217232!
    218233!--          If required, solve the energy balance for the surface and run soil
    219 !--          model. Call for horizontal as well as vertical surfaces
    220              IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
     234!--          model. Call for horizontal as well as vertical surfaces.
     235!--          The prognostic equation for soil moisure is switched off
     236             IF ( land_surface )  THEN
    221237
    222238                CALL cpu_log( log_point(54), 'land_surface', 'start' )
     
    224240!--             Call for horizontal upward-facing surfaces
    225241                CALL lsm_energy_balance( .TRUE., -1 )
    226                 CALL lsm_soil_model( .TRUE., -1 )
     242                CALL lsm_soil_model( .TRUE., -1, .FALSE. )
    227243!
    228244!--             Call for northward-facing surfaces
    229245                CALL lsm_energy_balance( .FALSE., 0 )
    230                 CALL lsm_soil_model( .FALSE., 0 )
     246                CALL lsm_soil_model( .FALSE., 0, .FALSE. )
    231247!
    232248!--             Call for southward-facing surfaces
    233249                CALL lsm_energy_balance( .FALSE., 1 )
    234                 CALL lsm_soil_model( .FALSE., 1 )
     250                CALL lsm_soil_model( .FALSE., 1, .FALSE. )
    235251!
    236252!--             Call for eastward-facing surfaces
    237253                CALL lsm_energy_balance( .FALSE., 2 )
    238                 CALL lsm_soil_model( .FALSE., 2 )
     254                CALL lsm_soil_model( .FALSE., 2, .FALSE. )
    239255!
    240256!--             Call for westward-facing surfaces
    241257                CALL lsm_energy_balance( .FALSE., 3 )
    242                 CALL lsm_soil_model( .FALSE., 3 )
     258                CALL lsm_soil_model( .FALSE., 3, .FALSE. )
    243259
    244260                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
     
    262278!--       If required, calculate radiative fluxes and heating rates
    263279          IF ( radiation .AND. intermediate_timestep_count                     &
    264                == intermediate_timestep_count_max .AND. simulated_time >       &
    265                skip_time_do_radiation )  THEN
     280               == intermediate_timestep_count_max )  THEN
    266281
    267282               time_radiation = time_radiation + dt_spinup
     
    292307!
    293308!--    Increase simulation time and output times
    294        nr_timesteps_this_run      = nr_timesteps_this_run + 1
    295        current_timestep_number    = current_timestep_number + 1
     309       current_timestep_number_spinup = current_timestep_number_spinup + 1
    296310       simulated_time             = simulated_time   + dt_spinup
    297311       simulated_time_chr         = time_to_string( simulated_time )
     
    350364!--    Computation and output of run control parameters.
    351365!--    This is also done whenever perturbations have been imposed
    352        IF ( time_run_control >= dt_run_control  .OR.                           &
    353             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
    354        THEN
    355           CALL run_control
    356           IF ( time_run_control >= dt_run_control )  THEN
    357              time_run_control = MOD( time_run_control,                         &
    358                                      MAX( dt_run_control, dt_spinup ) )
    359           ENDIF
     366!        IF ( time_run_control >= dt_run_control  .OR.                           &
     367!             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
     368!        THEN
     369!           CALL run_control
     370!           IF ( time_run_control >= dt_run_control )  THEN
     371!              time_run_control = MOD( time_run_control,                         &
     372!                                      MAX( dt_run_control, dt_spinup ) )
     373!           ENDIF
     374!        ENDIF
     375
     376       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
     377
     378
     379!
     380!--    Run control output
     381       IF ( myid == 0 )  THEN
     382!
     383!--       If necessary, write header
     384          IF ( .NOT. run_control_header_spinup )  THEN
     385             CALL check_open( 15 )
     386             WRITE ( 15, 100 )
     387             run_control_header_spinup = .TRUE.
     388          ENDIF
     389!
     390!--       Write some general information about the spinup in run control file
     391          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_spinup, pt_spinup, rad_sw_in(0,nysg,nxlg)
     392!
     393!--       Write buffer contents to disc immediately
     394          FLUSH( 15 )
    360395       ENDIF
    361396
    362        CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
    363 
    364        IF ( myid == 0 )  THEN
    365           PRINT*, time_since_reference_point
    366        ENDIF
     397
    367398
    368399    ENDDO   ! time loop
     
    375406    DEALLOCATE(pt_save)
    376407
    377     CALL location_message( 'finished time-stepping spinup', .TRUE. )
     408    CALL location_message( 'finished spinup-sequence', .TRUE. )
     409
     410
     411!
     412!-- Formats
     413100 FORMAT (///'Spinup control output:'/  &
     414            '----------------------------------------'// &
     415            'ITER.  HH:MM:SS    DT   PT(z_MO)     SWD'/   &
     416            '----------------------------------------')
     417101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
     418
     419 CONTAINS
     420
     421!
     422!-- Returns the cosine of the solar zenith angle at a given time. This routine
     423!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
     424    FUNCTION solar_angle( local_time )
     425
     426       USE constants,                                                          &
     427           ONLY:  pi
     428       
     429       USE kinds
     430
     431       USE radiation_model_mod,                                                &
     432           ONLY:  day_init, decl_1, decl_2, decl_3, lat, lon
     433
     434       IMPLICIT NONE
     435
     436
     437       REAL(wp) ::  solar_angle     !< cosine of the solar zenith angle
     438
     439       REAL(wp) ::  day             !< day of the year
     440       REAL(wp) ::  declination     !< solar declination angle
     441       REAL(wp) ::  hour_angle      !< solar hour angle
     442       REAL(wp) ::  time_utc        !< current time in UTC
     443       REAL(wp), INTENT(IN) :: local_time
     444!
     445!--    Calculate current day and time based on the initial values and simulation
     446!--    time
     447
     448       day = day_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
     449       time_utc = MOD(local_time, 86400.0_wp)
     450
     451
     452!
     453!--    Calculate solar declination and hour angle   
     454       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
     455       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
     456
     457!
     458!--    Calculate cosine of solar zenith angle
     459       solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
     460                     * COS(hour_angle)
     461
     462
     463    END FUNCTION solar_angle
     464
    378465
    379466 END SUBROUTINE time_integration_spinup
Note: See TracChangeset for help on using the changeset viewer.