Changeset 2723 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jan 5, 2018 9:27:03 AM (7 years ago)
Author:
maronga
Message:

bugfixes for spinup mechanism to work with lsm+usm+radiation

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r2718 r2723  
    2525! -----------------
    2626! $Id$
     27! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
     28!
     29! 2718 2018-01-02 08:49:38Z maronga
    2730! Corrected "Former revisions" section
    2831!
     
    19761979
    19771980
    1978 !        IF ( surf_t_surface_p%var_1d(m) < 270.0_wp .OR. surf_t_surface_p%var_1d(m) > 350.0_wp )  THEN
    1979 !           PRINT*, "myid: ", myid
    1980 !           PRINT*, "m: ", m
    1981 !           PRINT*, "i,j,k: ", i, j, k
    1982 !           PRINT*, "pt_p: ", surf_t_surface_p%var_1d(m)
    1983 !           PRINT*, "f_shf: ", f_shf
    1984 !           PRINT*, "f_qsws: ", f_qsws
    1985 !         ENDIF
    1986 
    1987 
    19881981!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
    19891982       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
     
    45274520       IF ( spinup_time > 0.0_wp )  THEN
    45284521          coupling_start_time = spinup_time
    4529           end_time = end_time + spinup_time
    45304522          IF ( spinup_pt_mean == 9999999.9_wp )  THEN
    45314523             spinup_pt_mean = pt_surface
    45324524          ENDIF
    4533           spinup = .TRUE.
     4525          IF ( .NOT. spinup )  THEN
     4526             end_time = end_time + spinup_time
     4527             spinup = .TRUE.
     4528          ENDIF
    45344529       ENDIF
    45354530
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2718 r2723  
    2525! -----------------
    2626! $Id$
     27! Bugfix in calculation of rad_lw_out (clear-sky). First grid level was used
     28! instead of the surface value
     29!
     30! 2718 2018-01-02 08:49:38Z maronga
    2731! Corrected "Former revisions" section
    2832!
     
    23702374
    23712375             INTEGER(iwp) ::  i         !< index x-direction
    2372              INTEGER(iwp) ::  ioff      !< offset between surface element and adjacent grid point along x
    23732376             INTEGER(iwp) ::  j         !< index y-direction
    2374              INTEGER(iwp) ::  joff      !< offset between surface element and adjacent grid point along y
    23752377             INTEGER(iwp) ::  k         !< index z-direction
    2376              INTEGER(iwp) ::  koff      !< offset between surface element and adjacent grid point along z
    23772378             INTEGER(iwp) ::  m         !< running index for surface elements
    23782379
     
    23852386
    23862387                k = nzut
    2387 !
    2388 !-- MS: Why k+1 ?
    2389 !-- MS: @Mohamed: emissivity belongs now to surface type with 3 different values for each
    2390 !--               surface element (due to tile approach).
     2388
    23912389                exn1 = ( hyp(k+1) / 100000.0_wp )**0.286_wp
    23922390
     
    24092407!--          element.
    24102408             ELSE
    2411 !
    2412 !--             Determine index offset between surface element and adjacent
    2413 !--             atmospheric grid point (depends on surface orientation).
    2414                 ioff = surf%ioff
    2415                 joff = surf%joff
    2416                 koff = surf%koff
    24172409
    24182410                DO  m = 1, surf%ns
     
    24292421!--                no weighted averaging is performed ( only one surface type per
    24302422!--                surface element ).
     2423!--                ATTENTION: when radiation interactions are switched on the
     2424!--                calculated fluxes below are not actually used as they are
     2425!--                overwritten in radiation_interaction.
    24312426                   IF ( ALLOCATED( surf%frac ) )  THEN
    24322427
     
    24412436                                           )                                      &
    24422437                                           * sigma_sb                             &
    2443                                            * ( pt(k+koff,j+joff,i+ioff) * exn )**4
    2444 
     2438                                           * ( surf%pt_surface(m) * exn )**4
    24452439
    24462440                      surf%rad_lw_out_change_0(m) =                               &
     
    24492443                                         + surf%frac(2,m) * surf%emissivity(2,m)  &
    24502444                                         ) * 3.0_wp * sigma_sb                    &
    2451                                          * ( pt(k+koff,j+joff,i+ioff) * exn )** 3
     2445                                         * ( surf%pt_surface(m) * exn )** 3
    24522446
    24532447                   ELSE
     
    24572451                      surf%rad_lw_out(m) = surf%emissivity(0,m)                   &
    24582452                                           * sigma_sb                             &
    2459                                            * ( pt(k+koff,j+joff,i+ioff) * exn )**4
    2460 
     2453                                           * ( surf%pt_surface(m) * exn )**4
    24612454
    24622455                      surf%rad_lw_out_change_0(m) = surf%emissivity(0,m)          &
    24632456                                           * 3.0_wp * sigma_sb                    &
    2464                                            * ( pt(k+koff,j+joff,i+ioff) * exn )** 3
     2457                                           * ( surf%pt_surface(m) * exn )** 3
    24652458
    24662459                   ENDIF
     
    26352628                                           )                                      &
    26362629                                         * sigma_sb                               &
    2637                                          * ( pt(k+koff,j+joff,i+ioff) * exn )**4
     2630                                         * ( surf%pt_surface(m) * exn )**4
    26382631
    26392632                      surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m)   &
     
    26532646                      surf%rad_lw_out(m) = surf%emissivity(0,m)                   &
    26542647                                         * sigma_sb                               &
    2655                                          * ( pt(k+koff,j+joff,i+ioff) * exn )**4
     2648                                         * ( surf%pt_surface(m) * exn )**4
    26562649
    26572650                      surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m)   &
     
    53635356                                            surfinlw(i) - surfoutlw(i)
    53645357!
    5365 !--        westward-facding
     5358!--        westward-facing
    53665359           ELSEIF ( surfl(1,i) == iwest_l )  THEN
    53675360              surf_lsm_v(3)%rad_sw_in(m)  = surfinsw(i)
  • palm/trunk/SOURCE/time_integration_spinup.f90

    r2718 r2723  
    2525! -----------------
    2626! $Id$
     27! Bugfix: array rad_sw_in no longer exists and is thus removed from RUN_CONTROL
     28! output.
     29! Added output of XY and 3D data during spinup.
     30! Bugfix: time step in LSM and USM was set to dt_3d instead of dt_spinup
     31!
     32! 2718 2018-01-02 08:49:38Z maronga
    2733! Corrected "Former revisions" section
    2834!
     
    6268               coupling_start_time, current_timestep_number,                   &
    6369               data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
    64                dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control,        &
    65                dt_spinup, humidity, intermediate_timestep_count,               &
     70               dt_averaging_input_pr, dt_dopr, dt_dots, dt_do2d_xy, dt_do3d, dt_run_control,        &
     71               dt_spinup, dt_3d, humidity, intermediate_timestep_count,               &
    6672               intermediate_timestep_count_max, land_surface,                  &
    67                simulated_time, simulated_time_chr,      &
    68                skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
    69                spinup_time, timestep_count, timestep_scheme, time_dopr,        &
    70                time_dopr_av, time_dots, time_run_control,                      &
     73               simulated_time, simulated_time_chr,                             &
     74               skip_time_dopr, skip_time_do2d_xy, skip_time_do3d, spinup, spinup_pt_amplitude, &
     75               spinup_pt_mean, spinup_time, timestep_count, timestep_scheme,   &
     76               time_dopr, time_dopr_av, time_dots, time_do2d_xy, time_do3d, time_run_control,           &
    7177               time_since_reference_point, urban_surface
    7278
     
    130136
    131137    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
     138    REAL(wp) ::  dt_save     !< temporary storage for time step
    132139                 
    133140    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
     
    137144    CALL exchange_horiz( pt, nbgp )   
    138145    pt_save = pt
     146
     147    dt_save = dt_3d
     148    dt_3d   = dt_spinup
    139149
    140150    CALL location_message( 'starting spinup-sequence', .TRUE. )
     
    299309               == intermediate_timestep_count_max )  THEN
    300310
    301                time_radiation = time_radiation + dt_spinup
     311               time_radiation = time_radiation + dt_3d
    302312
    303313             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
     
    327337!--    Increase simulation time and output times
    328338       current_timestep_number_spinup = current_timestep_number_spinup + 1
    329        simulated_time             = simulated_time   + dt_spinup
     339       simulated_time             = simulated_time   + dt_3d
    330340       simulated_time_chr         = time_to_string( simulated_time )
    331341       time_since_reference_point = simulated_time - coupling_start_time
    332342
    333343       IF ( data_output_during_spinup )  THEN
    334           time_dots          = time_dots        + dt_spinup
     344          IF ( simulated_time >= skip_time_do2d_xy )  THEN
     345             time_do2d_xy       = time_do2d_xy     + dt_3d
     346          ENDIF
     347          IF ( simulated_time >= skip_time_do3d    )  THEN
     348             time_do3d          = time_do3d        + dt_3d
     349          ENDIF
     350          time_dots          = time_dots        + dt_3d
    335351          IF ( simulated_time >= skip_time_dopr )  THEN
    336              time_dopr       = time_dopr        + dt_spinup
    337           ENDIF
    338           time_run_control   = time_run_control + dt_spinup
     352             time_dopr       = time_dopr        + dt_3d
     353          ENDIF
     354          time_run_control   = time_run_control + dt_3d
    339355
    340356!
     
    354370             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
    355371             simulated_time >= skip_time_dopr )  THEN
    356              time_dopr_av = time_dopr_av + dt_spinup
     372             time_dopr_av = time_dopr_av + dt_3d
    357373             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
    358374                do_sum = .TRUE.
    359375                time_dopr_av = MOD( time_dopr_av,                              &
    360                                MAX( dt_averaging_input_pr, dt_spinup ) )
     376                               MAX( dt_averaging_input_pr, dt_3d ) )
    361377             ENDIF
    362378          ENDIF
     
    367383          IF ( time_dopr >= dt_dopr )  THEN
    368384             IF ( dopr_n /= 0 )  CALL data_output_profiles
    369              time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
     385             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
    370386             time_dopr_av = 0.0_wp    ! due to averaging (see above)
    371387          ENDIF
     
    375391          IF ( time_dots >= dt_dots )  THEN
    376392             CALL data_output_tseries
    377              time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
    378           ENDIF
     393             time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
     394          ENDIF
     395
     396!
     397!--       2d-data output (cross-sections)
     398          IF ( time_do2d_xy >= dt_do2d_xy )  THEN
     399             CALL data_output_2d( 'xy', 0 )
     400             time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
     401          ENDIF
     402
     403!
     404!--       3d-data output (volume data)
     405          IF ( time_do3d >= dt_do3d )  THEN
     406             CALL data_output_3d( 0 )
     407             time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
     408          ENDIF
     409
    379410
    380411       ENDIF
     
    389420!           IF ( time_run_control >= dt_run_control )  THEN
    390421!              time_run_control = MOD( time_run_control,                         &
    391 !                                      MAX( dt_run_control, dt_spinup ) )
     422!                                      MAX( dt_run_control, dt_3d ) )
    392423!           ENDIF
    393424!        ENDIF
     
    408439!
    409440!--       Write some general information about the spinup in run control file
    410           WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_spinup, pt_spinup, rad_sw_in(0,nysg,nxlg)
     441          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_3d, pt_spinup
    411442!
    412443!--       Write buffer contents to disc immediately
     
    423454    pt_p(:,:,:) = pt_save
    424455
     456!
     457!-- Reset time step
     458    dt_3d = dt_save
     459
    425460    DEALLOCATE(pt_save)
    426461
     
    431466!-- Formats
    432467100 FORMAT (///'Spinup control output:'/  &
    433             '----------------------------------------'// &
    434             'ITER.  HH:MM:SS    DT   PT(z_MO)     SWD'/   &
    435             '----------------------------------------')
     468            '--------------------------------'// &
     469            'ITER.  HH:MM:SS    DT   PT(z_MO)'/   &
     470            '--------------------------------')
    436471101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
    437472
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r2720 r2723  
    2626! -----------------
    2727! $Id$
     28! Bugfix for spinups (end_time was increased twice in case of LSM + USM runs)
     29!
     30! 2720 2018-01-02 16:27:15Z kanani
    2831! Correction of comment
    2932!
     
    51075110       IF ( spinup_time > 0.0_wp )  THEN
    51085111          coupling_start_time = spinup_time
    5109           end_time = end_time + spinup_time
    51105112          IF ( spinup_pt_mean == 9999999.9_wp )  THEN
    51115113             spinup_pt_mean = pt_surface
    51125114          ENDIF
    5113           spinup = .TRUE.
     5115          IF ( .NOT. spinup )  THEN
     5116             end_time = end_time + spinup_time
     5117             spinup = .TRUE.
     5118          ENDIF
    51145119       ENDIF
    51155120
Note: See TracChangeset for help on using the changeset viewer.