Changeset 4540 for palm/trunk/SOURCE/time_integration_spinup.f90
- Timestamp:
- May 18, 2020 3:23:29 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/time_integration_spinup.f90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/time_integration_spinup.f90
r4457 r4540 1 1 !> @file time_integration_spinup.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 4 ! 5 ! PALM is free software: you can redistribute it and/or modify it under the 6 ! terms of the GNU General Public License as published by the Free Software 7 ! Foundation, either version 3 of the License, or (at your option) any later 8 ! version. 9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 12 ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 13 ! 14 ! You should have received a copy of the GNU General Public License along with 15 ! PALM. If not, see <http://www.gnu.org/licenses/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 !------------------------------------------------------------------------------! 17 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: 21 ! ----------------- -21 ! ----------------- 22 22 ! 23 23 ! … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! use statement for exchange horiz added 28 ! 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 4457 2020-03-11 14:20:43Z raasch 30 ! Use statement for exchange horiz added 31 ! 29 32 ! 4444 2020-03-05 15:59:50Z raasch 30 ! bugfix: cpp-directives for serial mode added31 ! 33 ! Bugfix: cpp-directives for serial mode added 34 ! 32 35 ! 4360 2020-01-07 11:25:50Z suehring 33 36 ! Enable output of diagnostic quantities, e.g. 2-m temperature 34 ! 37 ! 35 38 ! 4227 2019-09-10 18:04:34Z gronemeier 36 ! implement new palm_date_time_mod37 ! 39 ! Implement new palm_date_time_mod 40 ! 38 41 ! 4223 2019-09-10 09:20:47Z gronemeier 39 42 ! Corrected "Former revisions" section 40 ! 43 ! 41 44 ! 4064 2019-07-01 05:33:33Z gronemeier 42 45 ! Moved call to radiation module out of intermediate time loop 43 ! 46 ! 44 47 ! 4023 2019-06-12 13:20:01Z maronga 45 48 ! Time stamps are now negative in run control output 46 ! 49 ! 47 50 ! 3885 2019-04-11 11:29:34Z kanani 48 ! Changes related to global restructuring of location messages and introduction 49 ! of additional debugmessages50 ! 51 ! Changes related to global restructuring of location messages and introduction of additional debug 52 ! messages 53 ! 51 54 ! 3766 2019-02-26 16:23:41Z raasch 52 ! unused variable removed53 ! 55 ! Unused variable removed 56 ! 54 57 ! 3719 2019-02-06 13:10:18Z kanani 55 ! Removed log_point(19,54,74,50,75), since they count together with same log 56 ! points in time_integration, impossible to separate the contributions.57 ! Instead, the entire spinup gets anindividual log_point in palm.f9058 ! 58 ! Removed log_point(19,54,74,50,75), since they count together with same log points in 59 ! time_integration, impossible to separate the contributions. Instead, the entire spinup gets an 60 ! individual log_point in palm.f90 61 ! 59 62 ! 3655 2019-01-07 16:51:22Z knoop 60 ! Removed call to calculation of near air (10 cm) potential temperature (now in 61 ! surface layer fluxes) 62 ! 63 ! Removed call to calculation of near air (10 cm) potential temperature (now in surface layer fluxes) 64 ! 63 65 ! 2296 2017-06-28 07:53:56Z maronga 64 66 ! Initial revision … … 67 69 ! Description: 68 70 ! ------------ 69 !> Integration in time of the non-atmospheric model components such as land 70 !> surface model and urban surface model71 !------------------------------------------------------------------------------ !71 !> Integration in time of the non-atmospheric model components such as land surface model and urban 72 !> surface model 73 !--------------------------------------------------------------------------------------------------! 72 74 SUBROUTINE time_integration_spinup 73 74 USE arrays_3d, & 75 ONLY: pt, pt_p, u, u_init, v, v_init 76 77 USE control_parameters, & 78 ONLY: averaging_interval_pr, calc_soil_moisture_during_spinup, & 79 constant_diffusion, constant_flux_layer, coupling_start_time, & 80 data_output_during_spinup, dopr_n, do_sum, & 81 dt_averaging_input_pr, dt_dopr, dt_dots, dt_do2d_xy, dt_do3d, & 82 dt_spinup, dt_3d, humidity, intermediate_timestep_count, & 83 intermediate_timestep_count_max, land_surface, & 84 simulated_time, simulated_time_chr, skip_time_dopr, & 85 skip_time_do2d_xy, skip_time_do3d, spinup_pt_amplitude, & 86 spinup_pt_mean, spinup_time, timestep_count, time_dopr, & 87 time_dopr_av, time_dots, time_do2d_xy, time_do3d, & 88 time_run_control, time_since_reference_point, urban_surface 89 90 USE cpulog, & 91 ONLY: cpu_log, log_point_s 92 93 USE diagnostic_output_quantities_mod, & 75 76 USE arrays_3d, & 77 ONLY: pt, & 78 pt_p, & 79 u, & 80 u_init, & 81 v, & 82 v_init 83 84 USE control_parameters, & 85 ONLY: averaging_interval_pr, & 86 calc_soil_moisture_during_spinup, & 87 constant_diffusion, & 88 constant_flux_layer, & 89 coupling_start_time, & 90 data_output_during_spinup, & 91 dopr_n, & 92 do_sum, & 93 dt_averaging_input_pr, & 94 dt_dopr, & 95 dt_dots, & 96 dt_do2d_xy, & 97 dt_do3d, & 98 dt_spinup, & 99 dt_3d, & 100 humidity, & 101 intermediate_timestep_count, & 102 intermediate_timestep_count_max, & 103 land_surface, & 104 simulated_time, & 105 simulated_time_chr, & 106 skip_time_dopr, & 107 skip_time_do2d_xy, & 108 skip_time_do3d, & 109 spinup_pt_amplitude, & 110 spinup_pt_mean, & 111 spinup_time, & 112 timestep_count, & 113 time_dopr, & 114 time_dopr_av, & 115 time_dots, & 116 time_do2d_xy, & 117 time_do3d, & 118 time_run_control, & 119 time_since_reference_point, & 120 urban_surface 121 122 USE cpulog, & 123 ONLY: cpu_log, & 124 log_point_s 125 126 USE diagnostic_output_quantities_mod, & 94 127 ONLY: doq_calculate 95 128 96 USE exchange_horiz_mod, &129 USE exchange_horiz_mod, & 97 130 ONLY: exchange_horiz 98 131 99 USE indices, & 100 ONLY: nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg 101 102 USE land_surface_model_mod, & 103 ONLY: lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel 132 USE indices, & 133 ONLY: nbgp, & 134 nzb, & 135 nzt, & 136 nysg, & 137 nyng, & 138 nxlg, & 139 nxrg 140 141 USE land_surface_model_mod, & 142 ONLY: lsm_energy_balance, & 143 lsm_soil_model, & 144 lsm_swap_timelevel 104 145 105 146 USE pegrid 106 147 107 148 #if defined( __parallel ) 108 USE pmc_interface, &149 USE pmc_interface, & 109 150 ONLY: nested_run 110 151 #endif … … 112 153 USE kinds 113 154 114 USE palm_date_time_mod, & 115 ONLY: get_date_time, seconds_per_hour 116 117 USE radiation_model_mod, & 118 ONLY: force_radiation_call, radiation, radiation_control, & 119 radiation_interaction, radiation_interactions, time_radiation 120 121 USE statistics, & 155 USE palm_date_time_mod, & 156 ONLY: get_date_time, & 157 seconds_per_hour 158 159 USE radiation_model_mod, & 160 ONLY: force_radiation_call, & 161 radiation, & 162 radiation_control, & 163 radiation_interaction, & 164 radiation_interactions, & 165 time_radiation 166 167 USE statistics, & 122 168 ONLY: flow_statistics_called 123 169 124 USE surface_layer_fluxes_mod, &170 USE surface_layer_fluxes_mod, & 125 171 ONLY: surface_layer_fluxes 126 172 127 USE surface_mod, & 128 ONLY : surf_lsm_h, surf_lsm_v, surf_usm_h, & 173 USE surface_mod, & 174 ONLY : surf_lsm_h, & 175 surf_lsm_v, surf_usm_h, & 129 176 surf_usm_v 130 177 131 USE urban_surface_mod, & 132 ONLY: usm_material_heat_model, usm_material_model, & 133 usm_surface_energy_balance, usm_swap_timelevel, & 178 USE urban_surface_mod, & 179 ONLY: usm_material_heat_model, & 180 usm_material_model, & 181 usm_surface_energy_balance, & 182 usm_swap_timelevel, & 134 183 usm_green_heat_model 135 184 … … 139 188 IMPLICIT NONE 140 189 141 CHARACTER (LEN=9) :: time_to_string !< 142 143 144 CHARACTER (LEN=1) :: sign_chr !< String containing '-' or ' ' 145 CHARACTER (LEN=9) :: time_since_reference_point_chr !< time since reference point, i.e., negative during spinup 146 147 INTEGER(iwp) :: i !< running index 148 INTEGER(iwp) :: j !< running index 149 INTEGER(iwp) :: k !< running index 150 INTEGER(iwp) :: l !< running index 151 INTEGER(iwp) :: m !< running index 152 153 INTEGER(iwp) :: current_timestep_number_spinup = 0 !< number if timestep during spinup 154 INTEGER(iwp) :: day_of_year !< day of the year 155 156 LOGICAL :: run_control_header_spinup = .FALSE. !< flag parameter for steering whether the header information must be output 157 190 CHARACTER(LEN=1) :: sign_chr !< String containing '-' or ' ' 191 CHARACTER(LEN=9) :: time_since_reference_point_chr !< time since reference point, i.e., negative during spinup 192 CHARACTER(LEN=9) :: time_to_string !< 193 194 195 INTEGER(iwp) :: current_timestep_number_spinup = 0 !< number if timestep during spinup 196 INTEGER(iwp) :: day_of_year !< day of the year 197 198 INTEGER(iwp) :: i !< running index 199 INTEGER(iwp) :: j !< running index 200 INTEGER(iwp) :: k !< running index 201 INTEGER(iwp) :: l !< running index 202 INTEGER(iwp) :: m !< running index 203 204 205 LOGICAL :: run_control_header_spinup = .FALSE. !< flag parameter for steering whether the header information must be output 206 207 208 REAL(wp) :: dt_save !< temporary storage for time step 158 209 REAL(wp) :: pt_spinup !< temporary storage of temperature 159 REAL(wp) :: dt_save !< temporary storage for time step160 210 REAL(wp) :: second_of_day !< second of the day 161 211 162 212 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: pt_save !< temporary storage of temperature 163 213 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: u_save !< temporary storage of u wind component … … 171 221 ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 172 222 173 CALL exchange_horiz( pt, nbgp ) 174 CALL exchange_horiz( u, nbgp ) 175 CALL exchange_horiz( v, nbgp ) 176 223 CALL exchange_horiz( pt, nbgp ) 224 CALL exchange_horiz( u, nbgp ) 225 CALL exchange_horiz( v, nbgp ) 226 177 227 pt_save = pt 178 228 u_save = u … … 180 230 181 231 ! 182 !-- Set the same wall-adjacent velocity to all grid points. The sign of the 183 !-- original velocity field must be preserved because the surface schemes crash 184 !-- otherwise. The precise reason is still unknown. A minimum velocity of 0.1 185 !-- m/s is used to maintain turbulent transfer at the surface. 232 !-- Set the same wall-adjacent velocity to all grid points. The sign of the original velocity field 233 !-- must be preserved because the surface schemes crash otherwise. The precise reason is still 234 !-- unknown. A minimum velocity of 0.1 m/s is used to maintain turbulent transfer at the surface. 186 235 IF ( land_surface ) THEN 187 236 DO m = 1, surf_lsm_h%ns 188 i = surf_lsm_h%i(m) 237 i = surf_lsm_h%i(m) 189 238 j = surf_lsm_h%j(m) 190 239 k = surf_lsm_h%k(m) 191 u(k,j,i) = SIGN( 1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)192 v(k,j,i) = SIGN( 1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)240 u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp) 241 v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp) 193 242 ENDDO 194 243 195 244 DO l = 0, 3 196 245 DO m = 1, surf_lsm_v(l)%ns 197 i = surf_lsm_v(l)%i(m) 246 i = surf_lsm_v(l)%i(m) 198 247 j = surf_lsm_v(l)%j(m) 199 248 k = surf_lsm_v(l)%k(m) 200 u(k,j,i) = SIGN( 1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)201 v(k,j,i) = SIGN( 1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)249 u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp) 250 v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp) 202 251 ENDDO 203 252 ENDDO … … 206 255 IF ( urban_surface ) THEN 207 256 DO m = 1, surf_usm_h%ns 208 i = surf_usm_h%i(m) 257 i = surf_usm_h%i(m) 209 258 j = surf_usm_h%j(m) 210 259 k = surf_usm_h%k(m) 211 u(k,j,i) = SIGN( 1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)212 v(k,j,i) = SIGN( 1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)260 u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp) 261 v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp) 213 262 ENDDO 214 263 215 264 DO l = 0, 3 216 265 DO m = 1, surf_usm_v(l)%ns 217 i = surf_usm_v(l)%i(m) 266 i = surf_usm_v(l)%i(m) 218 267 j = surf_usm_v(l)%j(m) 219 268 k = surf_usm_v(l)%k(m) 220 u(k,j,i) = SIGN( 1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)221 v(k,j,i) = SIGN( 1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)269 u(k,j,i) = SIGN( 1.0_wp, u_init(k) ) * MAX( ABS( u_init(k) ), 0.1_wp) 270 v(k,j,i) = SIGN( 1.0_wp, v_init(k) ) * MAX( ABS( v_init(k) ), 0.1_wp) 222 271 ENDDO 223 272 ENDDO 224 273 ENDIF 225 274 226 CALL exchange_horiz( u, 227 CALL exchange_horiz( v, 275 CALL exchange_horiz( u, nbgp ) 276 CALL exchange_horiz( v, nbgp ) 228 277 229 278 dt_save = dt_3d … … 236 285 237 286 CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' ) 238 287 239 288 ! 240 289 !-- Start of intermediate step loop 241 290 intermediate_timestep_count = 0 242 DO WHILE ( intermediate_timestep_count < & 243 intermediate_timestep_count_max ) 291 DO WHILE ( intermediate_timestep_count < intermediate_timestep_count_max ) 244 292 245 293 intermediate_timestep_count = intermediate_timestep_count + 1 246 294 247 295 ! 248 !-- Set the steering factors for the prognostic equations which depend 249 !-- on the timestep scheme 296 !-- Set the steering factors for the prognostic equations which depend on the timestep scheme 250 297 CALL timestep_scheme_steering 251 298 252 299 253 300 ! 254 !-- Estimate a near-surface air temperature based on the position of the 255 !-- sun and user input about mean temperature and amplitude. The time is 256 !-- shifted by one hour to simulate a lag between air temperature and 257 !-- incoming radiation 258 CALL get_date_time( simulated_time - spinup_time - seconds_per_hour, & 259 day_of_year=day_of_year, & 260 second_of_day=second_of_day ) 261 262 pt_spinup = spinup_pt_mean + spinup_pt_amplitude & 263 * solar_angle(day_of_year, second_of_day) 264 265 ! 266 !-- Map air temperature to all grid points in the vicinity of a surface 267 !-- element 301 !-- Estimate a near-surface air temperature based on the position of the sun and user input 302 !-- about mean temperature and amplitude. The time is shifted by one hour to simulate a lag 303 !-- between air temperature and incoming radiation. 304 CALL get_date_time( simulated_time - spinup_time - seconds_per_hour, & 305 day_of_year = day_of_year, second_of_day = second_of_day ) 306 307 pt_spinup = spinup_pt_mean + spinup_pt_amplitude * & 308 solar_angle( day_of_year, second_of_day ) 309 310 ! 311 !-- Map air temperature to all grid points in the vicinity of a surface element 268 312 IF ( land_surface ) THEN 269 313 DO m = 1, surf_lsm_h%ns 270 i = surf_lsm_h%i(m) 314 i = surf_lsm_h%i(m) 271 315 j = surf_lsm_h%j(m) 272 316 k = surf_lsm_h%k(m) … … 276 320 DO l = 0, 3 277 321 DO m = 1, surf_lsm_v(l)%ns 278 i = surf_lsm_v(l)%i(m) 322 i = surf_lsm_v(l)%i(m) 279 323 j = surf_lsm_v(l)%j(m) 280 324 k = surf_lsm_v(l)%k(m) … … 286 330 IF ( urban_surface ) THEN 287 331 DO m = 1, surf_usm_h%ns 288 i = surf_usm_h%i(m) 332 i = surf_usm_h%i(m) 289 333 j = surf_usm_h%j(m) 290 334 k = surf_usm_h%k(m) … … 297 341 DO l = 0, 3 298 342 DO m = 1, surf_usm_v(l)%ns 299 i = surf_usm_v(l)%i(m) 343 i = surf_usm_v(l)%i(m) 300 344 j = surf_usm_v(l)%j(m) 301 345 k = surf_usm_v(l)%k(m) … … 308 352 ENDIF 309 353 310 CALL exchange_horiz( pt, nbgp )354 CALL exchange_horiz( pt, nbgp ) 311 355 312 356 … … 314 358 !-- Swap the time levels in preparation for the next time step. 315 359 timestep_count = timestep_count + 1 316 360 317 361 IF ( land_surface ) THEN 318 362 CALL lsm_swap_timelevel ( 0 ) … … 324 368 325 369 IF ( land_surface ) THEN 326 CALL lsm_swap_timelevel ( MOD( timestep_count, 2 ) )370 CALL lsm_swap_timelevel ( MOD( timestep_count, 2 ) ) 327 371 ENDIF 328 372 329 373 IF ( urban_surface ) THEN 330 CALL usm_swap_timelevel ( MOD( timestep_count, 2 ) )331 ENDIF 332 333 ! 334 !-- If required, compute virtual potential temperature 335 IF ( humidity ) THEN 336 CALL compute_vpt 337 ENDIF 374 CALL usm_swap_timelevel ( MOD( timestep_count, 2 ) ) 375 ENDIF 376 377 ! 378 !-- If required, compute virtual potential temperature 379 IF ( humidity ) THEN 380 CALL compute_vpt 381 ENDIF 338 382 339 383 ! … … 342 386 343 387 ! 344 !-- First the vertical (and horizontal) fluxes in the surface 345 !-- (constant flux) layer arecomputed388 !-- First the vertical (and horizontal) fluxes in the surface (constant flux) layer are 389 !-- computed 346 390 IF ( constant_flux_layer ) THEN 347 391 CALL surface_layer_fluxes … … 349 393 350 394 ! 351 !-- If required, solve the energy balance for the surface and run soil 352 !-- model. Call for horizontal as well as vertical surfaces.353 !-- The prognostic equation for soil moisure isswitched off395 !-- If required, solve the energy balance for the surface and run soil model. Call for 396 !-- horizontal as well as vertical surfaces. The prognostic equation for soil moisure is 397 !-- switched off 354 398 IF ( land_surface ) THEN 355 399 … … 378 422 379 423 ! 380 !-- If required, solve the energy balance for urban surfaces and run 381 !-- the material heat model 424 !-- If required, solve the energy balance for urban surfaces and run the material heat model 382 425 IF (urban_surface) THEN 383 426 … … 417 460 !-- Increase simulation time and output times 418 461 current_timestep_number_spinup = current_timestep_number_spinup + 1 419 simulated_time = simulated_time + dt_3d420 simulated_time_chr = time_to_string( simulated_time )421 time_since_reference_point = simulated_time - coupling_start_time422 time_since_reference_point_chr = time_to_string( ABS( time_since_reference_point) )423 462 simulated_time = simulated_time + dt_3d 463 simulated_time_chr = time_to_string( simulated_time ) 464 time_since_reference_point = simulated_time - coupling_start_time 465 time_since_reference_point_chr = time_to_string( ABS( time_since_reference_point ) ) 466 424 467 IF ( time_since_reference_point < 0.0_wp ) THEN 425 468 sign_chr = '-' … … 427 470 sign_chr = ' ' 428 471 ENDIF 429 430 472 473 431 474 IF ( data_output_during_spinup ) THEN 432 475 IF ( simulated_time >= skip_time_do2d_xy ) THEN 433 time_do2d_xy 476 time_do2d_xy = time_do2d_xy + dt_3d 434 477 ENDIF 435 478 IF ( simulated_time >= skip_time_do3d ) THEN 436 time_do3d 437 ENDIF 438 time_dots = time_dots + dt_3d479 time_do3d = time_do3d + dt_3d 480 ENDIF 481 time_dots = time_dots + dt_3d 439 482 IF ( simulated_time >= skip_time_dopr ) THEN 440 time_dopr = time_dopr + dt_3d441 ENDIF 442 time_run_control = time_run_control + dt_3d483 time_dopr = time_dopr + dt_3d 484 ENDIF 485 time_run_control = time_run_control + dt_3d 443 486 444 487 ! 445 488 !-- Carry out statistical analysis and output at the requested output times. 446 !-- The MOD function is used for calculating the output time counters (like 447 !-- time_dopr) in order to regard a possible decrease of the output time 448 !-- interval in case of restart runs 449 450 ! 451 !-- Set a flag indicating that so far no statistics have been created 452 !-- for this time step 489 !-- The MOD function is used for calculating the output time counters (like time_dopr) in 490 !-- order to regard a possible decrease of the output time interval in case of restart runs. 491 492 ! 493 !-- Set a flag indicating that so far no statistics have been created for this time step 453 494 flow_statistics_called = .FALSE. 454 495 455 496 ! 456 497 !-- If required, call flow_statistics for averaging in time 457 IF ( averaging_interval_pr /= 0.0_wp .AND. & 458 ( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. & 459 simulated_time >= skip_time_dopr ) THEN 498 IF ( averaging_interval_pr /= 0.0_wp .AND. & 499 ( dt_dopr - time_dopr ) <= averaging_interval_pr .AND. & 500 simulated_time >= skip_time_dopr ) & 501 THEN 460 502 time_dopr_av = time_dopr_av + dt_3d 461 503 IF ( time_dopr_av >= dt_averaging_input_pr ) THEN 462 504 do_sum = .TRUE. 463 time_dopr_av = MOD( time_dopr_av, & 464 MAX( dt_averaging_input_pr, dt_3d ) ) 505 time_dopr_av = MOD( time_dopr_av, MAX( dt_averaging_input_pr, dt_3d ) ) 465 506 ENDIF 466 507 ENDIF … … 472 513 IF ( dopr_n /= 0 ) CALL data_output_profiles 473 514 time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) ) 474 time_dopr_av = 0.0_wp ! due to averaging (see above)515 time_dopr_av = 0.0_wp ! Due to averaging (see above) 475 516 ENDIF 476 517 … … 502 543 503 544 ! 504 !-- Computation and output of run control parameters. 505 !-- This is also done whenever perturbations have been imposed 506 ! IF ( time_run_control >= dt_run_control .OR. & 507 ! timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) & 508 ! THEN 545 !-- Computation and output of run control parameters. This is also done whenever perturbations 546 !-- have been imposed 547 ! IF ( time_run_control >= dt_run_control .OR. & 548 ! timestep_scheme(1:5) /= 'runge' .OR. disturbance_created ) THEN 509 549 ! CALL run_control 510 550 ! IF ( time_run_control >= dt_run_control ) THEN 511 ! time_run_control = MOD( time_run_control, & 512 ! MAX( dt_run_control, dt_3d ) ) 551 ! time_run_control = MOD( time_run_control, MAX( dt_run_control, dt_3d ) ) 513 552 ! ENDIF 514 553 ! ENDIF … … 529 568 ! 530 569 !-- Write some general information about the spinup in run control file 531 WRITE ( 15, 101 ) current_timestep_number_spinup, sign_chr, time_since_reference_point_chr, dt_3d, pt_spinup 570 WRITE ( 15, 101 ) current_timestep_number_spinup, sign_chr, & 571 time_since_reference_point_chr, dt_3d, pt_spinup 532 572 ! 533 573 !-- Write buffer contents to disc immediately … … 537 577 538 578 539 ENDDO ! time loop579 ENDDO ! Time loop 540 580 541 581 ! … … 563 603 ! 564 604 !-- Formats 565 100 FORMAT (///'Spinup control output:'/ & 566 '---------------------------------'// & 567 'ITER. HH:MM:SS DT PT(z_MO)'/ & 568 '---------------------------------') 605 100 FORMAT (///'Spinup control output:---------------------------------'// & 606 'ITER. HH:MM:SS DT PT(z_MO)---------------------------------') 569 607 101 FORMAT (I5,2X,A1,A9,1X,F6.2,3X,F6.2,2X,F6.2) 570 608 … … 572 610 573 611 ! 574 !-- Returns the cosine of the solar zenith angle at a given time. This routine 575 !-- is similar to that for calculation zenith (see radiation_model_mod.f90) 576 !> @todo Load function calc_zenith of radiation model instead of 577 !> rewrite the function here. 578 FUNCTION solar_angle( day_of_year, second_of_day ) 579 580 USE basic_constants_and_equations_mod, & 612 !-- Returns the cosine of the solar zenith angle at a given time. This routine is similar to that 613 !-- for calculation zenith (see radiation_model_mod.f90) 614 !> @todo Load function calc_zenith of radiation model instead of rewrite the function here. 615 FUNCTION solar_angle( day_of_year, second_of_day ) 616 617 USE basic_constants_and_equations_mod, & 581 618 ONLY: pi 582 619 583 620 USE kinds 584 621 585 USE radiation_model_mod, & 586 ONLY: decl_1, decl_2, decl_3, lat, lon 587 588 IMPLICIT NONE 622 USE radiation_model_mod, & 623 ONLY: decl_1, & 624 decl_2, & 625 decl_3, & 626 lat, & 627 lon 628 629 IMPLICIT NONE 589 630 590 631 591 632 INTEGER(iwp), INTENT(IN) :: day_of_year !< day of the year 592 633 593 REAL(wp) :: declination 594 REAL(wp) :: hour_angle 595 REAL(wp), INTENT(IN) :: second_of_day 596 REAL(wp) :: solar_angle 597 ! 598 !-- Calculate solar declination and hour angle 599 declination = ASIN( decl_1 * SIN( decl_2 * REAL(day_of_year, KIND=wp) - decl_3) )600 hour_angle = 2.0_wp * pi * ( second_of_day / 86400.0_wp) + lon - pi634 REAL(wp) :: declination !< solar declination angle 635 REAL(wp) :: hour_angle !< solar hour angle 636 REAL(wp), INTENT(IN) :: second_of_day !< current time of the day in UTC 637 REAL(wp) :: solar_angle !< cosine of the solar zenith angle 638 ! 639 !-- Calculate solar declination and hour angle 640 declination = ASIN( decl_1 * SIN( decl_2 * REAL( day_of_year, KIND = wp) - decl_3 ) ) 641 hour_angle = 2.0_wp * pi * ( second_of_day / 86400.0_wp ) + lon - pi 601 642 602 643 ! 603 644 !-- Calculate cosine of solar zenith angle 604 solar_angle = SIN( lat) * SIN(declination) + COS(lat) * COS(declination)&605 * COS(hour_angle)645 solar_angle = SIN( lat ) * SIN( declination ) + COS( lat ) * COS( declination ) * & 646 COS( hour_angle ) 606 647 607 648 END FUNCTION solar_angle
Note: See TracChangeset
for help on using the changeset viewer.
