SCALE-RM
scale_land_dyn_bucket.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
19  use scale_debug
20  use scale_file_external_input, only: &
22  !-----------------------------------------------------------------------------
23  implicit none
24  private
25  !-----------------------------------------------------------------------------
26  !
27  !++ Public procedure
28  !
29  public :: land_dyn_bucket_setup
30  public :: land_dyn_bucket
31 
32  !-----------------------------------------------------------------------------
33  !
34  !++ Public parameters & variables
35  !
36  !-----------------------------------------------------------------------------
37  !
38  !++ Private procedure
39  !
40  !-----------------------------------------------------------------------------
41  !
42  !++ Private parameters & variables
43  !
44  logical, private :: land_dyn_bucket_update_bottom_temp = .false. ! Is LAND_TEMP updated in the lowest level?
45  logical, private :: land_dyn_bucket_update_bottom_water = .false. ! Is LAND_WATER updated in the lowest level?
46 
47  logical, private :: land_dyn_bucket_nudging = .false. ! Is nudging for land physics used?
48  real(DP), private :: land_dyn_bucket_nudging_tau = 0.0_dp ! time constant for nudging [sec]
49  character(len=H_SHORT), private :: land_dyn_bucket_nudging_tau_unit = "SEC"
50  character(len=H_LONG), private :: land_dyn_bucket_nudging_basename(file_external_input_file_limit) = ''
51  logical, private :: land_dyn_bucket_nudging_enable_periodic_year = .false.
52  logical, private :: land_dyn_bucket_nudging_enable_periodic_month = .false.
53  logical, private :: land_dyn_bucket_nudging_enable_periodic_day = .false.
54  integer, private :: land_dyn_bucket_nudging_step_fixed = 0
55  real(RP), private :: land_dyn_bucket_nudging_offset = 0.0_rp
56  real(RP), private :: land_dyn_bucket_nudging_defval ! = UNDEF
57  logical, private :: land_dyn_bucket_nudging_check_coordinates = .true.
58  integer, private :: land_dyn_bucket_nudging_step_limit = 0
59 
60  real(RP), private :: water_denscs
61  real(DP), private :: land_dyn_bucket_nudging_tausec
62  !-----------------------------------------------------------------------------
63 contains
64  !-----------------------------------------------------------------------------
66  subroutine land_dyn_bucket_setup
67  use scale_prc, only: &
68  prc_abort
69  use scale_const, only: &
70  dwatr => const_dwatr, &
71  cl => const_cl
72  use scale_calendar, only: &
74  use scale_file_external_input, only: &
76  use scale_const, only: &
77  undef => const_undef
78  implicit none
79 
80  namelist / param_land_dyn_bucket / &
81  land_dyn_bucket_nudging, &
82  land_dyn_bucket_nudging_tau, &
83  land_dyn_bucket_nudging_tau_unit, &
84  land_dyn_bucket_nudging_basename, &
85  land_dyn_bucket_nudging_enable_periodic_year, &
86  land_dyn_bucket_nudging_enable_periodic_month, &
87  land_dyn_bucket_nudging_enable_periodic_day, &
88  land_dyn_bucket_nudging_step_fixed, &
89  land_dyn_bucket_nudging_offset, &
90  land_dyn_bucket_nudging_defval, &
91  land_dyn_bucket_nudging_check_coordinates, &
92  land_dyn_bucket_nudging_step_limit, &
93  land_dyn_bucket_update_bottom_temp, &
94  land_dyn_bucket_update_bottom_water
95 
96  integer :: ierr
97  !---------------------------------------------------------------------------
98 
99  log_newline
100  log_info("LAND_DYN_BUCKET_setup",*) 'Setup'
101 
102  land_dyn_bucket_nudging_defval = undef
103 
104  !--- read namelist
105  rewind(io_fid_conf)
106  read(io_fid_conf,nml=param_land_dyn_bucket,iostat=ierr)
107  if( ierr < 0 ) then !--- missing
108  log_info("LAND_DYN_BUCKET_setup",*) 'Not found namelist. Default used.'
109  elseif( ierr > 0 ) then !--- fatal error
110  log_error("LAND_DYN_BUCKET_setup",*) 'Not appropriate names in namelist PARAM_LAND_DYN_BUCKET. Check!'
111  call prc_abort
112  endif
113  log_nml(param_land_dyn_bucket)
114 
115  if ( land_dyn_bucket_nudging ) then
116  call calendar_unit2sec( land_dyn_bucket_nudging_tausec, land_dyn_bucket_nudging_tau, land_dyn_bucket_nudging_tau_unit )
117 
118  log_info("LAND_DYN_BUCKET_setup",*) 'Use nudging for LAND physics : ON'
119  log_info("LAND_DYN_BUCKET_setup",*) 'Relaxation time Tau [sec] : ', land_dyn_bucket_nudging_tausec
120 
121  if ( land_dyn_bucket_nudging_tausec == 0.0_rp ) then
122  log_info("LAND_DYN_BUCKET_setup",*) 'Tau=0 means that LST is completely replaced by the external data.'
123  endif
124 
125  if ( land_dyn_bucket_nudging_basename(1) == '' ) then
126  log_error("LAND_DYN_BUCKET_setup",*) 'LAND_DYN_BUCKET_nudging_basename is necessary !!'
127  call prc_abort
128  end if
129 
130  call file_external_input_regist( land_dyn_bucket_nudging_basename(:), & ! [IN]
131  'LAND_TEMP', & ! [IN]
132  'LXY', & ! [IN]
133  land_dyn_bucket_nudging_enable_periodic_year, & ! [IN]
134  land_dyn_bucket_nudging_enable_periodic_month, & ! [IN]
135  land_dyn_bucket_nudging_enable_periodic_day, & ! [IN]
136  land_dyn_bucket_nudging_step_fixed, & ! [IN]
137  land_dyn_bucket_nudging_offset, & ! [IN]
138  land_dyn_bucket_nudging_defval, & ! [IN]
139  land_dyn_bucket_nudging_check_coordinates, & ! [IN]
140  land_dyn_bucket_nudging_step_limit ) ! [IN]
141 
142  call file_external_input_regist( land_dyn_bucket_nudging_basename(:), & ! [IN]
143  'LAND_WATER', & ! [IN]
144  'LXY', & ! [IN]
145  land_dyn_bucket_nudging_enable_periodic_year, & ! [IN]
146  land_dyn_bucket_nudging_enable_periodic_month, & ! [IN]
147  land_dyn_bucket_nudging_enable_periodic_day, & ! [IN]
148  land_dyn_bucket_nudging_step_fixed, & ! [IN]
149  land_dyn_bucket_nudging_offset, & ! [IN]
150  land_dyn_bucket_nudging_defval, & ! [IN]
151  land_dyn_bucket_nudging_check_coordinates, & ! [IN]
152  land_dyn_bucket_nudging_step_limit ) ! [IN]
153 
154  log_info("LAND_DYN_BUCKET_setup",*) 'Use nudging for Land physics: ON'
155  else
156  log_info("LAND_DYN_BUCKET_setup",*) 'Use nudging for Land physics: OFF'
157  end if
158 
159  water_denscs = dwatr * cl
160 
161  log_newline
162  log_info("LAND_DYN_BUCKET_setup",*) 'Update soil temperature of bottom layer? : ', land_dyn_bucket_update_bottom_temp
163  log_info("LAND_DYN_BUCKET_setup",*) 'Update soil moisture of bottom layer? : ', land_dyn_bucket_update_bottom_water
164 
165  return
166  end subroutine land_dyn_bucket_setup
167 
168  !-----------------------------------------------------------------------------
170  subroutine land_dyn_bucket( &
171  LKMAX, LKS, LKE, LIA, LIS, LIE, LJA, LJS, LJE, &
172  TEMP_t, WATER_t, &
173  WaterLimit, &
174  ThermalCond, &
175  HeatCapacity, &
176  WaterDiff, &
177  SFLX_GH, SFLX_water, SFLX_ice, &
178  fact_land, CDZ, &
179  dt, NOWDAYSEC, &
180  TEMP, WATER, &
181  RUNOFF )
182  use scale_const, only: &
183  tem00 => const_tem00, &
184  dwatr => const_dwatr, &
185  emelt => const_emelt
186  use scale_prc, only: &
187  prc_abort
188  use scale_file_external_input, only: &
189  file_external_input_update
190  use scale_matrix, only: &
191  matrix_solver_tridiagonal
192  implicit none
193  integer, intent(in) :: LKMAX, LKS, LKE
194  integer, intent(in) :: LIA, LIS, LIE
195  integer, intent(in) :: LJA, LJS, LJE
196 
197  real(RP), intent(in) :: TEMP_t (lkmax,lia,lja)
198  real(RP), intent(in) :: WATER_t (lkmax,lia,lja)
199  real(RP), intent(in) :: WaterLimit (lia,lja)
200  real(RP), intent(in) :: ThermalCond (lia,lja)
201  real(RP), intent(in) :: HeatCapacity(lia,lja)
202  real(RP), intent(in) :: WaterDiff (lia,lja)
203  real(RP), intent(in) :: SFLX_GH (lia,lja) ! positive for downward
204  real(RP), intent(in) :: SFLX_water (lia,lja) ! positive for downward
205  real(RP), intent(in) :: SFLX_ice (lia,lja) ! positive for downward
206  real(RP), intent(in) :: fact_land (lia,lja)
207  real(RP), intent(in) :: CDZ (lkmax)
208  real(DP), intent(in) :: dt
209  real(DP), intent(in) :: NOWDAYSEC
210 
211  real(RP), intent(inout) :: TEMP (lkmax,lia,lja)
212  real(RP), intent(inout) :: WATER(lkmax,lia,lja)
213 
214  real(RP), intent(out) :: RUNOFF(lia,lja)
215 
216  logical :: solve_matrix
217  logical :: error
218 
219  real(RP) :: TEMP1 (lkmax,lia,lja)
220  real(RP) :: WATER1(lkmax,lia,lja)
221 
222  real(RP) :: LAND_DENSCS(lkmax,lia,lja)
223  real(RP) :: ThermalDiff(lkmax,lia,lja)
224 
225  real(RP) :: U(lkmax,lia,lja)
226  real(RP) :: M(lkmax,lia,lja)
227  real(RP) :: L(lkmax,lia,lja)
228  real(RP) :: V(lkmax,lia,lja)
229 
230  real(RP) :: NDG_TEMP (lkmax,lia,lja)
231  real(RP) :: NDG_WATER(lkmax,lia,lja)
232 
233  real(RP) :: Q
234 
235  integer :: k, i, j
236  !---------------------------------------------------------------------------
237 
238  log_progress(*) 'land / dynamics / bucket'
239 
240  if( land_dyn_bucket_nudging ) then
241 
242  call file_external_input_update( &
243  'LAND_TEMP', & ! (in)
244  nowdaysec, & ! (in)
245  temp1, & ! (out)
246  error ) ! (out)
247  if ( error ) then
248  log_error("LAND_DYN_BUCKET",*) 'Requested data is not found!'
249  call prc_abort
250  end if
251 
252  call file_external_input_update( &
253  'LAND_WATER', & ! (in)
254  nowdaysec, & ! (in)
255  water1, & ! (out)
256  error ) ! (out)
257  if ( error ) then
258  log_error("LAND_DYN_BUCKET",*) 'Requested data is not found!'
259  call prc_abort
260  end if
261 
262  if( land_dyn_bucket_nudging_tau > 0.0_rp ) then
263  ! nudging is used
264  solve_matrix = .true.
265 
266  !$omp parallel do
267  do j = ljs,lje
268  do i = lis,lie
269  do k = lks,lke
270  ndg_temp(k,i,j) = ( temp1(k,i,j) - temp(k,i,j) ) / land_dyn_bucket_nudging_tausec * dt
271  ndg_water(k,i,j) = ( water1(k,i,j) - water(k,i,j) ) / land_dyn_bucket_nudging_tausec * dt
272  end do
273  end do
274  end do
275 
276  else
277  ! replace data to reference
278  solve_matrix = .false.
279 
280  end if
281 
282  else
283  ! nudging is NOT used
284  solve_matrix = .true.
285 
286  ndg_temp(:,:,:) = 0.0_rp
287  ndg_water(:,:,:) = 0.0_rp
288 
289  end if
290 
291  if( solve_matrix ) then
292 
293  ! Solve diffusion of soil moisture (tridiagonal matrix)
294  do j = ljs, lje
295  do i = lis, lie
296  l(lks,i,j) = 0.0_rp
297  u(lks,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(lks) * ( cdz(lks) + cdz(lks+1) ) ) * dt
298  l(lke,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(lke) * ( cdz(lke) + cdz(lke-1) ) ) * dt
299  u(lke,i,j) = 0.0_rp
300 
301  m(lks,i,j) = 1.0_rp - l(lks,i,j) - u(lks,i,j)
302  m(lke,i,j) = 1.0_rp - l(lke,i,j) - u(lke,i,j)
303  end do
304  end do
305 
306  do j = ljs, lje
307  do i = lis, lie
308  do k = lks+1, lke-1
309  l(k,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(k) * ( cdz(k) + cdz(k-1) ) ) * dt
310  u(k,i,j) = -2.0_rp * waterdiff(i,j) / ( cdz(k) * ( cdz(k) + cdz(k+1) ) ) * dt
311  m(k,i,j) = 1.0_rp - l(k,i,j) - u(k,i,j)
312  end do
313  end do
314  end do
315 
316  ! input from atmosphere
317  do j = ljs, lje
318  do i = lis, lie
319  v(lks,i,j) = water(lks,i,j) + ndg_water(lks,i,j) &
320  + ( sflx_water(i,j) + sflx_ice(i,j) ) / ( cdz(lks) * dwatr ) * dt
321  end do
322  end do
323 
324  do j = ljs, lje
325  do i = lis, lie
326  do k = lks+1, lke
327  v(k,i,j) = water(k,i,j) + ndg_water(k,i,j)
328  end do
329  end do
330  end do
331 
332  call matrix_solver_tridiagonal( lkmax, & ! [IN]
333  lia, lis, lie, & ! [IN]
334  lja, ljs, lje, & ! [IN]
335  u(:,:,:), & ! [IN]
336  m(:,:,:), & ! [IN]
337  l(:,:,:), & ! [IN]
338  v(:,:,:), & ! [IN]
339  water1(:,:,:) ) ! [OUT]
340 
341  if ( .not. land_dyn_bucket_update_bottom_water ) then
342  do j = ljs, lje
343  do i = lis, lie
344  water1(lke,i,j) = water(lke,i,j)
345  end do
346  end do
347  endif
348 
349  ! runoff of soil moisture (vertical sum)
350  do j = ljs, lje
351  do i = lis, lie
352  runoff(i,j) = 0.0_rp
353  do k = lks, lke
354  runoff(i,j) = runoff(i,j) + max( water1(k,i,j) - waterlimit(i,j), 0.0_rp ) * cdz(k) * dwatr
355  water1(k,i,j) = min( water1(k,i,j), waterlimit(i,j) )
356  end do
357  end do
358  end do
359 
360  ! estimate thermal diffusivity
361  do j = ljs, lje
362  do i = lis, lie
363  do k = lks, lke
364  land_denscs(k,i,j) = ( 1.0_rp - waterlimit(i,j) ) * heatcapacity(i,j) + water_denscs * water1(k,i,j)
365  thermaldiff(k,i,j) = thermalcond(i,j) / land_denscs(k,i,j)
366  end do
367  end do
368  end do
369 
370  ! Solve diffusion of soil temperature (tridiagonal matrix)
371  do j = ljs, lje
372  do i = lis, lie
373  l(lks,i,j) = 0.0_rp
374  u(lks,i,j) = -2.0_rp * thermaldiff(lks,i,j) / ( cdz(lks) * ( cdz(lks) + cdz(lks+1) ) ) * dt
375  l(lke,i,j) = -2.0_rp * thermaldiff(lke,i,j) / ( cdz(lke) * ( cdz(lke) + cdz(lke-1) ) ) * dt
376  u(lke,i,j) = 0.0_rp
377 
378  m(lks,i,j) = 1.0_rp - l(lks,i,j) - u(lks,i,j)
379  m(lke,i,j) = 1.0_rp - l(lke,i,j) - u(lke,i,j)
380  end do
381  end do
382 
383  do j = ljs, lje
384  do i = lis, lie
385  do k = lks+1, lke-1
386  l(k,i,j) = -2.0_rp * thermaldiff(k,i,j) / ( cdz(k) * ( cdz(k) + cdz(k-1) ) ) * dt
387  u(k,i,j) = -2.0_rp * thermaldiff(k,i,j) / ( cdz(k) * ( cdz(k) + cdz(k+1) ) ) * dt
388  m(k,i,j) = 1.0_rp - l(k,i,j) - u(k,i,j)
389  end do
390  end do
391  end do
392 
393  ! input from atmosphere
394  do j = ljs, lje
395  do i = lis, lie
396  ! tentative (MUST consider ice water explicitly)
397  q = sflx_gh(i,j)
398  if ( temp(lks,i,j) >= tem00 ) q = q - sflx_ice(i,j) * emelt
399  v(lks,i,j) = temp(lks,i,j) + ndg_temp(lks,i,j) + temp_t(lks,i,j) &
400  + q / ( land_denscs(lks,i,j) * cdz(lks) ) * dt
401  end do
402  end do
403 
404  do j = ljs, lje
405  do i = lis, lie
406  do k = lks+1, lke
407  v(k,i,j) = temp(k,i,j) + ndg_temp(k,i,j) + temp_t(k,i,j)
408  end do
409  end do
410  end do
411 
412  call matrix_solver_tridiagonal( lkmax, & ! [IN]
413  lia, lis, lie, & ! [IN]
414  lja, ljs, lje, & ! [IN]
415  u(:,:,:), & ! [IN]
416  m(:,:,:), & ! [IN]
417  l(:,:,:), & ! [IN]
418  v(:,:,:), & ! [IN]
419  temp1(:,:,:) ) ! [OUT]
420 
421  if ( .not. land_dyn_bucket_update_bottom_temp ) then
422  do j = ljs, lje
423  do i = lis, lie
424  temp1(lke,i,j) = temp(lke,i,j)
425  end do
426  end do
427  endif
428 
429  end if
430 
431  ! calculate tendency
432  do j = ljs, lje
433  do i = lis, lie
434  if( fact_land(i,j) > 0.0_rp ) then
435  do k = lks, lke
436  temp(k,i,j) = temp1(k,i,j)
437  water(k,i,j) = water1(k,i,j)
438  end do
439  end if
440  end do
441  end do
442 
443  return
444  end subroutine land_dyn_bucket
445 
446 end module scale_land_dyn_bucket
module DEBUG
Definition: scale_debug.F90:11
subroutine, public calendar_unit2sec(second, value, unit)
Convert several units to second.
real(rp), parameter, public const_dwatr
density of water [kg/m3]
Definition: scale_const.F90:82
real(rp), parameter, public const_cl
specific heat (liquid water) [J/kg/K]
Definition: scale_const.F90:66
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:90
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
module file / external_input
subroutine, public file_external_input_regist(basename, varname, axistype, enable_periodic_year, enable_periodic_month, enable_periodic_day, step_fixed, offset, defval, check_coordinates, step_limit, exist)
Regist data.
real(rp), public const_undef
Definition: scale_const.F90:41
module land / dynamics / bucket
integer, parameter, public file_external_input_file_limit
limit of file (for one item)
module PROCESS
Definition: scale_prc.F90:11
module MATRIX
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
module CONSTANT
Definition: scale_const.F90:11
real(rp), parameter, public const_emelt
Definition: scale_const.F90:72
subroutine, public land_dyn_bucket(LKMAX, LKS, LKE, LIA, LIS, LIE, LJA, LJS, LJE, TEMP_t, WATER_t, WaterLimit, ThermalCond, HeatCapacity, WaterDiff, SFLX_GH, SFLX_water, SFLX_ice, fact_land, CDZ, dt, NOWDAYSEC, TEMP, WATER, RUNOFF)
Physical processes for land submodel.
module profiler
Definition: scale_prof.F90:11
module PRECISION
module CALENDAR
module STDIO
Definition: scale_io.F90:10
subroutine, public land_dyn_bucket_setup
Setup.