MONC
gridmanager.F90
Go to the documentation of this file.
1 
2 module gridmanager_mod
3  use datadefn_mod, only : string_length
4  use monc_component_mod, only : component_descriptor_type
5  use state_mod, only : model_state_type
6  use optionsdatabase_mod, only : options_get_integer, options_get_logical, options_get_real, &
7  options_get_logical_array, options_get_real_array, options_get_string_array, options_get_array_size, options_get_string
8  use grids_mod, only : vertical_grid_configuration_type, x_index, y_index, z_index
9  use logging_mod, only : log_info, log_error, log_master_log, log_log
10  use conversions_mod, only : conv_to_string
11  use datadefn_mod, only : default_precision
12  use science_constants_mod, only : von_karman_constant, z0, z0th, cp, r_over_cp, r, g, rlvap_over_cp
13  use saturation_mod, only : qsaturation, dqwsatdt
14  use q_indices_mod, only: get_q_index, standard_q_names
15  use interpolation_mod, only: piecewise_linear_1d
16 
17  implicit none
18 
19 #ifndef TEST_MODE
20  private
21 #endif
22 
23  ! 1 = No adjustment
24  ! 2 = ensure P0=PSF by adjusting THREF profile by constant factor
25  ! 3 = ensure P0=PSF by adjusting PSF (not advised)
26  ! 4 = ensure P0=PSF by adjusting PTOP
27  integer, parameter :: anelastic_profile_mode=4
28  real, parameter :: default_spacing = 1.e9
29  real(kind=default_precision), dimension(:,:), allocatable :: qinit
30 
31  public gridmanager_get_descriptor
32 
33 contains
34 
37  type(component_descriptor_type) function gridmanager_get_descriptor()
38  gridmanager_get_descriptor%name="grid_manager"
39  gridmanager_get_descriptor%version=0.1
40  gridmanager_get_descriptor%initialisation=>initialise_callback
41  gridmanager_get_descriptor%finalisation=>finalise_callback
42  end function gridmanager_get_descriptor
43 
47  subroutine initialise_callback(current_state)
48  type(model_state_type), target, intent(inout) :: current_state
49 
50  integer :: dimensions
51 
52  if (.not. current_state%initialised) then
53  call log_log(log_error, "Must initialise the model state_mod before constructing the grid properties")
54  end if
55 
56  call initialise_horizontalgrid_configuration_types(current_state)
57  call initialise_verticalgrid_configuration_type(current_state)
58  dimensions=1
59  if (current_state%global_grid%active(x_index)) dimensions = dimensions+1
60  if (current_state%global_grid%active(y_index)) dimensions = dimensions+1
61  call log_master_log(log_info, trim(conv_to_string(dimensions))//"D system; z="//&
62  trim(conv_to_string(current_state%global_grid%size(z_index)))//", y="//&
63  trim(conv_to_string(current_state%global_grid%size(y_index)))//", x="//&
64  trim(conv_to_string(current_state%global_grid%size(x_index))))
65  end subroutine initialise_callback
66 
69  subroutine finalise_callback(current_state)
70  type(model_state_type), target, intent(inout) :: current_state
71 
72  type(vertical_grid_configuration_type) :: vertical_grid
73 
74  vertical_grid=current_state%global_grid%configuration%vertical
75 
76  deallocate(vertical_grid%z, vertical_grid%zn, vertical_grid%dz, vertical_grid%dzn, vertical_grid%czb, vertical_grid%cza, &
77  vertical_grid%czg, vertical_grid%czh, vertical_grid%rdz, vertical_grid%rdzn, vertical_grid%tzc1, vertical_grid%tzc2,&
78  vertical_grid%tzd1, vertical_grid%tzd2, vertical_grid%thref, vertical_grid%theta_init, vertical_grid%temp_init, &
79  vertical_grid%rh_init, vertical_grid%tref, &
80  vertical_grid%prefn, vertical_grid%pdiff, vertical_grid%prefrcp, vertical_grid%rprefrcp, vertical_grid%rho, &
81  vertical_grid%rhon, vertical_grid%tstarpr, vertical_grid%qsat, vertical_grid%dqsatdt, vertical_grid%qsatfac, &
82  vertical_grid%dthref, vertical_grid%rneutml, vertical_grid%rneutml_sq, vertical_grid%buoy_co, &
83  vertical_grid%u_init, vertical_grid%v_init, vertical_grid%q_init, &
84  vertical_grid%q_rand, vertical_grid%theta_rand, vertical_grid%w_subs, vertical_grid%w_rand, &
85  vertical_grid%q_force, vertical_grid%theta_force, vertical_grid%u_force, vertical_grid%v_force &
86  )
87  end subroutine finalise_callback
88 
91  subroutine initialise_verticalgrid_configuration_type(current_state)
92  type(model_state_type), intent(inout) :: current_state
93 
94  call allocate_vertical_grid_data(current_state%global_grid%configuration%vertical, &
95  current_state%global_grid%size(z_index), current_state%number_q_fields )
96  call set_up_and_smooth_grid(current_state%global_grid%configuration%vertical, &
97  current_state%global_grid%configuration%vertical%kgd, current_state%global_grid%configuration%vertical%hgd, &
98  size(current_state%global_grid%configuration%vertical%kgd), current_state%global_grid%size(z_index), &
99  current_state%global_grid%top(z_index), options_get_integer(current_state%options_database, "nsmth"), &
100  current_state%origional_vertical_grid_setup, current_state%continuation_run)
101  call set_vertical_reference_profile(current_state, current_state%global_grid%configuration%vertical, &
102  current_state%global_grid%size(z_index))
103 
104  end subroutine initialise_verticalgrid_configuration_type
105 
110  subroutine set_vertical_reference_profile(current_state, vertical_grid, kkp)
111  type(model_state_type), intent(inout) :: current_state
112  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
113  integer, intent(in) :: kkp
114 
115  integer :: k
116 
117  call calculate_initial_profiles(current_state, vertical_grid)
118  call set_up_vertical_reference_properties(current_state, vertical_grid, current_state%global_grid%size(z_index))
119  call set_anelastic_pressure(current_state)
120  !
121  call set_qv_init_from_rh(current_state)
122 
123  do k=2,kkp-1
124  ! for diffusion onto p-level from below
125  vertical_grid%czb(k)=(vertical_grid%rho(k-1)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k))
126  ! for diffusion onto p-level from above
127  vertical_grid%cza(k)=(vertical_grid%rho(k)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k+1))
128  vertical_grid%czg(k)=-vertical_grid%czb(k)-vertical_grid%cza(k)
129  if (k .gt. 2) vertical_grid%czh(k)=vertical_grid%czb(k)*vertical_grid%cza(k-1)
130  end do
131  do k=2,kkp-1
132  ! advection onto p-level from below
133  vertical_grid%tzc1(k)=0.25_default_precision*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k)
134  ! advection onto p-level from above
135  vertical_grid%tzc2(k)=0.25_default_precision*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k)
136  end do
137  do k=2,kkp-1
138  ! advection onto w-level (K) from below
139  vertical_grid%tzd1(k)=0.25_default_precision*vertical_grid%rdzn(k+1)*vertical_grid%rhon(k)/vertical_grid%rho(k)
140  ! advection onto w-level (K) from above
141  vertical_grid%tzd2(k)=0.25_default_precision*vertical_grid%rdzn(k+1)*vertical_grid%rhon(k+1)/vertical_grid%rho(k)
142  end do
143  k=kkp
144  vertical_grid%czb(k)=(vertical_grid%rho(k-1)/vertical_grid%rhon(k))/(vertical_grid%dz(k)*vertical_grid%dzn(k))
145  vertical_grid%cza(k)=0.0_default_precision
146  vertical_grid%czg(k)=-vertical_grid%czb(k)
147  vertical_grid%czh(k)=vertical_grid%czb(k)*vertical_grid%cza(k-1)
148  vertical_grid%tzc2(k)=0.25_default_precision*vertical_grid%rdz(k)*vertical_grid%rho(k)/vertical_grid%rhon(k)
149  vertical_grid%tzc1(k)=0.25_default_precision*vertical_grid%rdz(k)*vertical_grid%rho(k-1)/vertical_grid%rhon(k)
150  vertical_grid%czn=vertical_grid%dzn(2)*0.5_default_precision
151  vertical_grid%zlogm=log(1.0_default_precision+vertical_grid%zn(2)/z0)
152  vertical_grid%zlogth=log((vertical_grid%zn(2)+z0)/z0th)
153  vertical_grid%vk_on_zlogm=von_karman_constant/vertical_grid%zlogm
154  call setup_reference_state_liquid_water_temperature_and_saturation(&
155  current_state, vertical_grid, current_state%global_grid%size(z_index))
156  call calculate_mixing_length_for_neutral_case(current_state, vertical_grid, current_state%global_grid%size(z_index))
157  call set_buoyancy_coefficient(current_state, vertical_grid, current_state%global_grid%size(z_index))
158  end subroutine set_vertical_reference_profile
159 
164  subroutine calculate_initial_profiles(current_state, vertical_grid)
165  type(model_state_type), intent(inout) :: current_state
166  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
167 
168  integer :: nq_init ! The number of q fields to initialize
169  integer :: nzq ! The number of input levels for q_init
170  integer :: i,j,n, k ! loop counters
171  integer :: iq ! temporary q varible index
172 
173  real(kind=default_precision), dimension(:,:), allocatable :: f_init_pl_q ! Initial node values for q variables
174  real(kind=default_precision), dimension(:), allocatable :: z_init_pl_q ! Initial node height values for q variables
175  real(kind=default_precision), dimension(:), allocatable :: f_init_pl_theta ! Initial node values for potential temperature variable
176  real(kind=default_precision), dimension(:), allocatable :: z_init_pl_theta ! Initial node height values for potential temperature variable
177  real(kind=default_precision), dimension(:), allocatable :: f_init_pl_u ! Initial node values for u variable
178  real(kind=default_precision), dimension(:), allocatable :: z_init_pl_u ! Initial node height values for u variable
179  real(kind=default_precision), dimension(:), allocatable :: f_init_pl_v ! Initial node values for v variable
180  real(kind=default_precision), dimension(:), allocatable :: z_init_pl_v ! Initial node height values for v variable
181 
182  real(kind=default_precision), dimension(:), allocatable :: f_thref ! Initial node values for thref
183  real(kind=default_precision), dimension(:), allocatable :: z_thref ! Initial node height values for thref
184 
185  logical :: l_init_pl_u ! if .true. then initialize u field
186  logical :: l_init_pl_v ! if .true. then initialize v field
187  logical :: l_init_pl_theta ! if .true. then initialize potential temperature field
188  logical :: l_init_pl_rh ! if .true. then initialize relative humidity field
189  logical :: l_init_pl_q ! if .true. then initialize q fields
190  logical :: l_thref ! if .true. then initialize thref profile (overrides thref0)
191  logical :: l_matchthref ! if .true. then initialize thref to be the same as theta_init
192 
193  character(len=STRING_LENGTH), dimension(:), allocatable :: names_init_pl_q ! names of q variables to initialize
194 
195  real(kind=default_precision), allocatable :: f_init_pl_q_tmp(:) !temporary 1D storage of initial q field
196  real(kind=default_precision), allocatable :: zgrid(:) ! z grid to use in interpolation
197 
198  real(kind=default_precision) :: zztop ! top of the domain
199  real(kind=default_precision) :: qsat
200 
201  allocate(zgrid(current_state%local_grid%local_domain_end_index(z_index)))
202 
203  zztop = current_state%global_grid%top(z_index)
204 
205  ! Initialize everything to zero. This won't make sense for theta.
206  vertical_grid%q_init = 0.0_default_precision
207  vertical_grid%u_init = 0.0_default_precision
208  vertical_grid%v_init = 0.0_default_precision
209  vertical_grid%theta_init = 0.0_default_precision
210 
211  l_init_pl_theta=options_get_logical(current_state%options_database, "l_init_pl_theta")
212  l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh")
213  l_init_pl_q=options_get_logical(current_state%options_database, "l_init_pl_q")
214  if (l_init_pl_q) then
215  allocate(names_init_pl_q(options_get_array_size(current_state%options_database, "names_init_pl_q")))
216  call options_get_string_array(current_state%options_database, "names_init_pl_q", names_init_pl_q)
217  do n = 1,size(names_init_pl_q)
218  if (trim(names_init_pl_q(n)) .eq. 'vapour' .and. l_init_pl_rh) then
219  call log_master_log(log_error, "Initialisation of vapour and RH - STOP")
220  endif
221  enddo
222  end if
223  l_init_pl_u=options_get_logical(current_state%options_database, "l_init_pl_u")
224  l_init_pl_v=options_get_logical(current_state%options_database, "l_init_pl_v")
225 
226  l_thref=options_get_logical(current_state%options_database, "l_thref")
227  l_matchthref=options_get_logical(current_state%options_database, "l_matchthref")
228 
229  if (l_thref)then
230  if (.not. l_matchthref)then
231  allocate(z_thref(options_get_array_size(current_state%options_database, "z_thref")), &
232  f_thref(options_get_array_size(current_state%options_database, "f_thref")))
233  call options_get_real_array(current_state%options_database, "z_thref", z_thref)
234  call options_get_real_array(current_state%options_database, "f_thref", f_thref)
235  call check_top(zztop, z_thref(size(z_thref)), 'z_thref')
236  zgrid=current_state%global_grid%configuration%vertical%zn(:)
237  call piecewise_linear_1d(z_thref(1:size(z_thref)), f_thref(1:size(f_thref)), zgrid, &
238  current_state%global_grid%configuration%vertical%thref)
239  deallocate(z_thref, f_thref)
240  end if
241  else
242  current_state%global_grid%configuration%vertical%thref(:)=current_state%thref0
243  end if
244 
245  if (l_init_pl_theta)then
246  allocate(z_init_pl_theta(options_get_array_size(current_state%options_database, "z_init_pl_theta")), &
247  f_init_pl_theta(options_get_array_size(current_state%options_database, "f_init_pl_theta")))
248  call options_get_real_array(current_state%options_database, "z_init_pl_theta", z_init_pl_theta)
249  call options_get_real_array(current_state%options_database, "f_init_pl_theta", f_init_pl_theta)
250  call check_top(zztop, z_init_pl_theta(size(z_init_pl_theta)), 'z_init_pl_theta')
251  call check_input_levels(size(z_init_pl_theta), size(f_init_pl_theta), "f_init_pl_theta")
252  zgrid=current_state%global_grid%configuration%vertical%zn(:)
253  call piecewise_linear_1d(z_init_pl_theta(1:size(z_init_pl_theta)), f_init_pl_theta(1:size(f_init_pl_theta)), zgrid, &
254  current_state%global_grid%configuration%vertical%theta_init)
255  if (l_matchthref) then
256  if(.not. current_state%use_anelastic_equations) then
257  call log_master_log(log_error, "Non-anelastic equation set and l_maththref are incompatible")
258  end if
259  current_state%global_grid%configuration%vertical%thref = current_state%global_grid%configuration%vertical%theta_init
260  end if
261  if (.not. current_state%continuation_run) then
262  do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
263  do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
264  current_state%th%data(:,j,i) = current_state%global_grid%configuration%vertical%theta_init(:) - &
265  current_state%global_grid%configuration%vertical%thref(:)
266  end do
267  end do
268  end if
269  deallocate(z_init_pl_theta, f_init_pl_theta)
270  end if
271 
272  if (l_init_pl_u)then
273  allocate(z_init_pl_u(options_get_array_size(current_state%options_database, "z_init_pl_u")), &
274  f_init_pl_u(options_get_array_size(current_state%options_database, "f_init_pl_u")))
275  call options_get_real_array(current_state%options_database, "z_init_pl_u", z_init_pl_u)
276  call options_get_real_array(current_state%options_database, "f_init_pl_u", f_init_pl_u)
277  call check_top(zztop, z_init_pl_u(size(z_init_pl_u)), 'z_init_pl_u')
278  call check_input_levels(size(z_init_pl_u), size(f_init_pl_u), "f_init_pl_u")
279  zgrid=current_state%global_grid%configuration%vertical%zn(:)
280  call piecewise_linear_1d(z_init_pl_u(1:size(z_init_pl_u)), f_init_pl_u(1:size(f_init_pl_u)), &
281  zgrid, current_state%global_grid%configuration%vertical%u_init)
282  if (.not. current_state%continuation_run) then
283  do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
284  do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
285  current_state%u%data(:,j,i) = current_state%global_grid%configuration%vertical%u_init(:)
286  end do
287  end do
288  end if
289  deallocate(z_init_pl_u, f_init_pl_u)
290  end if
291 
292  if (l_init_pl_v)then
293  allocate(z_init_pl_v(options_get_array_size(current_state%options_database, "z_init_pl_v")), &
294  f_init_pl_v(options_get_array_size(current_state%options_database, "f_init_pl_v")))
295  call options_get_real_array(current_state%options_database, "z_init_pl_v", z_init_pl_v)
296  call options_get_real_array(current_state%options_database, "f_init_pl_v", f_init_pl_v)
297  call check_top(zztop, z_init_pl_v(size(z_init_pl_v)), 'z_init_pl_v')
298  call check_input_levels(size(z_init_pl_v), size(f_init_pl_v), "f_init_pl_v")
299  zgrid=current_state%global_grid%configuration%vertical%zn(:)
300  call piecewise_linear_1d(z_init_pl_v(1:size(z_init_pl_v)), f_init_pl_v(1:size(f_init_pl_v)), &
301  zgrid, current_state%global_grid%configuration%vertical%v_init)
302  if (.not. current_state%continuation_run) then
303  do i=current_state%local_grid%local_domain_start_index(x_index), current_state%local_grid%local_domain_end_index(x_index)
304  do j=current_state%local_grid%local_domain_start_index(y_index), current_state%local_grid%local_domain_end_index(y_index)
305  current_state%v%data(:,j,i) = current_state%global_grid%configuration%vertical%v_init(:)
306  end do
307  end do
308  end if
309  deallocate(z_init_pl_v, f_init_pl_v)
310  end if
311 
312  if (l_init_pl_q)then
313  nq_init=size(names_init_pl_q)
314  allocate(z_init_pl_q(options_get_array_size(current_state%options_database, "z_init_pl_q")))
315  call options_get_real_array(current_state%options_database, "z_init_pl_q", z_init_pl_q)
316  nzq=size(z_init_pl_q)
317  call check_top(zztop, z_init_pl_q(nzq), 'z_init_pl_q')
318  zgrid=current_state%global_grid%configuration%vertical%zn(:)
319  allocate(f_init_pl_q_tmp(nq_init*nzq))
320  call options_get_real_array(current_state%options_database, "f_init_pl_q", f_init_pl_q_tmp)
321  !call check_input_levels(size(z_init_pl_q), size(f_init_pl_q_tmp), "f_init_pl_q_tmp")
322  allocate(f_init_pl_q(nzq, nq_init))
323  f_init_pl_q(1:nzq, 1:nq_init)=reshape(f_init_pl_q_tmp, (/nzq, nq_init/))
324  do n=1, nq_init
325  iq=get_q_index(trim(names_init_pl_q(n)), 'piecewise_initialization')
326  call check_input_levels(size(z_init_pl_q), size(f_init_pl_q(1:nzq,n)), "f_init_pl_q")
327  call piecewise_linear_1d(z_init_pl_q(1:nzq), f_init_pl_q(1:nzq,n), zgrid, &
328  current_state%global_grid%configuration%vertical%q_init(:,iq))
329  if (.not. current_state%continuation_run) then
330  do i=current_state%local_grid%local_domain_start_index(x_index), &
331  current_state%local_grid%local_domain_end_index(x_index)
332  do j=current_state%local_grid%local_domain_start_index(y_index), &
333  current_state%local_grid%local_domain_end_index(y_index)
334  current_state%q(iq)%data(:,j,i) = current_state%global_grid%configuration%vertical%q_init(:, iq)
335  end do
336  end do
337  end if
338  end do
339  deallocate(f_init_pl_q_tmp, z_init_pl_q, f_init_pl_q, names_init_pl_q)
340  end if
341  deallocate(zgrid)
342  end subroutine calculate_initial_profiles
343 
348  subroutine calculate_mixing_length_for_neutral_case(current_state, vertical_grid, kkp)
349  type(model_state_type), intent(inout) :: current_state
350  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
351  integer, intent(in) :: kkp
352 
353  integer :: k
354 
355  do k=2, kkp-1
356  vertical_grid%rneutml(k)=sqrt(1.0_default_precision/(1.0_default_precision/(von_karman_constant*&
357  (vertical_grid%z(k)+z0))**2+1.0_default_precision/current_state%rmlmax**2) )
358  vertical_grid%rneutml_sq(k)=vertical_grid%rneutml(k)*vertical_grid%rneutml(k)
359  end do
360  end subroutine calculate_mixing_length_for_neutral_case
361 
366  subroutine set_buoyancy_coefficient(current_state, vertical_grid, kkp)
367  type(model_state_type), intent(inout), target :: current_state
368  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
369  integer, intent(in) :: kkp
370 
371  integer :: k
372 
373  if (.not. current_state%passive_th) then
374  if(current_state%use_anelastic_equations)then
375  do k=1, kkp-1
376  vertical_grid%buoy_co(k)=cp*(vertical_grid%prefn(k)**r_over_cp-vertical_grid%prefn(k+1)**r_over_cp)/&
377  ((current_state%surface_reference_pressure**r_over_cp)*vertical_grid%dzn(k+1))
378  end do
379  else
380  vertical_grid%buoy_co(1:kkp-1)=g/current_state%thref0 ! _Boussinesq
381  end if
382  ! Dummy value at top level
383  vertical_grid%buoy_co(kkp)=0.
384  else
385  vertical_grid%buoy_co(:)=0.
386  end if
387  end subroutine set_buoyancy_coefficient
388 
393  subroutine setup_reference_state_liquid_water_temperature_and_saturation(current_state, vertical_grid, kkp)
394  type(model_state_type), intent(inout) :: current_state
395  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
396  integer, intent(in) :: kkp
397 
398  real(kind=default_precision) :: delta_t=1.0_default_precision, qlinit, tinit, qsatin, dqsatdtin, dsatfacin
399  integer :: iter, k
400 
401  do k=1,kkp
402  vertical_grid%tref(k)=vertical_grid%thref(k)*(vertical_grid%prefn(k)/current_state%surface_reference_pressure)**r_over_cp
403  vertical_grid%tstarpr(k)=0.0_default_precision
404  end do
405  if (current_state%th%active) then
406  ! PREFRCP is used and hence calculated if theta is active
407  do k = 1,kkp
408  vertical_grid%prefrcp(k)=(current_state%surface_reference_pressure/vertical_grid%prefn(k))**r_over_cp
409  vertical_grid%rprefrcp(k)=1.0_default_precision/vertical_grid%prefrcp(k)
410  ! Denotion between setup run and chain run in LEM - need to consider here too
411  vertical_grid%qsat(k)=qsaturation(vertical_grid%tref(k), 0.01_default_precision*vertical_grid%prefn(k))
412  vertical_grid%dqsatdt(k)=(qsaturation(vertical_grid%tref(k)+delta_t, 0.01_default_precision*vertical_grid%prefn(k)) -&
413  qsaturation(vertical_grid%tref(k)-delta_t, 0.01_default_precision*vertical_grid%prefn(k)))/&
414  (2.0_default_precision*delta_t)
415  vertical_grid%qsatfac(k)=1.0_default_precision/(1.0_default_precision+rlvap_over_cp*vertical_grid%dqsatdt(k))
416  end do
417  if (current_state%calculate_th_and_q_init) then
418  do k=1,kkp
419  ! !Note that at this point THETA_INIT and QINIT(IQ=1) are still
420  ! !theta_l and q_t, as read in from configuration.
421  ! ! start from input QL profile
422  qlinit=qinit(k, current_state%liquid_water_mixing_ratio_index)
423  do iter=1,5
424  ! ! calculate T and thence new q_l from Taylor expansion
425  ! ! keeping theta_l and q_t fixed
426  ! ! Note theta_l = theta - (L/c_p)*q_l here
427  tinit = vertical_grid%theta_init(k)*vertical_grid%rprefrcp(k) + rlvap_over_cp*qlinit
428  qsatin = qsaturation(tinit, 0.01_default_precision*vertical_grid%prefn(k))
429  dqsatdtin = dqwsatdt(qsatin, tinit)
430  dsatfacin=( 1.0_default_precision/(1.0_default_precision + rlvap_over_cp*dqsatdtin*vertical_grid%rprefrcp(k)))
431  qlinit=max(0.0_default_precision, (qinit(k, current_state%water_vapour_mixing_ratio_index)-&
432  (qsatin+dqsatdtin*(vertical_grid%theta_init(k)*vertical_grid%rprefrcp(k)-tinit) ))*dsatfacin)
433  end do
434  qinit(k, current_state%liquid_water_mixing_ratio_index)=qlinit
435  qinit(k, current_state%water_vapour_mixing_ratio_index)=qinit(k,current_state%water_vapour_mixing_ratio_index)-qlinit
436  vertical_grid%theta_init(k)=vertical_grid%theta_init(k)+rlvap_over_cp*qlinit
437 
438  ! Denotion between setup run and chain run in LEM - need to consider here too
439  vertical_grid%tstarpr(k)= tinit-vertical_grid%tref(k)
440  vertical_grid%qsat(k)=qsatin
441  vertical_grid%dqsatdt(k)=dqsatdtin
442  vertical_grid%qsatfac(k)= ( 1.0_default_precision/ ( 1.0_default_precision + rlvap_over_cp*dqsatdtin ) )
443  end do
444  endif
445  end if
446  end subroutine setup_reference_state_liquid_water_temperature_and_saturation
447 
452  subroutine set_up_vertical_reference_properties(current_state, vertical_grid, kkp)
453  type(model_state_type), intent(inout) :: current_state
454  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
455  integer, intent(in) :: kkp
456 
457  integer :: k
458 
459  do k=1,kkp
460 ! vertical_grid%thref(k)=current_state%thref0
461 ! vertical_grid%theta_init(k)=vertical_grid%thref(k) ! In LEM this can also be set from configuration (TODO)
462  vertical_grid%prefn(k)=0.0_default_precision
463  vertical_grid%pdiff(k)=0.0_default_precision
464  vertical_grid%rho(k)=current_state%rhobous
465  vertical_grid%rhon(k)=current_state%rhobous
466  end do
467  do k=1,kkp-1
468  vertical_grid%dthref(k)=vertical_grid%thref(k+1)-vertical_grid%thref(k)
469  end do
470  vertical_grid%dthref(kkp)=0.0_default_precision
471  end subroutine set_up_vertical_reference_properties
472 
481  subroutine set_up_and_smooth_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop, nsmth, origional_setup, continuation_run)
482  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
483  integer, dimension(:), intent(in) :: kgd
484  real(kind=default_precision), dimension(:), intent(in) :: hgd
485  integer, intent(in) :: ninitp, kkp, nsmth
486  real(kind=default_precision),intent(in) :: zztop
487  logical, intent(in) :: origional_setup, continuation_run
488 
489  integer :: k
490 
491  if (.not. continuation_run) then
492  if (origional_setup) then
493  call original_vertical_grid_setup(vertical_grid, kgd, hgd, ninitp, kkp, zztop, nsmth)
494  else
495  call new_vertical_grid_setup(vertical_grid, kgd, kkp, zztop)
496  end if
497  end if
498 
499  ! Regardless of the vertical grid computation method, set the level deltas
500  do k=2,kkp
501  vertical_grid%dz(k)=vertical_grid%z(k)-vertical_grid%z(k-1)
502  vertical_grid%dzn(k)= vertical_grid%zn(k)-vertical_grid%zn(k-1)
503  vertical_grid%rdz(k)=1./vertical_grid%dz(k)
504  vertical_grid%rdzn(k)=1./vertical_grid%dzn(k)
505  end do
506  vertical_grid%dzn(1)=0.d0
507  end subroutine set_up_and_smooth_grid
508 
516  subroutine original_vertical_grid_setup(vertical_grid, kgd, hgd, ninitp, kkp, zztop, nsmth)
517  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
518  integer, dimension(:), intent(in) :: kgd
519  real(kind=default_precision), dimension(:), intent(in) :: hgd
520  integer, intent(in) :: ninitp, kkp, nsmth
521  real(kind=default_precision), intent(in) :: zztop
522 
523  integer :: n, k
524 
525  call create_linear_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop)
526  ! Smooth grid
527  vertical_grid%z(1)=0.0_default_precision
528  vertical_grid%z(kkp)=zztop
529  do n=1,nsmth
530  do k=2,kkp
531  vertical_grid%zn(k)=0.5_default_precision*(vertical_grid%z(k)+vertical_grid%z(k-1))
532  end do
533  do k=2,kkp-1
534  vertical_grid%z(k)=0.5_default_precision*(vertical_grid%zn(k)+vertical_grid%zn(k+1))
535  end do
536  end do
537  ! Fourth order interpolation
538  do k=3,kkp-1
539  vertical_grid%zn(k)=0.0625_default_precision*(9.0_default_precision*&
540  (vertical_grid%z(k-1)+vertical_grid%z(k))-vertical_grid%z(k+1)-vertical_grid%z(k-2))
541  end do
542  vertical_grid%zn(2)=0.5_default_precision*(vertical_grid%z(1)+vertical_grid%z(2))
543  vertical_grid%zn(1)=-vertical_grid%zn(2)
544  vertical_grid%zn(kkp)=0.5_default_precision*(vertical_grid%z(kkp-1)+vertical_grid%z(kkp))
545  end subroutine original_vertical_grid_setup
546 
552  subroutine new_vertical_grid_setup(vertical_grid, kgd, kkp, zztop)
553  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
554  integer, dimension(:), intent(in) :: kgd
555  integer, intent(in) :: kkp
556  real(kind=default_precision), intent(in) :: zztop
557 
558  real(kind=default_precision) :: a(2*kkp), r1, d1, dd, d0, a0
559  logical :: first_gt=.true.
560  integer :: k, k0
561 
562  r1=1.10_default_precision
563  d1=10.0_default_precision
564 
565  dd=0.5_default_precision*d1
566  d0=dd
567  a(1)=-dd
568  a(2)=0.0_default_precision
569  do k=3, kkp*2
570  if (d0 .gt. dd .or. k==3) then
571  a(k)=a(k-1)+dd
572  if (.not. (dd .gt. 25.0_default_precision .and. a(k) .lt. 2000.0_default_precision)) then
573  if (a(k) .lt. 2000.0_default_precision) then
574  dd=dd*r1
575  else
576  dd=dd*(1.0_default_precision+(r1-1.0_default_precision)/1.5_default_precision)
577  end if
578  end if
579  d0=(zztop-a(k))/real(kkp*2-k, kind=default_precision)
580  else
581  if (first_gt) then
582  k0=k
583  a0=a(k-1)+d0
584  first_gt=.false.
585  end if
586  a(k)=a0+real(k-k0, kind=default_precision)*d0
587  end if
588  end do
589 
590  do k=1, kkp
591  vertical_grid%z(k)=a(k*2)
592  vertical_grid%zn(k)=a(2*k-1)
593  end do
594  end subroutine new_vertical_grid_setup
595 
606  subroutine create_linear_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop)
607  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
608  integer, dimension(:), intent(in) :: kgd
609  real(kind=default_precision), dimension(:), intent(in) :: hgd
610  integer, intent(in) :: ninitp, kkp
611  real(kind=default_precision), intent(in) :: zztop
612 
613  integer :: kmax, k, i
614  real(kind=default_precision) :: zmax
615 
616  vertical_grid%z(1) = 0.0_default_precision
617  kmax=kgd(1)
618  zmax=0.0_default_precision
619  if (kgd(1) .gt. 1) then
620  do k=1,kgd(1)
621  ! Loop up to first division point
622  vertical_grid%z(k)=real(k-1, kind=default_precision)*hgd(1)/real(kgd(1)-1, kind=default_precision)
623  end do
624  do i=2,ninitp
625  if(kgd(i) .gt. 0) then
626  kmax=kgd(i)
627  zmax=hgd(i)
628 
629  do k=kgd(i-1)+1,kgd(i)
630  vertical_grid%z(k)=hgd(i-1)+(hgd(i)-hgd(i-1))*real(k-kgd(i-1), kind=default_precision)&
631  /real(kgd(i)-kgd(i-1), kind=default_precision)
632  end do
633  end if
634  end do
635  end if
636  if(kmax .lt. kkp)then
637  do k=kmax,kkp
638  ! Handle any points above the kth max division
639  vertical_grid%z(k)=zmax+(zztop-zmax)*real(k-kmax, kind=default_precision)/real(kkp-kmax, kind=default_precision)
640  end do
641  end if
642  end subroutine create_linear_grid
643 
648  subroutine allocate_vertical_grid_data(vertical_grid, n, nq)
649  type(vertical_grid_configuration_type), intent(inout) :: vertical_grid
650  integer, intent(in) :: n
651  integer, intent(in) :: nq
652 
653  allocate(vertical_grid%dz(n), vertical_grid%dzn(n),&
654  vertical_grid%czb(n), vertical_grid%cza(n), vertical_grid%czg(n), vertical_grid%czh(n),&
655  vertical_grid%rdz(n), vertical_grid%rdzn(n), vertical_grid%tzc1(n), vertical_grid%tzc2(n),&
656  vertical_grid%tzd1(n), vertical_grid%tzd2(n), vertical_grid%theta_init(n), vertical_grid%temp_init(n), &
657  vertical_grid%rh_init(n), &
658  vertical_grid%tref(n), vertical_grid%prefn(n), vertical_grid%pdiff(n), vertical_grid%prefrcp(n), &
659  vertical_grid%rprefrcp(n), vertical_grid%rho(n), vertical_grid%rhon(n), vertical_grid%tstarpr(n), &
660  vertical_grid%qsat(n), vertical_grid%dqsatdt(n), vertical_grid%qsatfac(n), vertical_grid%dthref(n), &
661  vertical_grid%rneutml(n), vertical_grid%rneutml_sq(n), vertical_grid%buoy_co(n), &
662  vertical_grid%u_init(n), vertical_grid%v_init(n), vertical_grid%theta_rand(n), vertical_grid%w_rand(n), &
663  vertical_grid%w_subs(n), vertical_grid%u_force(n), vertical_grid%v_force(n), vertical_grid%theta_force(n))
664 
665  if (.not. allocated(vertical_grid%thref)) allocate(vertical_grid%thref(n))
666  if (.not. allocated(vertical_grid%z)) allocate(vertical_grid%z(n))
667  if (.not. allocated(vertical_grid%zn)) allocate(vertical_grid%zn(n))
668 
669  allocate(vertical_grid%q_rand(n,nq), vertical_grid%q_init(n,nq), vertical_grid%q_force(n,nq))
670  end subroutine allocate_vertical_grid_data
671 
674  subroutine initialise_horizontalgrid_configuration_types(current_state)
675  type(model_state_type), intent(inout) :: current_state
676 
677  current_state%global_grid%configuration%horizontal%dx = merge(real(current_state%global_grid%resolution(x_index)), &
678  default_spacing, current_state%global_grid%active(x_index))
679  current_state%global_grid%configuration%horizontal%dy = merge(real(current_state%global_grid%resolution(y_index)), &
680  default_spacing, current_state%global_grid%active(y_index))
681 
682  current_state%global_grid%configuration%horizontal%cx=1./current_state%global_grid%configuration%horizontal%dx
683  current_state%global_grid%configuration%horizontal%cy=1./current_state%global_grid%configuration%horizontal%dy
684  current_state%global_grid%configuration%horizontal%cx2=current_state%global_grid%configuration%horizontal%cx ** 2
685  current_state%global_grid%configuration%horizontal%cy2=current_state%global_grid%configuration%horizontal%cy ** 2
686  current_state%global_grid%configuration%horizontal%cxy=current_state%global_grid%configuration%horizontal%cx * &
687  current_state%global_grid%configuration%horizontal%cy
688  current_state%global_grid%configuration%horizontal%tcx=&
689  0.25_default_precision/current_state%global_grid%configuration%horizontal%dx
690  current_state%global_grid%configuration%horizontal%tcy=&
691  0.25_default_precision/current_state%global_grid%configuration%horizontal%dy
692  end subroutine initialise_horizontalgrid_configuration_types
693 
699  subroutine set_anelastic_pressure(current_state)
700  type(model_state_type), intent(inout) :: current_state
701 
702  if (current_state%use_anelastic_equations) then
703  call compute_anelastic_pressure_profile_and_density(current_state)
704  else
705  if (current_state%passive_th) then
706  current_state%global_grid%configuration%vertical%prefn=0.0_default_precision
707  else
708  call compute_anelastic_pressure_profile_only(current_state)
709  end if
710  current_state%global_grid%configuration%vertical%rho=current_state%rhobous
711  current_state%global_grid%configuration%vertical%rhon=current_state%rhobous
712  current_state%global_grid%configuration%vertical%pdiff=0.0_default_precision
713  end if
714  end subroutine set_anelastic_pressure
715 
718  subroutine compute_anelastic_pressure_profile_only(current_state)
719  type(model_state_type), intent(inout) :: current_state
720 
721  integer :: ipass, k
722  real(kind=default_precision) :: p0 &!pressure at z=0 adjustments made after 1st iteration so P0=PSF after 2nd iteration
723  , ptop &!pressure at z=ZN(KKP)
724  , thprof(current_state%local_grid%size(z_index))
725 
726 
727  ! TODO: NOTE - we are mocking in thprof at the moment, this should be read from a configuration and used here instead
728  thprof=0.0_default_precision
729  ptop=0.0_default_precision
730  current_state%global_grid%configuration%vertical%pdiff(current_state%local_grid%size(z_index))=0.0_default_precision
731 
732  do ipass=1,2 ! _after first pass adjust PTOP
733  current_state%global_grid%configuration%vertical%prefn(current_state%local_grid%size(z_index))=&
734  (ptop/current_state%surface_reference_pressure)**r_over_cp
735  do k=current_state%local_grid%size(z_index)-1,1,-1
736  current_state%global_grid%configuration%vertical%pdiff(k)=g*&
737  current_state%global_grid%configuration%vertical%dzn(k+1)/(0.5_default_precision*cp*&
738  (current_state%global_grid%configuration%vertical%thref(k)+&
739  current_state%global_grid%configuration%vertical%thref(k+1)))
740  end do
741  do k=current_state%local_grid%size(z_index)-1,1,-1
742  current_state%global_grid%configuration%vertical%prefn(k)=&
743  current_state%global_grid%configuration%vertical%prefn(k+1)+&
744  current_state%global_grid%configuration%vertical%pdiff(k)
745  end do
746  do k=current_state%local_grid%size(z_index),1,-1
747  current_state%global_grid%configuration%vertical%prefn(k)=current_state%surface_reference_pressure*&
748  current_state%global_grid%configuration%vertical%prefn(k)**(1.0_default_precision/r_over_cp)
749  end do
750  p0=0.5_default_precision*(current_state%global_grid%configuration%vertical%prefn(1)+&
751  current_state%global_grid%configuration%vertical%prefn(2))
752  if (ipass .eq. 1) then
753  ptop=current_state%surface_pressure**r_over_cp+ptop**r_over_cp-p0**r_over_cp
754  if (ptop .le. 0.0_default_precision .and. current_state%parallel%my_rank==0) then
755  call log_log(log_error, "Negative ptop in setup of anelastic. Need a warmer THREF or different setup options")
756  end if
757  ptop=ptop**(1.0_default_precision/r_over_cp)
758  end if
759  end do
760  end subroutine compute_anelastic_pressure_profile_only
761 
764  subroutine compute_anelastic_pressure_profile_and_density(current_state)
765  type(model_state_type), intent(inout) :: current_state
766 
767  integer :: ipass, k
768  real(kind=default_precision) :: p0 &!pressure at z=0 adjustments made after 1st iteration so P0=PSF after 2nd iteration
769  , ptop &!pressure at z=ZN(KKP)
770  , thfactor !factor for multiplying TH profile (if IADJANELP=2)
771 
772  ptop=0.0_default_precision
773  current_state%global_grid%configuration%vertical%pdiff(current_state%local_grid%size(z_index))=0.0_default_precision
774  do ipass=1,2 ! _after first pass, may adjust basic states
775  if (ipass .eq. 1 .or. anelastic_profile_mode .gt. 1) then
776  current_state%global_grid%configuration%vertical%prefn(current_state%local_grid%size(z_index))=&
777  (ptop/current_state%surface_reference_pressure)**r_over_cp
778  do k=current_state%local_grid%size(z_index)-1,1,-1
779  current_state%global_grid%configuration%vertical%pdiff(k)=g*&
780  current_state%global_grid%configuration%vertical%dzn(k+1)/(0.5_default_precision*cp*&
781  (current_state%global_grid%configuration%vertical%thref(k)+&
782  current_state%global_grid%configuration%vertical%thref(k+1)))
783  end do
784  do k=current_state%local_grid%size(z_index)-1,1,-1
785  current_state%global_grid%configuration%vertical%prefn(k)=&
786  current_state%global_grid%configuration%vertical%prefn(k+1)+&
787  current_state%global_grid%configuration%vertical%pdiff(k)
788  end do
789  do k=current_state%local_grid%size(z_index),1,-1
790  current_state%global_grid%configuration%vertical%prefn(k)=current_state%surface_reference_pressure*&
791  current_state%global_grid%configuration%vertical%prefn(k)**(1.0_default_precision/r_over_cp)
792  end do
793  end if
794  p0=0.5_default_precision*(current_state%global_grid%configuration%vertical%prefn(1)+&
795  current_state%global_grid%configuration%vertical%prefn(2))
796  ! !-------------------------------------------------------------
797  ! ! _If IADJANELP>1 we adjust the basic states to ensure P0=PSF,
798  ! ! as follows:
799  ! ! IADJANELP=2 adjust THREF profile by constant factor
800  ! ! IADJANELP=3 adjust PSF
801  ! ! IADJANELP=4 adjust PTOP
802  ! ! _Option 3 tends to give rather large changes in PSF, so
803  ! ! I prefer 2 or 4 for most purposes
804  ! !-------------------------------------------------------------
805  if (ipass .eq. 1 .and. anelastic_profile_mode .eq. 2) then
806  ! ! _adjust THREF profile by constant factor to enforce
807  ! ! P0 = (fixed) PSF
808  thfactor=((p0/current_state%surface_reference_pressure)**r_over_cp-(ptop/current_state%surface_reference_pressure)**&
809  r_over_cp)/((current_state%surface_pressure/current_state%surface_reference_pressure)**r_over_cp-&
810  (ptop/current_state%surface_reference_pressure)**r_over_cp)
811  do k=1,current_state%local_grid%size(z_index)
812  current_state%global_grid%configuration%vertical%thref(k)=&
813  current_state%global_grid%configuration%vertical%thref(k)*thfactor
814  end do
815  end if
816  if (ipass .eq. 1 .and. anelastic_profile_mode .eq. 4) then
817  ! ! _adjust PTOP so that P0 = (fixed) PSF
818  ptop=current_state%surface_pressure**r_over_cp+ptop**r_over_cp-p0**r_over_cp
819  if (ptop .le. 0.0_default_precision .and. current_state%parallel%my_rank==0) then
820  call log_log(log_error, "Negative ptop in setup of anelastic. Need a warmer THREF or different setup options")
821  end if
822  ptop=ptop**(1.0_default_precision/r_over_cp)
823  end if
824  end do
825  ! !---------------------------------------
826  ! ! _Finally compute density from pressure
827  ! !---------------------------------------
828  do k=1,current_state%local_grid%size(z_index)
829  current_state%global_grid%configuration%vertical%rhon(k)=current_state%global_grid%configuration%vertical%prefn(k)&
830  /(r*current_state%global_grid%configuration%vertical%thref(k)*&
831  (current_state%global_grid%configuration%vertical%prefn(k)/current_state%surface_reference_pressure)**r_over_cp)
832  end do
833  do k=1,current_state%local_grid%size(z_index)-1
834  current_state%global_grid%configuration%vertical%rho(k)=sqrt(current_state%global_grid%configuration%vertical%rhon(k)*&
835  current_state%global_grid%configuration%vertical%rhon(k+1))
836  end do
837  current_state%global_grid%configuration%vertical%rho(current_state%local_grid%size(z_index))=&
838  current_state%global_grid%configuration%vertical%rhon(current_state%local_grid%size(z_index))**2/&
839  current_state%global_grid%configuration%vertical%rhon(current_state%local_grid%size(z_index)-1)
840  end subroutine compute_anelastic_pressure_profile_and_density
841 
842 
843  subroutine check_top(zztop, z, info)
844  real(kind=default_precision), intent(in) :: zztop
845  real(kind=default_precision), intent(in) :: z
846  character(*), intent(in) :: info
847 
848  if (z<zztop)then
849  call log_master_log(log_error, "Top of input profile is below the top of the domain:"//trim(info))
850  end if
851 
852  end subroutine check_top
853 
854  subroutine check_input_levels(z_levels, field_levels, field)
855  integer, intent(in) :: z_levels
856  integer, intent(in) :: field_levels
857  character(*), intent(in) :: field
858 
859  if (z_levels /= field_levels)then
860  call log_master_log(log_error, "Input levels not equal for "//trim(field)//", z_levels = "// &
861  trim(conv_to_string(z_levels))//" field_levels = "//conv_to_string(field_levels))
862  end if
863 
864  end subroutine check_input_levels
865 
866  subroutine set_qv_init_from_rh(current_state)
867 
868  type(model_state_type), intent(inout) :: current_state
869 
870  logical :: l_init_pl_rh ! if .true. then initialize relative humidity field
871  real(kind=default_precision) :: zztop ! top of the domain
872  real(kind=default_precision) :: qsat
873  real(kind=default_precision), dimension(:), allocatable :: f_init_pl_rh ! Initial node values for relative humidity variable
874  real(kind=default_precision), dimension(:), allocatable :: z_init_pl_rh ! Initial node height values for relative humidity variable
875  real(kind=default_precision), allocatable :: zgrid(:) ! z grid to use in interpolation
876  real(kind=default_precision), allocatable :: tdegk(:) ! temperature in Kelvin
877  integer :: i,j,n, k ! loop counters
878  integer :: iq ! temporary q varible index
879 
880  type(vertical_grid_configuration_type) :: vertical_grid
881 
882  vertical_grid=current_state%global_grid%configuration%vertical
883 
884  allocate(zgrid(current_state%local_grid%local_domain_end_index(z_index)))
885 
886  zztop = current_state%global_grid%top(z_index)
887 
888  l_init_pl_rh=options_get_logical(current_state%options_database, "l_init_pl_rh")
889 
890  if (l_init_pl_rh)then
891  allocate(z_init_pl_rh(options_get_array_size(current_state%options_database, "z_init_pl_rh")), &
892  f_init_pl_rh(options_get_array_size(current_state%options_database, "f_init_pl_rh")))
893  call options_get_real_array(current_state%options_database, "z_init_pl_rh", z_init_pl_rh)
894  call options_get_real_array(current_state%options_database, "f_init_pl_rh", f_init_pl_rh)
895  call check_top(zztop, z_init_pl_rh(size(z_init_pl_rh)), 'z_init_pl_rh')
896  call check_input_levels(size(z_init_pl_rh), size(f_init_pl_rh), "f_init_pl_rh")
897  zgrid=current_state%global_grid%configuration%vertical%zn(:)
898  call piecewise_linear_1d(z_init_pl_rh(1:size(z_init_pl_rh)), f_init_pl_rh(1:size(f_init_pl_rh)), zgrid, &
899  current_state%global_grid%configuration%vertical%rh_init)
900 
901  if (.not. current_state%passive_q .and. current_state%th%active) then
902  iq=get_q_index('vapour', 'piecewise_initialization')
903  allocate(tdegk(current_state%local_grid%local_domain_end_index(z_index)))
904  tdegk(:) = current_state%global_grid%configuration%vertical%theta_init(:)* &
905  (vertical_grid%prefn(:)/current_state%surface_reference_pressure)**r_over_cp
906  do k = current_state%local_grid%local_domain_start_index(z_index), &
907  current_state%local_grid%local_domain_end_index(z_index)
908  qsat=qsaturation(tdegk(k), current_state%global_grid%configuration%vertical%prefn(k)/100.)
909  current_state%global_grid%configuration%vertical%q_init(k, iq) = &
910  (current_state%global_grid%configuration%vertical%rh_init(k)/100.0)*qsat
911  !print *, current_state%global_grid%configuration%vertical%rh_init(k), &
912  ! current_state%global_grid%configuration%vertical%q_init(k, iq), &
913  ! TdegK(k)
914  enddo
915  if (.not. current_state%continuation_run) then
916  do i=current_state%local_grid%local_domain_start_index(x_index), &
917  current_state%local_grid%local_domain_end_index(x_index)
918  do j=current_state%local_grid%local_domain_start_index(y_index), &
919  current_state%local_grid%local_domain_end_index(y_index)
920  current_state%q(iq)%data(:,j,i) = current_state%global_grid%configuration%vertical%q_init(:, iq)
921  end do
922  end do
923  end if
924 
925  deallocate(tdegk)
926  else
927  call log_master_log(log_error, "Initialising with RH but q and/or theta passive")
928  end if
929 
930  deallocate(z_init_pl_rh, f_init_pl_rh)
931  end if
932 
933  end subroutine set_qv_init_from_rh
934 
935 end module gridmanager_mod
logging_mod::log_error
integer, parameter, public log_error
Only log ERROR messages.
Definition: logging.F90:11
conversions_mod
Conversion between common inbuilt FORTRAN data types.
Definition: conversions.F90:5
science_constants_mod::r_over_cp
real(kind=default_precision), public r_over_cp
Definition: scienceconstants.F90:13
saturation_mod
Saturation physics functionality which is used throughout the code.
Definition: saturation.F90:2
gridmanager_mod::calculate_mixing_length_for_neutral_case
subroutine calculate_mixing_length_for_neutral_case(current_state, vertical_grid, kkp)
Calculates the mixing length for the neutral case.
Definition: gridmanager.F90:349
optionsdatabase_mod::options_get_array_size
integer function, public options_get_array_size(options_database, key)
Gets the size of the array held in the options database corresponding to a specific key.
Definition: optionsdatabase.F90:342
gridmanager_mod::qinit
real(kind=default_precision), dimension(:,:), allocatable qinit
Definition: gridmanager.F90:29
gridmanager_mod::compute_anelastic_pressure_profile_and_density
subroutine compute_anelastic_pressure_profile_and_density(current_state)
Computes the anelastic pressure and density - if we are using Anelastic approximation.
Definition: gridmanager.F90:765
interpolation_mod::piecewise_linear_1d
subroutine piecewise_linear_1d(zvals, vals, zgrid, field)
Does a simple 1d piecewise linear interpolation.
Definition: interpolation.F90:16
grids_mod::x_index
integer, parameter, public x_index
Definition: grids.F90:14
saturation_mod::qsaturation
real(kind=default_precision) function, public qsaturation(temperature, pressure)
Function to return the saturation mixing ratio over water based on tetans formular QS=3....
Definition: saturation.F90:22
optionsdatabase_mod::options_get_integer
integer function, public options_get_integer(options_database, key, index)
Retrieves an integer value from the database that matches the provided key.
Definition: optionsdatabase.F90:217
grids_mod::y_index
integer, parameter, public y_index
Definition: grids.F90:14
logging_mod::log_info
integer, parameter, public log_info
Log INFO, WARNING and ERROR messages.
Definition: logging.F90:13
gridmanager_mod::create_linear_grid
subroutine create_linear_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop)
Creates the linear vertical grid based upon the configuration properties.
Definition: gridmanager.F90:607
gridmanager_mod::check_top
subroutine check_top(zztop, z, info)
Definition: gridmanager.F90:844
logging_mod::log_log
subroutine, public log_log(level, message, str)
Logs a message at the specified level. If the level is above the current level then the message is ig...
Definition: logging.F90:75
science_constants_mod::cp
real(kind=default_precision), public cp
Definition: scienceconstants.F90:13
optionsdatabase_mod::options_get_logical_array
subroutine, public options_get_logical_array(options_database, key, array_data, from, to)
Retrieves an entire (or subset) logical array.
Definition: optionsdatabase.F90:176
gridmanager_mod::initialise_verticalgrid_configuration_type
subroutine initialise_verticalgrid_configuration_type(current_state)
Will initialise the vertical grid configuration.
Definition: gridmanager.F90:92
gridmanager_mod::new_vertical_grid_setup
subroutine new_vertical_grid_setup(vertical_grid, kgd, kkp, zztop)
The newer vertical grid setup routine.
Definition: gridmanager.F90:553
monc_component_mod
Interfaces and types that MONC components must specify.
Definition: monc_component.F90:6
gridmanager_mod::finalise_callback
subroutine finalise_callback(current_state)
Called as MONC exits, for good practice this will deallocate the memory used for the grids.
Definition: gridmanager.F90:70
optionsdatabase_mod::options_get_string
character(len=string_length) function, public options_get_string(options_database, key, index)
Retrieves a string value from the database that matches the provided key.
Definition: optionsdatabase.F90:280
science_constants_mod
Scientific constant values used throughout simulations. Each has a default value and this can be over...
Definition: scienceconstants.F90:3
science_constants_mod::r
real(kind=default_precision), public r
Definition: scienceconstants.F90:13
grids_mod::z_index
integer, parameter, public z_index
Grid index parameters.
Definition: grids.F90:14
conversions_mod::conv_to_string
Converts data types to strings.
Definition: conversions.F90:38
gridmanager_mod::allocate_vertical_grid_data
subroutine allocate_vertical_grid_data(vertical_grid, n, nq)
Allocates the data required for the vertical grid configuration.
Definition: gridmanager.F90:649
q_indices_mod::standard_q_names
type(standard_q_names_type), public standard_q_names
Definition: q_indices.F90:59
state_mod::model_state_type
The ModelState which represents the current state of a run.
Definition: state.F90:39
gridmanager_mod::initialise_callback
subroutine initialise_callback(current_state)
Called during initialisation and will initialise the horizontal and vertical grid configurations Note...
Definition: gridmanager.F90:48
optionsdatabase_mod::options_get_logical
logical function, public options_get_logical(options_database, key, index)
Retrieves a logical value from the database that matches the provided key.
Definition: optionsdatabase.F90:154
gridmanager_mod::set_qv_init_from_rh
subroutine set_qv_init_from_rh(current_state)
Definition: gridmanager.F90:867
gridmanager_mod::set_vertical_reference_profile
subroutine set_vertical_reference_profile(current_state, vertical_grid, kkp)
Sets up the vertical grid reference profile at each point.
Definition: gridmanager.F90:111
q_indices_mod::get_q_index
integer function, public get_q_index(name, assigning_component)
Add in a new entry into the register if the name does not already exist or return the index of the pr...
Definition: q_indices.F90:112
logging_mod
Logging utility.
Definition: logging.F90:2
saturation_mod::dqwsatdt
real(kind=default_precision) function, public dqwsatdt(saturation_mixing_ratio, temperature)
Calculated the rate of change with temperature of saturation mixing ratio over liquid water....
Definition: saturation.F90:33
gridmanager_mod::setup_reference_state_liquid_water_temperature_and_saturation
subroutine setup_reference_state_liquid_water_temperature_and_saturation(current_state, vertical_grid, kkp)
Setting up reference state liquid water temperature and saturation mixing ratio on main levels.
Definition: gridmanager.F90:394
grids_mod::vertical_grid_configuration_type
The configuration of the grid vertically.
Definition: grids.F90:28
datadefn_mod
Contains common definitions for the data and datatypes used by MONC.
Definition: datadefn.F90:2
gridmanager_mod::set_buoyancy_coefficient
subroutine set_buoyancy_coefficient(current_state, vertical_grid, kkp)
Sets the buoyancy coefficient from the grid configuration and configuration.
Definition: gridmanager.F90:367
science_constants_mod::z0th
real(kind=default_precision), public z0th
Definition: scienceconstants.F90:13
gridmanager_mod::gridmanager_get_descriptor
type(component_descriptor_type) function, public gridmanager_get_descriptor()
Provides the descriptor back to the caller and is used in component registration.
Definition: gridmanager.F90:38
datadefn_mod::string_length
integer, parameter, public string_length
Default length of strings.
Definition: datadefn.F90:10
gridmanager_mod::compute_anelastic_pressure_profile_only
subroutine compute_anelastic_pressure_profile_only(current_state)
Computes the anelastic pressure only - if we are using Boussinesq approximation.
Definition: gridmanager.F90:719
logging_mod::log_master_log
subroutine, public log_master_log(level, message)
Will log just from the master process.
Definition: logging.F90:47
q_indices_mod
This manages the Q variables and specifically the mapping between names and the index that they are s...
Definition: q_indices.F90:2
gridmanager_mod::default_spacing
real, parameter default_spacing
The default spacing used if no grid is active in a specific dimension.
Definition: gridmanager.F90:28
gridmanager_mod::initialise_horizontalgrid_configuration_types
subroutine initialise_horizontalgrid_configuration_types(current_state)
Initialises the horizontal grid configurations.
Definition: gridmanager.F90:675
gridmanager_mod::set_anelastic_pressure
subroutine set_anelastic_pressure(current_state)
Set reference profile of potential temperature for the Boussinesq/Anelastic approximation Note that t...
Definition: gridmanager.F90:700
science_constants_mod::z0
real(kind=default_precision), public z0
Definition: scienceconstants.F90:13
optionsdatabase_mod::options_get_string_array
subroutine, public options_get_string_array(options_database, key, array_data, from, to)
Retrieves an entire (or subset) string array.
Definition: optionsdatabase.F90:302
gridmanager_mod::set_up_vertical_reference_properties
subroutine set_up_vertical_reference_properties(current_state, vertical_grid, kkp)
Sets up the reference properties for the vertical grid at each point.
Definition: gridmanager.F90:453
grids_mod
Functionality to support the different types of grid and abstraction between global grids and local o...
Definition: grids.F90:5
science_constants_mod::von_karman_constant
real(kind=default_precision), public von_karman_constant
Definition: scienceconstants.F90:13
gridmanager_mod::anelastic_profile_mode
integer, parameter anelastic_profile_mode
Definition: gridmanager.F90:27
optionsdatabase_mod
Manages the options database. Contains administration functions and deduce runtime options from the c...
Definition: optionsdatabase.F90:7
interpolation_mod
Definition: interpolation.F90:2
gridmanager_mod::set_up_and_smooth_grid
subroutine set_up_and_smooth_grid(vertical_grid, kgd, hgd, ninitp, kkp, zztop, nsmth, origional_setup, continuation_run)
Sets up and smooths the vertical grid. This is based upon the grid configuration already read in.
Definition: gridmanager.F90:482
optionsdatabase_mod::options_get_real_array
subroutine, public options_get_real_array(options_database, key, array_data, from, to)
Retrieves an entire (or subset) real array.
Definition: optionsdatabase.F90:113
gridmanager_mod
Manages the grid based upon the model state_mod.
Definition: gridmanager.F90:2
gridmanager_mod::calculate_initial_profiles
subroutine calculate_initial_profiles(current_state, vertical_grid)
Calculates the initial profiles for U, V, TH & Q if required.
Definition: gridmanager.F90:165
monc_component_mod::component_descriptor_type
Description of a component.
Definition: monc_component.F90:42
gridmanager_mod::original_vertical_grid_setup
subroutine original_vertical_grid_setup(vertical_grid, kgd, hgd, ninitp, kkp, zztop, nsmth)
The original vertical grid setup and smoothing as per the LEM.
Definition: gridmanager.F90:517
optionsdatabase_mod::options_get_real
real(kind=default_precision) function, public options_get_real(options_database, key, index)
Retrieves a real value from the database that matches the provided key.
Definition: optionsdatabase.F90:91
datadefn_mod::default_precision
integer, parameter, public default_precision
MPI communication type which we use for the prognostic and calculation data.
Definition: datadefn.F90:17
state_mod
The model state which represents the current state of a run.
Definition: state.F90:2
gridmanager_mod::check_input_levels
subroutine check_input_levels(z_levels, field_levels, field)
Definition: gridmanager.F90:855
science_constants_mod::rlvap_over_cp
real(kind=default_precision), public rlvap_over_cp
Definition: scienceconstants.F90:13
science_constants_mod::g
real(kind=default_precision), public g
Definition: scienceconstants.F90:13
Generated by doxygen 1.8.20