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

Last change on this file since 3302 was 3274, checked in by knoop, 6 years ago

Modularization of all bulk cloud physics code components

  • Property svn:keywords set to Id
File size: 21.7 KB
Line 
1!> @file time_integration_spinup.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: time_integration_spinup.f90 3274 2018-09-24 15:42:55Z raasch $
27! Modularization of all bulk cloud physics code components
28!
29! 3241 2018-09-12 15:02:00Z raasch
30! unused variables removed
31!
32! 2983 2018-04-18 10:43:40Z suehring
33! Revise limitation of wall-adjacent velocity.
34!
35! 2934 2018-03-26 19:13:22Z suehring
36! Synchronize parent and child models after spinup.
37!
38! 2881 2018-03-13 16:24:40Z maronga
39! Added flag for switching on/off calculation of soil moisture
40!
41! 2818 2018-02-19 16:42:36Z maronga
42! Velocity components near walls/ground are now set to the profiles stored in
43! u_init and v_init. Activated soil moisture calculation during spinup.
44!
45! 2782 2018-02-02 11:51:10Z maronga
46! Bugfix and re-activation of homogeneous setting of velocity components
47! during spinup
48!
49! 2758 2018-01-17 12:55:21Z suehring
50! Comment out homogeneous setting of wind velocity as this will lead to zero
51! friction velocity and cause problems in MOST relationships.
52!
53! 2728 2018-01-09 07:03:53Z maronga
54! Set velocity componenets to homogeneous values during spinup
55!
56! 2724 2018-01-05 12:12:38Z maronga
57! Use dt_spinup for all active components during spinup
58!
59! 2723 2018-01-05 09:27:03Z maronga
60! Bugfix: array rad_sw_in no longer exists and is thus removed from RUN_CONTROL
61! output.
62! Added output of XY and 3D data during spinup.
63! Bugfix: time step in LSM and USM was set to dt_3d instead of dt_spinup
64!
65! 2718 2018-01-02 08:49:38Z maronga
66! Corrected "Former revisions" section
67!
68! 2696 2017-12-14 17:12:51Z kanani
69! Change in file header (GPL part)
70! Added radiation interactions (moved from USM) (MS)
71!
72! 2544 2017-10-13 18:09:32Z maronga
73! Date and time quantities are now read from date_and_time_mod
74!
75! 2299 2017-06-29 10:14:38Z maronga
76! Call of soil model adjusted to avoid prognostic equation for soil moisture
77! during spinup.
78! Better representation of diurnal cycle of near-surface temperature.
79! Excluded prognostic equation for soil moisture during spinup.
80! Added output of run control data for spinup.
81!
82! 2297 2017-06-28 14:35:57Z scharf
83! bugfixes
84!
85! 2296 2017-06-28 07:53:56Z maronga
86! Initial revision
87!
88!
89! Description:
90! ------------
91!> Integration in time of the non-atmospheric model components such as land
92!> surface model and urban surface model
93!------------------------------------------------------------------------------!
94 SUBROUTINE time_integration_spinup
95 
96    USE arrays_3d,                                                             &
97        ONLY:  pt, pt_p, u, u_init, v, v_init
98
99    USE control_parameters,                                                    &
100        ONLY:  averaging_interval_pr, calc_soil_moisture_during_spinup,        &
101               constant_diffusion, constant_flux_layer, coupling_start_time,   &
102               data_output_during_spinup, dopr_n, do_sum,                      &
103               dt_averaging_input_pr, dt_dopr, dt_dots, dt_do2d_xy, dt_do3d,   &
104               dt_spinup, dt_3d, humidity, intermediate_timestep_count,        &
105               intermediate_timestep_count_max, land_surface,                  &
106               simulated_time, simulated_time_chr, skip_time_dopr,             &
107               skip_time_do2d_xy, skip_time_do3d, spinup_pt_amplitude,         &
108               spinup_pt_mean, spinup_time, timestep_count, time_dopr,         &
109               time_dopr_av, time_dots, time_do2d_xy, time_do3d,               &
110               time_run_control, time_since_reference_point, urban_surface
111
112    USE cpulog,                                                                &
113        ONLY:  cpu_log, log_point, log_point_s
114
115    USE date_and_time_mod,                                                     &
116        ONLY: day_of_year_init, time_utc_init
117
118    USE indices,                                                               &
119        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
120
121
122    USE land_surface_model_mod,                                                &
123        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
124
125    USE pegrid
126
127    USE pmc_interface,                                                         &
128        ONLY:  nested_run
129
130    USE kinds
131
132    USE radiation_model_mod,                                                   &
133        ONLY:  force_radiation_call, radiation, radiation_control,             &
134               radiation_interaction, radiation_interactions, time_radiation
135
136    USE statistics,                                                            &
137        ONLY:  flow_statistics_called
138
139    USE surface_layer_fluxes_mod,                                              &
140        ONLY:  surface_layer_fluxes
141
142    USE surface_mod,                                                           &
143        ONLY :  surf_lsm_h, surf_lsm_v, surf_usm_h,    &
144                surf_usm_v
145
146    USE urban_surface_mod,                                                     &
147        ONLY:  usm_material_heat_model, usm_material_model,                    &
148               usm_surface_energy_balance, usm_swap_timelevel,                 &
149               usm_green_heat_model, usm_temperature_near_surface
150
151
152
153
154    IMPLICIT NONE
155
156    CHARACTER (LEN=9) ::  time_to_string          !<
157 
158    INTEGER(iwp) ::  i !< running index
159    INTEGER(iwp) ::  j !< running index
160    INTEGER(iwp) ::  k !< running index
161    INTEGER(iwp) ::  l !< running index
162    INTEGER(iwp) ::  m !< running index
163
164    INTEGER(iwp) :: current_timestep_number_spinup = 0  !< number if timestep during spinup
165 
166    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
167
168    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
169    REAL(wp) ::  dt_save     !< temporary storage for time step
170                 
171    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save  !< temporary storage of temperature
172    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_save   !< temporary storage of u wind component
173    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_save   !< temporary storage of v wind component
174
175
176!
177!-- Save 3D arrays because they are to be changed for spinup purpose
178    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
179    ALLOCATE( u_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
180    ALLOCATE( v_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
181
182    CALL exchange_horiz( pt, nbgp )   
183    CALL exchange_horiz( u,  nbgp ) 
184    CALL exchange_horiz( v,  nbgp ) 
185 
186    pt_save = pt
187    u_save  = u
188    v_save  = v
189
190!
191!-- Set the same wall-adjacent velocity to all grid points. The sign of the
192!-- original velocity field must be preserved because the surface schemes crash
193!-- otherwise. The precise reason is still unknown. A minimum velocity of 0.1
194!-- m/s is used to maintain turbulent transfer at the surface.
195    IF ( land_surface )  THEN
196       DO  m = 1, surf_lsm_h%ns
197          i   = surf_lsm_h%i(m)           
198          j   = surf_lsm_h%j(m)
199          k   = surf_lsm_h%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)
202       ENDDO
203
204       DO  l = 0, 3
205          DO  m = 1, surf_lsm_v(l)%ns
206             i   = surf_lsm_v(l)%i(m)           
207             j   = surf_lsm_v(l)%j(m)
208             k   = surf_lsm_v(l)%k(m)
209             u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
210             v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
211          ENDDO
212       ENDDO
213    ENDIF
214
215    IF ( urban_surface )  THEN
216       DO  m = 1, surf_usm_h%ns
217          i   = surf_usm_h%i(m)           
218          j   = surf_usm_h%j(m)
219          k   = surf_usm_h%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)
222       ENDDO
223
224       DO  l = 0, 3
225          DO  m = 1, surf_usm_v(l)%ns
226             i   = surf_usm_v(l)%i(m)           
227             j   = surf_usm_v(l)%j(m)
228             k   = surf_usm_v(l)%k(m)
229             u(k,j,i) = SIGN(1.0_wp,u_init(k)) * MAX( ABS( u_init(k) ),0.1_wp)
230             v(k,j,i) = SIGN(1.0_wp,v_init(k)) * MAX( ABS( v_init(k) ),0.1_wp)
231          ENDDO
232       ENDDO
233    ENDIF
234
235    CALL exchange_horiz( u,  nbgp )
236    CALL exchange_horiz( v,  nbgp )
237
238    dt_save = dt_3d
239    dt_3d   = dt_spinup
240
241    CALL location_message( 'starting spinup-sequence', .TRUE. )
242!
243!-- Start of the time loop
244    DO  WHILE ( simulated_time < spinup_time )
245
246       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
247   
248!
249!--    Start of intermediate step loop
250       intermediate_timestep_count = 0
251       DO  WHILE ( intermediate_timestep_count < &
252                   intermediate_timestep_count_max )
253
254          intermediate_timestep_count = intermediate_timestep_count + 1
255
256!
257!--       Set the steering factors for the prognostic equations which depend
258!--       on the timestep scheme
259          CALL timestep_scheme_steering
260
261
262!
263!--       Estimate a near-surface air temperature based on the position of the
264!--       sun and user input about mean temperature and amplitude. The time is
265!--       shifted by one hour to simulate a lag between air temperature and
266!--       incoming radiation
267          pt_spinup = spinup_pt_mean + spinup_pt_amplitude                     &
268             * solar_angle (time_utc_init + time_since_reference_point - 3600.0)
269
270!
271!--       Map air temperature to all grid points in the vicinity of a surface
272!--       element
273          IF ( land_surface )  THEN
274             DO  m = 1, surf_lsm_h%ns
275                i   = surf_lsm_h%i(m)           
276                j   = surf_lsm_h%j(m)
277                k   = surf_lsm_h%k(m)
278                pt(k,j,i) = pt_spinup
279             ENDDO
280
281             DO  l = 0, 3
282                DO  m = 1, surf_lsm_v(l)%ns
283                   i   = surf_lsm_v(l)%i(m)           
284                   j   = surf_lsm_v(l)%j(m)
285                   k   = surf_lsm_v(l)%k(m)
286                   pt(k,j,i) = pt_spinup
287                ENDDO
288             ENDDO
289          ENDIF
290
291          IF ( urban_surface )  THEN
292             DO  m = 1, surf_usm_h%ns
293                i   = surf_usm_h%i(m)           
294                j   = surf_usm_h%j(m)
295                k   = surf_usm_h%k(m)
296                pt(k,j,i) = pt_spinup
297             ENDDO
298
299             DO  l = 0, 3
300                DO  m = 1, surf_usm_v(l)%ns
301                   i   = surf_usm_v(l)%i(m)           
302                   j   = surf_usm_v(l)%j(m)
303                   k   = surf_usm_v(l)%k(m)
304                   pt(k,j,i) = pt_spinup
305                ENDDO
306             ENDDO
307          ENDIF
308
309          CALL exchange_horiz( pt,  nbgp )   
310
311
312!
313!--       Swap the time levels in preparation for the next time step.
314          timestep_count = timestep_count + 1
315     
316          IF ( land_surface )  THEN
317              CALL lsm_swap_timelevel ( 0 )
318          ENDIF
319
320          IF ( urban_surface )  THEN
321             CALL usm_swap_timelevel ( 0 )
322          ENDIF
323
324          IF ( land_surface )  THEN
325             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
326          ENDIF
327
328          IF ( urban_surface )  THEN
329             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
330          ENDIF
331         
332!
333!--       If required, compute virtual potential temperature
334          IF ( humidity )  THEN
335             CALL compute_vpt 
336          ENDIF 
337
338!
339!--       Compute the diffusion quantities
340          IF ( .NOT. constant_diffusion )  THEN
341
342!
343!--          First the vertical (and horizontal) fluxes in the surface
344!--          (constant flux) layer are computed
345             IF ( constant_flux_layer )  THEN
346                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
347                CALL surface_layer_fluxes
348                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
349             ENDIF
350
351!
352!--          If required, solve the energy balance for the surface and run soil
353!--          model. Call for horizontal as well as vertical surfaces.
354!--          The prognostic equation for soil moisure is switched off
355             IF ( land_surface )  THEN
356
357                CALL cpu_log( log_point(54), 'land_surface', 'start' )
358!
359!--             Call for horizontal upward-facing surfaces
360                CALL lsm_energy_balance( .TRUE., -1 )
361                CALL lsm_soil_model( .TRUE., -1, calc_soil_moisture_during_spinup )
362!
363!--             Call for northward-facing surfaces
364                CALL lsm_energy_balance( .FALSE., 0 )
365                CALL lsm_soil_model( .FALSE., 0, calc_soil_moisture_during_spinup )
366!
367!--             Call for southward-facing surfaces
368                CALL lsm_energy_balance( .FALSE., 1 )
369                CALL lsm_soil_model( .FALSE., 1, calc_soil_moisture_during_spinup )
370!
371!--             Call for eastward-facing surfaces
372                CALL lsm_energy_balance( .FALSE., 2 )
373                CALL lsm_soil_model( .FALSE., 2, calc_soil_moisture_during_spinup )
374!
375!--             Call for westward-facing surfaces
376                CALL lsm_energy_balance( .FALSE., 3 )
377                CALL lsm_soil_model( .FALSE., 3, calc_soil_moisture_during_spinup )
378
379                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
380             ENDIF
381
382!
383!--          If required, solve the energy balance for urban surfaces and run
384!--          the material heat model
385             IF (urban_surface) THEN
386                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
387                CALL usm_surface_energy_balance
388                IF ( usm_material_model )  THEN
389                   CALL usm_green_heat_model
390                   CALL usm_material_heat_model
391                ENDIF
392                IF ( urban_surface ) THEN
393                   CALL usm_temperature_near_surface
394                ENDIF
395                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
396             ENDIF
397
398          ENDIF
399
400!
401!--       If required, calculate radiative fluxes and heating rates
402          IF ( radiation .AND. intermediate_timestep_count                     &
403               == intermediate_timestep_count_max )  THEN
404
405               time_radiation = time_radiation + dt_3d
406
407             IF ( time_radiation >= dt_3d .OR. force_radiation_call )          &
408             THEN
409
410                CALL cpu_log( log_point(50), 'radiation', 'start' )
411
412                IF ( .NOT. force_radiation_call )  THEN
413                   time_radiation = time_radiation - dt_3d
414                ENDIF
415
416                CALL radiation_control
417
418                CALL cpu_log( log_point(50), 'radiation', 'stop' )
419
420                IF ( radiation_interactions )  THEN
421                   CALL cpu_log( log_point(75), 'radiation_interaction', 'start' )
422                   CALL radiation_interaction
423                   CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
424                ENDIF
425             ENDIF
426          ENDIF
427
428       ENDDO   ! Intermediate step loop
429
430!
431!--    Increase simulation time and output times
432       current_timestep_number_spinup = current_timestep_number_spinup + 1
433       simulated_time             = simulated_time   + dt_3d
434       simulated_time_chr         = time_to_string( simulated_time )
435       time_since_reference_point = simulated_time - coupling_start_time
436
437       IF ( data_output_during_spinup )  THEN
438          IF ( simulated_time >= skip_time_do2d_xy )  THEN
439             time_do2d_xy       = time_do2d_xy     + dt_3d
440          ENDIF
441          IF ( simulated_time >= skip_time_do3d    )  THEN
442             time_do3d          = time_do3d        + dt_3d
443          ENDIF
444          time_dots          = time_dots        + dt_3d
445          IF ( simulated_time >= skip_time_dopr )  THEN
446             time_dopr       = time_dopr        + dt_3d
447          ENDIF
448          time_run_control   = time_run_control + dt_3d
449
450!
451!--       Carry out statistical analysis and output at the requested output times.
452!--       The MOD function is used for calculating the output time counters (like
453!--       time_dopr) in order to regard a possible decrease of the output time
454!--       interval in case of restart runs
455
456!
457!--       Set a flag indicating that so far no statistics have been created
458!--       for this time step
459          flow_statistics_called = .FALSE.
460
461!
462!--       If required, call flow_statistics for averaging in time
463          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
464             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
465             simulated_time >= skip_time_dopr )  THEN
466             time_dopr_av = time_dopr_av + dt_3d
467             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
468                do_sum = .TRUE.
469                time_dopr_av = MOD( time_dopr_av,                              &
470                               MAX( dt_averaging_input_pr, dt_3d ) )
471             ENDIF
472          ENDIF
473          IF ( do_sum )  CALL flow_statistics
474
475!
476!--       Output of profiles
477          IF ( time_dopr >= dt_dopr )  THEN
478             IF ( dopr_n /= 0 )  CALL data_output_profiles
479             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_3d ) )
480             time_dopr_av = 0.0_wp    ! due to averaging (see above)
481          ENDIF
482
483!
484!--       Output of time series
485          IF ( time_dots >= dt_dots )  THEN
486             CALL data_output_tseries
487             time_dots = MOD( time_dots, MAX( dt_dots, dt_3d ) )
488          ENDIF
489
490!
491!--       2d-data output (cross-sections)
492          IF ( time_do2d_xy >= dt_do2d_xy )  THEN
493             CALL data_output_2d( 'xy', 0 )
494             time_do2d_xy = MOD( time_do2d_xy, MAX( dt_do2d_xy, dt_3d ) )
495          ENDIF
496
497!
498!--       3d-data output (volume data)
499          IF ( time_do3d >= dt_do3d )  THEN
500             CALL data_output_3d( 0 )
501             time_do3d = MOD( time_do3d, MAX( dt_do3d, dt_3d ) )
502          ENDIF
503
504
505       ENDIF
506
507!
508!--    Computation and output of run control parameters.
509!--    This is also done whenever perturbations have been imposed
510!        IF ( time_run_control >= dt_run_control  .OR.                           &
511!             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created )       &
512!        THEN
513!           CALL run_control
514!           IF ( time_run_control >= dt_run_control )  THEN
515!              time_run_control = MOD( time_run_control,                         &
516!                                      MAX( dt_run_control, dt_3d ) )
517!           ENDIF
518!        ENDIF
519
520       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
521
522
523!
524!--    Run control output
525       IF ( myid == 0 )  THEN
526!
527!--       If necessary, write header
528          IF ( .NOT. run_control_header_spinup )  THEN
529             CALL check_open( 15 )
530             WRITE ( 15, 100 )
531             run_control_header_spinup = .TRUE.
532          ENDIF
533!
534!--       Write some general information about the spinup in run control file
535          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_3d, pt_spinup
536!
537!--       Write buffer contents to disc immediately
538          FLUSH( 15 )
539       ENDIF
540
541
542
543    ENDDO   ! time loop
544
545!
546!-- Write back saved arrays to the 3D arrays
547    pt   = pt_save
548    pt_p = pt_save
549    u    = u_save
550    v    = v_save
551
552!
553!-- Reset time step
554    dt_3d = dt_save
555
556    DEALLOCATE(pt_save)
557    DEALLOCATE(u_save)
558    DEALLOCATE(v_save)
559
560#if defined( __parallel )
561    IF ( nested_run )  CALL MPI_BARRIER( MPI_COMM_WORLD, ierr )
562#endif
563
564    CALL location_message( 'finished spinup-sequence', .TRUE. )
565
566
567!
568!-- Formats
569100 FORMAT (///'Spinup control output:'/  &
570            '--------------------------------'// &
571            'ITER.  HH:MM:SS    DT   PT(z_MO)'/   &
572            '--------------------------------')
573101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
574
575 CONTAINS
576
577!
578!-- Returns the cosine of the solar zenith angle at a given time. This routine
579!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
580    FUNCTION solar_angle( local_time ) 
581
582       USE basic_constants_and_equations_mod,                                  &
583       ONLY:  pi
584     
585       USE kinds
586
587       USE radiation_model_mod,                                                &
588           ONLY:  decl_1, decl_2, decl_3, lat, lon
589
590       IMPLICIT NONE
591
592
593       REAL(wp) ::  solar_angle     !< cosine of the solar zenith angle
594
595       REAL(wp) ::  day             !< day of the year
596       REAL(wp) ::  declination     !< solar declination angle
597       REAL(wp) ::  hour_angle      !< solar hour angle
598       REAL(wp) ::  time_utc        !< current time in UTC
599       REAL(wp), INTENT(IN) :: local_time
600!
601!--    Calculate current day and time based on the initial values and simulation
602!--    time
603
604       day = day_of_year_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
605       time_utc = MOD(local_time, 86400.0_wp)
606
607
608!
609!--    Calculate solar declination and hour angle   
610       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
611       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
612
613!
614!--    Calculate cosine of solar zenith angle
615       solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
616                     * COS(hour_angle)
617
618
619    END FUNCTION solar_angle
620
621
622 END SUBROUTINE time_integration_spinup
Note: See TracBrowser for help on using the repository browser.