Changeset 2299 for palm/trunk/SOURCE/time_integration_spinup.f90
- Timestamp:
- Jun 29, 2017 10:14:38 AM (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/time_integration_spinup.f90
r2297 r2299 25 25 ! ----------------- 26 26 ! $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 27 34 ! bugfixes 28 35 ! 29 36 ! 2296 2017-06-28 07:53:56Z maronga 30 37 ! Initial revision 31 !32 !33 38 ! 34 39 ! … … 50 55 dt_spinup, humidity, intermediate_timestep_count, & 51 56 intermediate_timestep_count_max, land_surface, & 52 nr_timesteps_this_run,simulated_time, simulated_time_chr, &57 simulated_time, simulated_time_chr, & 53 58 skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean, & 54 59 spinup_time, timestep_count, timestep_scheme, time_dopr, & … … 67 72 68 73 USE land_surface_model_mod, & 69 ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel , &70 skip_time_do_lsm 71 72 USE pegrid74 ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel 75 76 USE pegrid, & 77 ONLY: myid 73 78 74 79 USE kinds … … 76 81 USE radiation_model_mod, & 77 82 ONLY: dt_radiation, force_radiation_call, radiation, & 78 radiation_control, skip_time_do_radiation, time_radiation83 radiation_control, rad_sw_in, time_radiation, time_utc_init 79 84 80 85 USE statistics, & … … 99 104 CHARACTER (LEN=9) :: time_to_string !< 100 105 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 106 113 114 LOGICAL :: run_control_header_spinup = .FALSE. !< flag parameter for steering whether the header information must be output 115 107 116 REAL(wp) :: pt_spinup !< temporary storage of temperature 108 117 109 118 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_save !< temporary storage of temperature 110 119 111 REAL(wp) :: day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude112 113 120 ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 114 121 115 CALL exchange_horiz( pt _p, nbgp )116 pt_save = pt _p122 CALL exchange_horiz( pt, nbgp ) 123 pt_save = pt 117 124 118 125 CALL location_message( 'starting spinup-sequence', .TRUE. ) … … 137 144 138 145 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 142 157 IF ( land_surface ) THEN 143 158 DO m = 1, surf_lsm_h%ns … … 145 160 j = surf_lsm_h%j(m) 146 161 k = surf_lsm_h%k(m) 147 pt _p(k,j,i) = pt_spinup162 pt(k,j,i) = pt_spinup 148 163 ENDDO 149 164 … … 153 168 j = surf_lsm_v(l)%j(m) 154 169 k = surf_lsm_v(l)%k(m) 155 pt _p(k,j,i) = pt_spinup170 pt(k,j,i) = pt_spinup 156 171 ENDDO 157 172 ENDDO … … 163 178 j = surf_usm_h%j(m) 164 179 k = surf_usm_h%k(m) 165 pt _p(k,j,i) = pt_spinup180 pt(k,j,i) = pt_spinup 166 181 ENDDO 167 182 … … 171 186 j = surf_usm_v(l)%j(m) 172 187 k = surf_usm_v(l)%k(m) 173 pt _p(k,j,i) = pt_spinup188 pt(k,j,i) = pt_spinup 174 189 ENDDO 175 190 ENDDO … … 217 232 ! 218 233 !-- 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 221 237 222 238 CALL cpu_log( log_point(54), 'land_surface', 'start' ) … … 224 240 !-- Call for horizontal upward-facing surfaces 225 241 CALL lsm_energy_balance( .TRUE., -1 ) 226 CALL lsm_soil_model( .TRUE., -1 )242 CALL lsm_soil_model( .TRUE., -1, .FALSE. ) 227 243 ! 228 244 !-- Call for northward-facing surfaces 229 245 CALL lsm_energy_balance( .FALSE., 0 ) 230 CALL lsm_soil_model( .FALSE., 0 )246 CALL lsm_soil_model( .FALSE., 0, .FALSE. ) 231 247 ! 232 248 !-- Call for southward-facing surfaces 233 249 CALL lsm_energy_balance( .FALSE., 1 ) 234 CALL lsm_soil_model( .FALSE., 1 )250 CALL lsm_soil_model( .FALSE., 1, .FALSE. ) 235 251 ! 236 252 !-- Call for eastward-facing surfaces 237 253 CALL lsm_energy_balance( .FALSE., 2 ) 238 CALL lsm_soil_model( .FALSE., 2 )254 CALL lsm_soil_model( .FALSE., 2, .FALSE. ) 239 255 ! 240 256 !-- Call for westward-facing surfaces 241 257 CALL lsm_energy_balance( .FALSE., 3 ) 242 CALL lsm_soil_model( .FALSE., 3 )258 CALL lsm_soil_model( .FALSE., 3, .FALSE. ) 243 259 244 260 CALL cpu_log( log_point(54), 'land_surface', 'stop' ) … … 262 278 !-- If required, calculate radiative fluxes and heating rates 263 279 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 266 281 267 282 time_radiation = time_radiation + dt_spinup … … 292 307 ! 293 308 !-- 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 296 310 simulated_time = simulated_time + dt_spinup 297 311 simulated_time_chr = time_to_string( simulated_time ) … … 350 364 !-- Computation and output of run control parameters. 351 365 !-- 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 ) 360 395 ENDIF 361 396 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 367 398 368 399 ENDDO ! time loop … … 375 406 DEALLOCATE(pt_save) 376 407 377 CALL location_message( 'finished time-stepping spinup', .TRUE. ) 408 CALL location_message( 'finished spinup-sequence', .TRUE. ) 409 410 411 ! 412 !-- Formats 413 100 FORMAT (///'Spinup control output:'/ & 414 '----------------------------------------'// & 415 'ITER. HH:MM:SS DT PT(z_MO) SWD'/ & 416 '----------------------------------------') 417 101 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 378 465 379 466 END SUBROUTINE time_integration_spinup
Note: See TracChangeset
for help on using the changeset viewer.