Ignore:
Timestamp:
Jun 29, 2017 10:14:38 AM (7 years ago)
Author:
maronga
Message:

improvements for spinup mechanism

File:
1 edited

Legend:

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

    r2297 r2299  
    2525! -----------------
    2626! $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
    2734! bugfixes
    2835!
    2936! 2296 2017-06-28 07:53:56Z maronga
    3037! Initial revision
    31 !
    32 !
    3338!
    3439!
     
    5055               dt_spinup, humidity, intermediate_timestep_count,               &
    5156               intermediate_timestep_count_max, land_surface,                  &
    52                nr_timesteps_this_run, simulated_time, simulated_time_chr,      &
     57               simulated_time, simulated_time_chr,      &
    5358               skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
    5459               spinup_time, timestep_count, timestep_scheme, time_dopr,        &
     
    6772
    6873    USE land_surface_model_mod,                                                &
    69         ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
    70                skip_time_do_lsm
    71 
    72     USE pegrid
     74        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
     75
     76    USE pegrid,                                                                &
     77        ONLY:  myid
    7378
    7479    USE kinds
     
    7681    USE radiation_model_mod,                                                   &
    7782        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
    78                radiation_control, skip_time_do_radiation, time_radiation
     83               radiation_control, rad_sw_in, time_radiation, time_utc_init
    7984
    8085    USE statistics,                                                            &
     
    99104    CHARACTER (LEN=9) ::  time_to_string          !<
    100105 
    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
    106113 
     114    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
     115
    107116    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
    108117                 
    109118    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
    110119
    111     REAL(wp) ::  day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
    112 
    113120    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    114121
    115     CALL exchange_horiz( pt_p, nbgp )   
    116     pt_save = pt_p
     122    CALL exchange_horiz( pt, nbgp )   
     123    pt_save = pt
    117124
    118125    CALL location_message( 'starting spinup-sequence', .TRUE. )
     
    137144
    138145
    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
    142157          IF ( land_surface )  THEN
    143158             DO  m = 1, surf_lsm_h%ns
     
    145160                j   = surf_lsm_h%j(m)
    146161                k   = surf_lsm_h%k(m)
    147                 pt_p(k,j,i) = pt_spinup
     162                pt(k,j,i) = pt_spinup
    148163             ENDDO
    149164
     
    153168                   j   = surf_lsm_v(l)%j(m)
    154169                   k   = surf_lsm_v(l)%k(m)
    155                    pt_p(k,j,i) = pt_spinup
     170                   pt(k,j,i) = pt_spinup
    156171                ENDDO
    157172             ENDDO
     
    163178                j   = surf_usm_h%j(m)
    164179                k   = surf_usm_h%k(m)
    165                 pt_p(k,j,i) = pt_spinup
     180                pt(k,j,i) = pt_spinup
    166181             ENDDO
    167182
     
    171186                   j   = surf_usm_v(l)%j(m)
    172187                   k   = surf_usm_v(l)%k(m)
    173                    pt_p(k,j,i) = pt_spinup
     188                   pt(k,j,i) = pt_spinup
    174189                ENDDO
    175190             ENDDO
     
    217232!
    218233!--          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
    221237
    222238                CALL cpu_log( log_point(54), 'land_surface', 'start' )
     
    224240!--             Call for horizontal upward-facing surfaces
    225241                CALL lsm_energy_balance( .TRUE., -1 )
    226                 CALL lsm_soil_model( .TRUE., -1 )
     242                CALL lsm_soil_model( .TRUE., -1, .FALSE. )
    227243!
    228244!--             Call for northward-facing surfaces
    229245                CALL lsm_energy_balance( .FALSE., 0 )
    230                 CALL lsm_soil_model( .FALSE., 0 )
     246                CALL lsm_soil_model( .FALSE., 0, .FALSE. )
    231247!
    232248!--             Call for southward-facing surfaces
    233249                CALL lsm_energy_balance( .FALSE., 1 )
    234                 CALL lsm_soil_model( .FALSE., 1 )
     250                CALL lsm_soil_model( .FALSE., 1, .FALSE. )
    235251!
    236252!--             Call for eastward-facing surfaces
    237253                CALL lsm_energy_balance( .FALSE., 2 )
    238                 CALL lsm_soil_model( .FALSE., 2 )
     254                CALL lsm_soil_model( .FALSE., 2, .FALSE. )
    239255!
    240256!--             Call for westward-facing surfaces
    241257                CALL lsm_energy_balance( .FALSE., 3 )
    242                 CALL lsm_soil_model( .FALSE., 3 )
     258                CALL lsm_soil_model( .FALSE., 3, .FALSE. )
    243259
    244260                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
     
    262278!--       If required, calculate radiative fluxes and heating rates
    263279          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
    266281
    267282               time_radiation = time_radiation + dt_spinup
     
    292307!
    293308!--    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
    296310       simulated_time             = simulated_time   + dt_spinup
    297311       simulated_time_chr         = time_to_string( simulated_time )
     
    350364!--    Computation and output of run control parameters.
    351365!--    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 )
    360395       ENDIF
    361396
    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
    367398
    368399    ENDDO   ! time loop
     
    375406    DEALLOCATE(pt_save)
    376407
    377     CALL location_message( 'finished time-stepping spinup', .TRUE. )
     408    CALL location_message( 'finished spinup-sequence', .TRUE. )
     409
     410
     411!
     412!-- Formats
     413100 FORMAT (///'Spinup control output:'/  &
     414            '----------------------------------------'// &
     415            'ITER.  HH:MM:SS    DT   PT(z_MO)     SWD'/   &
     416            '----------------------------------------')
     417101 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
    378465
    379466 END SUBROUTINE time_integration_spinup
Note: See TracChangeset for help on using the changeset viewer.