source: palm/trunk/SOURCE/time_integration_spinup.f90 @ 2296

Last change on this file since 2296 was 2296, checked in by maronga, 4 years ago

added new spinup mechanism for surface/radiation models

  • Property svn:keywords set to Id
File size: 13.0 KB
Line 
1!> @file time_integration_spinup.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration_spinup.f90 2296 2017-06-28 07:53:56Z maronga $
27! Initial revision
28!
29!
30!
31!
32! Description:
33! ------------
34!> Integration in time of the non-atmospheric model components such as land
35!> surface model and urban surface model
36!------------------------------------------------------------------------------!
37 SUBROUTINE time_integration_spinup
38 
39    USE arrays_3d,                                                             &
40        ONLY:  pt, pt_p
41
42    USE control_parameters,                                                    &
43        ONLY:  averaging_interval_pr, constant_diffusion, constant_flux_layer, coupling_start_time, current_timestep_number, data_output_during_spinup, disturbance_created, dopr_n, do_sum, dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control, dt_spinup, humidity, intermediate_timestep_count, intermediate_timestep_count_max, land_surface, nr_timesteps_this_run, simulated_time, simulated_time_chr, skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean, spinup_time, timestep_count, timestep_scheme, time_dopr, time_dopr_av, time_dots, time_run_control, time_since_reference_point, urban_surface
44
45    USE constants,                                                             &
46        ONLY:  pi
47
48    USE cpulog,                                                                &
49        ONLY:  cpu_log, log_point, log_point_s
50
51    USE indices,                                                               &
52        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
53
54
55    USE land_surface_model_mod,                                                &
56        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
57               skip_time_do_lsm
58
59    USE pegrid
60
61    USE kinds
62
63    USE radiation_model_mod,                                                   &
64        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
65               radiation_control, skip_time_do_radiation, time_radiation
66
67    USE statistics,                                                            &
68        ONLY:  flow_statistics_called
69
70    USE surface_layer_fluxes_mod,                                              &
71        ONLY:  surface_layer_fluxes
72
73    USE surface_mod,                                                        &
74        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
75                surf_usm_v
76
77    USE urban_surface_mod,                                                     &
78        ONLY:  usm_material_heat_model, usm_material_model,                    &
79               usm_radiation, usm_surface_energy_balance, usm_swap_timelevel
80
81
82
83
84    IMPLICIT NONE
85
86    CHARACTER (LEN=9) ::  time_to_string          !<
87 
88    INTEGER(iwp) :: i
89    INTEGER(iwp) :: j
90    INTEGER(iwp) :: k
91    INTEGER(iwp) :: l
92    INTEGER(iwp) :: m
93 
94    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
95                 
96    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
97
98    REAL(wp) ::  day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
99
100    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
101
102    CALL exchange_horiz( pt_p, nbgp )   
103    pt_save = pt_p
104
105    CALL location_message( 'starting spinup-sequence', .TRUE. )
106!
107!-- Start of the time loop
108    DO  WHILE ( simulated_time < spinup_time )
109
110       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
111   
112!
113!--    Start of intermediate step loop
114       intermediate_timestep_count = 0
115       DO  WHILE ( intermediate_timestep_count < &
116                   intermediate_timestep_count_max )
117
118          intermediate_timestep_count = intermediate_timestep_count + 1
119
120!
121!--       Set the steering factors for the prognostic equations which depend
122!--       on the timestep scheme
123          CALL timestep_scheme_steering
124
125
126!!!!      Set new values of pt_p here
127          pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
128
129          IF ( land_surface )  THEN
130             DO  m = 1, surf_lsm_h%ns
131                i   = surf_lsm_h%i(m)           
132                j   = surf_lsm_h%j(m)
133                k   = surf_lsm_h%k(m)
134                pt_p(k,j,i) = pt_spinup
135             ENDDO
136
137             DO  l = 0, 3
138                DO  m = 1, surf_lsm_v(l)%ns
139                   i   = surf_lsm_v(l)%i(m)           
140                   j   = surf_lsm_v(l)%j(m)
141                   k   = surf_lsm_v(l)%k(m)
142                   pt_p(k,j,i) = pt_spinup
143                ENDDO
144             ENDDO
145          ENDIF
146
147          IF ( urban_surface )  THEN
148             DO  m = 1, surf_usm_h%ns
149                i   = surf_usm_h%i(m)           
150                j   = surf_usm_h%j(m)
151                k   = surf_usm_h%k(m)
152                pt_p(k,j,i) = pt_spinup
153             ENDDO
154
155             DO  l = 0, 3
156                DO  m = 1, surf_usm_v(l)%ns
157                   i   = surf_usm_v(l)%i(m)           
158                   j   = surf_usm_v(l)%j(m)
159                   k   = surf_usm_v(l)%k(m)
160                   pt_p(k,j,i) = pt_spinup
161                ENDDO
162             ENDDO
163          ENDIF
164
165!
166!--       Swap the time levels in preparation for the next time step.
167          timestep_count = timestep_count + 1
168     
169          IF ( land_surface )  THEN
170              CALL lsm_swap_timelevel ( 0 )
171          ENDIF
172
173          IF ( urban_surface )  THEN
174             CALL usm_swap_timelevel ( 0 )
175          ENDIF
176
177          IF ( land_surface )  THEN
178             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
179          ENDIF
180
181          IF ( urban_surface )  THEN
182             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
183          ENDIF
184         
185!
186!--       If required, compute virtual potential temperature
187          IF ( humidity )  THEN
188             CALL compute_vpt
189          ENDIF
190
191!
192!--       Compute the diffusion quantities
193          IF ( .NOT. constant_diffusion )  THEN
194
195!
196!--          First the vertical (and horizontal) fluxes in the surface
197!--          (constant flux) layer are computed
198             IF ( constant_flux_layer )  THEN
199                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
200                CALL surface_layer_fluxes
201                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
202             ENDIF
203
204!
205!--          If required, solve the energy balance for the surface and run soil
206!--          model. Call for horizontal as well as vertical surfaces
207             IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
208
209                CALL cpu_log( log_point(54), 'land_surface', 'start' )
210!
211!--             Call for horizontal upward-facing surfaces
212                CALL lsm_energy_balance( .TRUE., -1 )
213                CALL lsm_soil_model( .TRUE., -1 )
214!
215!--             Call for northward-facing surfaces
216                CALL lsm_energy_balance( .FALSE., 0 )
217                CALL lsm_soil_model( .FALSE., 0 )
218!
219!--             Call for southward-facing surfaces
220                CALL lsm_energy_balance( .FALSE., 1 )
221                CALL lsm_soil_model( .FALSE., 1 )
222!
223!--             Call for eastward-facing surfaces
224                CALL lsm_energy_balance( .FALSE., 2 )
225                CALL lsm_soil_model( .FALSE., 2 )
226!
227!--             Call for westward-facing surfaces
228                CALL lsm_energy_balance( .FALSE., 3 )
229                CALL lsm_soil_model( .FALSE., 3 )
230
231                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
232             ENDIF
233
234!
235!--          If required, solve the energy balance for urban surfaces and run
236!--          the material heat model
237             IF (urban_surface) THEN
238                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
239                CALL usm_surface_energy_balance
240                IF ( usm_material_model )  THEN
241                   CALL usm_material_heat_model
242                ENDIF
243                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
244             ENDIF
245
246          ENDIF
247
248!
249!--       If required, calculate radiative fluxes and heating rates
250          IF ( radiation .AND. intermediate_timestep_count                     &
251               == intermediate_timestep_count_max .AND. simulated_time >    &
252               skip_time_do_radiation )  THEN
253
254               time_radiation = time_radiation + dt_spinup
255
256             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
257             THEN
258
259                CALL cpu_log( log_point(50), 'radiation', 'start' )
260
261                IF ( .NOT. force_radiation_call )  THEN
262                   time_radiation = time_radiation - dt_radiation
263                ENDIF
264
265                CALL radiation_control
266
267                CALL cpu_log( log_point(50), 'radiation', 'stop' )
268
269                IF (urban_surface)  THEN
270                   CALL cpu_log( log_point(75), 'usm_radiation', 'start' )
271                   CALL usm_radiation
272                   CALL cpu_log( log_point(75), 'usm_radiation', 'stop' )
273                ENDIF
274             ENDIF
275          ENDIF
276
277       ENDDO   ! Intermediate step loop
278
279!
280!--    Increase simulation time and output times
281       nr_timesteps_this_run      = nr_timesteps_this_run + 1
282       current_timestep_number    = current_timestep_number + 1
283       simulated_time             = simulated_time   + dt_spinup
284       simulated_time_chr         = time_to_string( simulated_time )
285       time_since_reference_point = simulated_time - coupling_start_time
286
287       IF ( data_output_during_spinup )  THEN
288          time_dots          = time_dots        + dt_spinup
289          IF ( simulated_time >= skip_time_dopr )  THEN
290             time_dopr       = time_dopr        + dt_spinup
291          ENDIF
292          time_run_control   = time_run_control + dt_spinup
293
294!
295!--       Carry out statistical analysis and output at the requested output times.
296!--       The MOD function is used for calculating the output time counters (like
297!--       time_dopr) in order to regard a possible decrease of the output time
298!--       interval in case of restart runs
299
300!
301!--       Set a flag indicating that so far no statistics have been created
302!--       for this time step
303          flow_statistics_called = .FALSE.
304
305!
306!--       If required, call flow_statistics for averaging in time
307          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
308             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
309             simulated_time >= skip_time_dopr )  THEN
310             time_dopr_av = time_dopr_av + dt_spinup
311             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
312                do_sum = .TRUE.
313                time_dopr_av = MOD( time_dopr_av,                              &
314                               MAX( dt_averaging_input_pr, dt_spinup ) )
315             ENDIF
316          ENDIF
317          IF ( do_sum )  CALL flow_statistics
318
319!
320!--       Output of profiles
321          IF ( time_dopr >= dt_dopr )  THEN
322             IF ( dopr_n /= 0 )  CALL data_output_profiles
323             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
324             time_dopr_av = 0.0_wp    ! due to averaging (see above)
325          ENDIF
326
327!
328!--       Output of time series
329          IF ( time_dots >= dt_dots )  THEN
330             CALL data_output_tseries
331             time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
332          ENDIF
333
334       ENDIF
335
336!
337!--    Computation and output of run control parameters.
338!--    This is also done whenever perturbations have been imposed
339       IF ( time_run_control >= dt_run_control  .OR.                     &
340            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
341       THEN
342          CALL run_control
343          IF ( time_run_control >= dt_run_control )  THEN
344             time_run_control = MOD( time_run_control, &
345                                     MAX( dt_run_control, dt_spinup ) )
346          ENDIF
347       ENDIF
348
349       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
350
351       IF ( myid == 0 )  THEN
352          PRINT*, time_since_reference_point
353       ENDIF
354
355    ENDDO   ! time loop
356
357!
358!-- Write back saved temperature to the 3D arrays
359    pt(:,:,:)   = pt_save
360    pt_p(:,:,:) = pt_save
361
362    DEALLOCATE(pt_save)
363
364    CALL location_message( 'finished time-stepping spinup', .TRUE. )
365
366 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.