host.F90 53.7 KB
Newer Older
Jorn Bruggeman's avatar
Jorn Bruggeman committed
1 2 3
#include "fabm_driver.h"
#include "fabm_private.h"

4 5
module host_hooks
   use fabm_driver
Jorn Bruggeman's avatar
Jorn Bruggeman committed
6

7
   implicit none
Jorn Bruggeman's avatar
Jorn Bruggeman committed
8

9
   type,extends(type_base_driver) :: type_test_driver
Jorn Bruggeman's avatar
Jorn Bruggeman committed
10
   contains
11 12 13
      procedure :: fatal_error => test_driver_fatal_error
      procedure :: log_message => test_driver_log_message
   end type
Jorn Bruggeman's avatar
Jorn Bruggeman committed
14

15
contains
Jorn Bruggeman's avatar
Jorn Bruggeman committed
16

17 18 19
   subroutine test_driver_fatal_error(self,location,message)
      class (type_test_driver), intent(inout) :: self
      character(len=*),         intent(in)    :: location,message
Jorn Bruggeman's avatar
Jorn Bruggeman committed
20

21 22 23
      write (*,*) trim(location)//': '//trim(message)
      stop 1
   end subroutine
Jorn Bruggeman's avatar
Jorn Bruggeman committed
24

25 26 27
   subroutine test_driver_log_message(self,message)
      class (type_test_driver), intent(inout) :: self
      character(len=*),         intent(in)    :: message
Jorn Bruggeman's avatar
Jorn Bruggeman committed
28

29 30
      write (*,*) trim(message)
   end subroutine
Jorn Bruggeman's avatar
Jorn Bruggeman committed
31

32
end module host_hooks
Jorn Bruggeman's avatar
Jorn Bruggeman committed
33 34 35 36 37
   
program test_host

use fabm
use fabm_config
Jorn Bruggeman's avatar
Jorn Bruggeman committed
38
use fabm_driver
Jorn Bruggeman's avatar
Jorn Bruggeman committed
39 40

use test_models
41
use host_hooks
Jorn Bruggeman's avatar
Jorn Bruggeman committed
42 43 44 45 46

implicit none

#if _FABM_DIMENSION_COUNT_>0
integer :: _LOCATION_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
47 48
#endif

49 50 51 52
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
!  No mask but variable bottom index. Index of depth dimension must be 1.
!  All loops over inner dimension should skip points below bottom.
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
53 54
#    define _IMIN_ bottom_index _INDEX_HORIZONTAL_LOCATION_
#    define _IMAX_ domain_extent(1)
55
#  else
56 57
#    define _IMIN_ 1
#    define _IMAX_ bottom_index _INDEX_HORIZONTAL_LOCATION_
58 59 60
#  endif
#else
!  Loops over inner dimension should span full domain
61 62
#  define _IMIN_ 1
#  define _IMAX_ domain_extent(1)
63
#endif
64
#define _IRANGE_ _IMIN_,_IMAX_
65

66 67 68 69
#if _FABM_DIMENSION_COUNT_==0
#  define _BEGIN_GLOBAL_LOOP_
#  define _END_GLOBAL_LOOP_
#elif _FABM_DIMENSION_COUNT_==1
70 71 72
#  define _BEGIN_GLOBAL_LOOP_ do i__=_IRANGE_
#  define _END_GLOBAL_LOOP_ end do;i__=domain_extent(1)
#  ifdef _FABM_DEPTH_DIMENSION_INDEX_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
73
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_
74
#    define _END_GLOBAL_HORIZONTAL_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
75
#  endif
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
#elif _FABM_DIMENSION_COUNT_==2
#  define _BEGIN_GLOBAL_LOOP_ do j__=1,domain_extent(2);do i__=_IRANGE_
#  define _END_GLOBAL_LOOP_ end do;end do;i__=domain_extent(1);j__=domain_extent(2)
#  if _FABM_DEPTH_DIMENSION_INDEX_==1
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ do j__=1,domain_extent(2)
#    define _END_GLOBAL_HORIZONTAL_LOOP_ end do;j__=domain_extent(2)
#  elif _FABM_DEPTH_DIMENSION_INDEX_==2
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ do i__=_IRANGE_
#    define _END_GLOBAL_HORIZONTAL_LOOP_ end do;i__=domain_extent(1)
#  endif
#elif _FABM_DIMENSION_COUNT_==3
#  define _BEGIN_GLOBAL_LOOP_ do k__=1,domain_extent(3);do j__=1,domain_extent(2);do i__=_IRANGE_
#  define _END_GLOBAL_LOOP_ end do;end do;end do;i__=domain_extent(1);j__=domain_extent(2);k__=domain_extent(3)
#  if _FABM_DEPTH_DIMENSION_INDEX_==1
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ do k__=1,domain_extent(3);do j__=1,domain_extent(2)
#    define _END_GLOBAL_HORIZONTAL_LOOP_ end do;end do;j__=domain_extent(2);k__=domain_extent(3)
#  elif _FABM_DEPTH_DIMENSION_INDEX_==2
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ do k__=1,domain_extent(3);do i__=_IRANGE_
#    define _END_GLOBAL_HORIZONTAL_LOOP_ end do;end do;i__=domain_extent(1);k__=domain_extent(3)
#  elif _FABM_DEPTH_DIMENSION_INDEX_==3
#    define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ do j__=1,domain_extent(2);do i__=_IRANGE_
#    define _END_GLOBAL_HORIZONTAL_LOOP_ end do;end do;i__=domain_extent(1);j__=domain_extent(2)
#  endif
#endif

! If there is no depth dimension, horizontal = global
#ifndef _FABM_DEPTH_DIMENSION_INDEX_
#  define _BEGIN_GLOBAL_HORIZONTAL_LOOP_ _BEGIN_GLOBAL_LOOP_
#  define _END_GLOBAL_HORIZONTAL_LOOP_ _END_GLOBAL_LOOP_
#endif

#ifndef _FABM_VECTORIZED_DIMENSION_INDEX_
   ! No vectorization: outer loops are global loops
#  define _BEGIN_OUTER_INTERIOR_LOOP_ _BEGIN_GLOBAL_LOOP_
#  define _END_OUTER_INTERIOR_LOOP_ _END_GLOBAL_LOOP_
#  define _BEGIN_OUTER_HORIZONTAL_LOOP_ _BEGIN_GLOBAL_HORIZONTAL_LOOP_
#  define _END_OUTER_HORIZONTAL_LOOP_ _END_GLOBAL_HORIZONTAL_LOOP_
#elif _FABM_DIMENSION_COUNT_==1
   ! Entire domain is vectorized; no outer loops needed
#  define _BEGIN_OUTER_INTERIOR_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
116
#  define _END_OUTER_INTERIOR_LOOP_
117 118
#  define _BEGIN_OUTER_HORIZONTAL_LOOP_
#  define _END_OUTER_HORIZONTAL_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
119
#elif _FABM_DIMENSION_COUNT_==2
120 121
#  define _BEGIN_OUTER_INTERIOR_LOOP_ do j__=1,domain_extent(2)
#  define _END_OUTER_INTERIOR_LOOP_ end do;j__=domain_extent(2)
122
#  if _FABM_DEPTH_DIMENSION_INDEX_==2
123
     ! The entire horizontal is already vectorized; no outer loop necessary
124
#    define _BEGIN_OUTER_HORIZONTAL_LOOP_
125 126 127 128 129
#    define _END_OUTER_HORIZONTAL_LOOP_
#  else
     ! No horizontal dimension vectorized; do full outer loop.
#    define _BEGIN_OUTER_HORIZONTAL_LOOP_ _BEGIN_GLOBAL_HORIZONTAL_LOOP_
#    define _END_OUTER_HORIZONTAL_LOOP_ _END_GLOBAL_HORIZONTAL_LOOP_
130 131
#  endif
#elif _FABM_DIMENSION_COUNT_==3
132 133 134 135 136 137
#  define _BEGIN_OUTER_INTERIOR_LOOP_ do k__=1,domain_extent(3);do j__=1,domain_extent(2)
#  define _END_OUTER_INTERIOR_LOOP_ end do;end do;j__=domain_extent(2);k__=domain_extent(3)
#  if _FABM_DEPTH_DIMENSION_INDEX_==2
#    define _BEGIN_OUTER_HORIZONTAL_LOOP_ do k__=1,domain_extent(3)
#    define _END_OUTER_HORIZONTAL_LOOP_ end do;k__=domain_extent(3)
#  elif _FABM_DEPTH_DIMENSION_INDEX_==3
138
#    define _BEGIN_OUTER_HORIZONTAL_LOOP_ do j__=1,domain_extent(2)
139 140 141 142 143
#    define _END_OUTER_HORIZONTAL_LOOP_ end do;j__=domain_extent(2)
#  else
     ! No horizontal dimension vectorized; do full outer loop.
#    define _BEGIN_OUTER_HORIZONTAL_LOOP_ _BEGIN_GLOBAL_HORIZONTAL_LOOP_
#    define _END_OUTER_HORIZONTAL_LOOP_ _END_GLOBAL_HORIZONTAL_LOOP_
144 145 146
#  endif
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
147
#ifdef _INTERIOR_IS_VECTORIZED_
148
integer :: _START_, _STOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
149 150
#endif

151 152 153
real(rk),allocatable _DIMENSION_GLOBAL_ :: tmp
real(rk),allocatable _DIMENSION_GLOBAL_HORIZONTAL_ :: tmp_hz

Jorn Bruggeman's avatar
Jorn Bruggeman committed
154
#ifdef _HAS_MASK_
155
_FABM_MASK_TYPE_,allocatable,target _DIMENSION_GLOBAL_HORIZONTAL_ :: mask_hz
Jorn Bruggeman's avatar
Jorn Bruggeman committed
156 157 158
#  ifndef _FABM_HORIZONTAL_MASK_
_FABM_MASK_TYPE_,allocatable,target _DIMENSION_GLOBAL_ :: mask
#  endif
159 160 161
#  ifndef _FABM_MASKED_VALUE_
#    define _FABM_MASKED_VALUE_ _FABM_UNMASKED_VALUE_+1
#  endif
162 163 164
#  ifndef _FABM_UNMASKED_VALUE_
#    define _FABM_UNMASKED_VALUE_ _FABM_MASKED_VALUE_+1
#  endif
165 166
#endif

167 168 169
integer :: interior_count
integer :: horizontal_count

170 171
#if _FABM_BOTTOM_INDEX_==-1
integer,allocatable,target _DIMENSION_GLOBAL_HORIZONTAL_ :: bottom_index
Jorn Bruggeman's avatar
Jorn Bruggeman committed
172 173 174 175 176 177
#endif

real(rk),allocatable,target _DIMENSION_GLOBAL_PLUS_1_            :: interior_state
real(rk),allocatable,target _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: surface_state
real(rk),allocatable,target _DIMENSION_GLOBAL_HORIZONTAL_PLUS_1_ :: bottom_state

Jorn Bruggeman's avatar
Jorn Bruggeman committed
178 179
real(rk),allocatable _DIMENSION_SLICE_PLUS_1_        :: dy
real(rk),allocatable _DIMENSION_SLICE_PLUS_1_        :: w
Jorn Bruggeman's avatar
Jorn Bruggeman committed
180
real(rk),allocatable _DIMENSION_HORIZONTAL_SLICE_PLUS_1_ :: flux,sms_sf,sms_bt
Jorn Bruggeman's avatar
Jorn Bruggeman committed
181

Jorn Bruggeman's avatar
Jorn Bruggeman committed
182 183
real(rk),allocatable,target _DIMENSION_GLOBAL_            :: temperature
real(rk),allocatable,target _DIMENSION_GLOBAL_            :: depth
Jorn Bruggeman's avatar
Jorn Bruggeman committed
184 185 186
real(rk),allocatable,target _DIMENSION_GLOBAL_HORIZONTAL_ :: wind_speed

type (type_model) :: model
Jorn Bruggeman's avatar
Jorn Bruggeman committed
187

Jorn Bruggeman's avatar
Jorn Bruggeman committed
188 189
class (type_test_model), pointer :: test_model

Jorn Bruggeman's avatar
Jorn Bruggeman committed
190 191
integer :: domain_extent(_FABM_DIMENSION_COUNT_)

192
character(len=20) :: arg
Jorn Bruggeman's avatar
Jorn Bruggeman committed
193 194
integer :: ivar
integer :: i
195
integer :: mode = 1
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
integer :: ntest = -1

! Parse command line arguments
i = 1
do
   call get_command_argument(i, arg)
   if (arg == '') exit
   select case (arg)
   case ('-s', '--simulate')
      mode = 2
      i = i + 1
   case ('-n')
      call get_command_argument(i + 1, arg)
      read (arg,*) ntest
      i = i + 2
   case default
      write (*,'(a)') 'Unknown command line argument: ' // trim(arg)
      stop 2
   end select
end do
216

Jorn Bruggeman's avatar
Jorn Bruggeman committed
217
#if _FABM_DIMENSION_COUNT_>0
Jorn Bruggeman's avatar
Jorn Bruggeman committed
218
i__=50
Jorn Bruggeman's avatar
Jorn Bruggeman committed
219 220
#endif
#if _FABM_DIMENSION_COUNT_>1
Jorn Bruggeman's avatar
Jorn Bruggeman committed
221
j__=40
Jorn Bruggeman's avatar
Jorn Bruggeman committed
222 223
#endif
#if _FABM_DIMENSION_COUNT_>2
Jorn Bruggeman's avatar
Jorn Bruggeman committed
224
k__=45
Jorn Bruggeman's avatar
Jorn Bruggeman committed
225 226 227
#endif

#if _FABM_DIMENSION_COUNT_>0
Jorn Bruggeman's avatar
Jorn Bruggeman committed
228 229 230
domain_extent = (/ _LOCATION_ /)
#endif

231 232 233 234 235 236 237 238 239
! Set defaults
if (ntest == -1) then
   if (mode == 1) then
      ntest = 1
   else
      ntest = 50000000/product(domain_extent)
   end if
end if

240 241 242 243 244 245 246
interior_count = product(domain_extent)
#ifdef _FABM_DEPTH_DIMENSION_INDEX_
horizontal_count = interior_count / domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)
#else
horizontal_count = interior_count
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
247
#ifdef _INTERIOR_IS_VECTORIZED_
248 249
_START_ = 1
_STOP_ = domain_extent(_FABM_VECTORIZED_DIMENSION_INDEX_)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
250 251
#endif

252 253 254 255
allocate(tmp _INDEX_LOCATION_)
allocate(tmp_hz _INDEX_HORIZONTAL_LOCATION_)

allocate(type_test_driver::driver)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
256 257
call fabm_initialize_library()

258 259 260 261 262 263 264 265 266
select case (mode)
case (1)
    ! Unit testing with built-in model
    allocate(test_model)
    call model%root%add_child(test_model,'test_model','test model',configunit=-1)
case (2)
    ! Test with user-provided fabm.yaml
    call fabm_create_model_from_yaml_file(model, do_not_initialize=.true.)
end select
267

Jorn Bruggeman's avatar
Jorn Bruggeman committed
268
call start_test('fabm_initialize')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
269
call fabm_initialize(model)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
270
call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
271

Jorn Bruggeman's avatar
Jorn Bruggeman committed
272
! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
273
! Provide extents of the spatial domain.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
274 275 276
! ======================================================================

call start_test('fabm_set_domain')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
277
call fabm_set_domain(model _ARGUMENTS_LOCATION_)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
278 279 280 281 282
call report_test_result()

! ======================================================================
! Set up spatial mask.
! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
283 284

#ifdef _HAS_MASK_
285 286
allocate(mask_hz _INDEX_HORIZONTAL_LOCATION_)
call start_test('fabm_set_mask')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
287
#  ifdef _FABM_HORIZONTAL_MASK_
288
call fabm_set_mask(model,mask_hz)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
289 290
#  else
allocate(mask _INDEX_LOCATION_)
291
call fabm_set_mask(model,mask,mask_hz)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
292
#  endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
293
call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
294 295
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
296 297 298 299
! ======================================================================
! Specify vertical indices of surface and bottom.
! ======================================================================

Jorn Bruggeman's avatar
Jorn Bruggeman committed
300
#ifdef _FABM_DEPTH_DIMENSION_INDEX_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
301
call start_test('set_surface_index')
302
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
303
call model%set_surface_index(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
304 305
#  else
call model%set_surface_index(1)
306
#  endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
307 308 309
call report_test_result()

call start_test('set_bottom_index')
310
#  if _FABM_BOTTOM_INDEX_==-1
311 312
allocate(bottom_index _INDEX_HORIZONTAL_LOCATION_)
call model%set_bottom_index(bottom_index)
313 314 315
#  elif defined(_FABM_VERTICAL_BOTTOM_TO_SURFACE_)
call model%set_bottom_index(1)
#  else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
316
call model%set_bottom_index(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
317
#  endif
318
call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
319 320 321 322 323 324
#endif

allocate(interior_state(_PREARG_LOCATION_ size(model%state_variables)))
allocate(surface_state(_PREARG_HORIZONTAL_LOCATION_ size(model%surface_state_variables)))
allocate(bottom_state(_PREARG_HORIZONTAL_LOCATION_ size(model%bottom_state_variables)))

Jorn Bruggeman's avatar
Jorn Bruggeman committed
325 326 327
! ======================================================================
! Send pointers to state variable data to FABM.
! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
328

Jorn Bruggeman's avatar
Jorn Bruggeman committed
329
call start_test('link_interior_state_data')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
330 331 332
do ivar=1,size(model%state_variables)
   call model%link_interior_state_data(ivar,interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar))
end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
333 334 335
call report_test_result()

call start_test('link_surface_state_data')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
336 337 338
do ivar=1,size(model%surface_state_variables)
   call model%link_surface_state_data(ivar,surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar))
end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
339 340 341
call report_test_result()

call start_test('link_bottom_state_data')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
342 343 344
do ivar=1,size(model%bottom_state_variables)
   call model%link_bottom_state_data(ivar,bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar))
end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
345
call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
346

Jorn Bruggeman's avatar
Jorn Bruggeman committed
347 348 349 350
! ======================================================================
! Transfer pointers to environmental data
! ======================================================================

351 352 353 354 355
select case (mode)
case (1)
    allocate(depth _INDEX_LOCATION_)
    allocate(temperature _INDEX_LOCATION_)
    allocate(wind_speed _INDEX_HORIZONTAL_LOCATION_)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
356

357 358 359 360 361 362 363 364 365 366 367
    call start_test('link_interior_data')
    call model%link_interior_data(standard_variables%temperature,temperature)
    call model%link_interior_data(standard_variables%depth,depth)
    call report_test_result()

    call start_test('link_horizontal_data')
    call model%link_horizontal_data(standard_variables%wind_speed,wind_speed)
    call report_test_result()
case (2)
    call read_environment
end select
Jorn Bruggeman's avatar
Jorn Bruggeman committed
368

Jorn Bruggeman's avatar
Jorn Bruggeman committed
369
! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
370 371
! Check whether FABM has all dependencies fulfilled
! (i.e., whether all required calls for fabm_link_*_data have been made)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
372
! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
373

Jorn Bruggeman's avatar
Jorn Bruggeman committed
374 375 376
call start_test('fabm_check_ready')
call fabm_check_ready(model)
call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
377 378

#ifdef _FABM_VECTORIZED_DIMENSION_INDEX_
379 380
allocate(dy(_START_:_STOP_,size(model%state_variables)))
allocate(w(_START_:_STOP_,size(model%state_variables)))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
381 382 383 384 385 386
#else
allocate(dy(size(model%state_variables)))
allocate(w(size(model%state_variables)))
#endif

#if defined(_FABM_VECTORIZED_DIMENSION_INDEX_)&&(_FABM_DEPTH_DIMENSION_INDEX_!=_FABM_VECTORIZED_DIMENSION_INDEX_)
387 388 389
allocate(flux(_START_:_STOP_,size(model%state_variables)))
allocate(sms_sf(_START_:_STOP_,size(model%surface_state_variables)))
allocate(sms_bt(_START_:_STOP_,size(model%bottom_state_variables)))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
390
#else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
391 392 393
allocate(flux(size(model%state_variables)))
allocate(sms_sf(size(model%surface_state_variables)))
allocate(sms_bt(size(model%bottom_state_variables)))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
394 395
#endif

396 397
select case (mode)
case (1)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
398 399 400
   do i=1,ntest
      call test_update
   end do
401
case(2)
402
   call simulate(ntest)
403
end select
Jorn Bruggeman's avatar
Jorn Bruggeman committed
404

405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473
contains

   subroutine read_environment
      use yaml, only: yaml_parse => parse, yaml_error_length => error_length
      use yaml_types, only: type_yaml_node => type_node, type_yaml_dictionary => type_dictionary, type_yaml_scalar => type_scalar, type_yaml_key_value_pair => type_key_value_pair

      integer, parameter :: yaml_unit = 100
      character(yaml_error_length) :: yaml_error
      class (type_yaml_node),pointer :: yaml_root
      type (type_yaml_key_value_pair), pointer :: yaml_pair
      real(rk) :: value
      logical :: success
      type type_input
         type (type_bulk_variable_id)                        :: interior_id
         type (type_horizontal_variable_id)                  :: horizontal_id
         type (type_scalar_variable_id)                      :: scalar_id
         real(rk), allocatable _DIMENSION_GLOBAL_            :: interior_data
         real(rk), allocatable _DIMENSION_GLOBAL_HORIZONTAL_ :: horizontal_data
         real(rk)                                            :: scalar_data
      end type
      type (type_input), pointer :: input

      yaml_root => yaml_parse('environment.yaml', yaml_unit, yaml_error)
      if (yaml_error /= '') then
         call driver%log_message(yaml_error)
         stop 2
      end if
      select type (yaml_root)
      class is (type_yaml_dictionary)
          yaml_pair => yaml_root%first
          do while (associated(yaml_pair))
              select type (node => yaml_pair%value)
              class is (type_yaml_scalar)
                  call driver%log_message('Setting '//trim(yaml_pair%key)//' to '//trim(node%string))
                  value = node%to_real(0._rk, success)
                  if (.not. success) then
                     call driver%log_message('Cannot parse '//trim(node%string)//' as real.')
                     stop 2
                  end if
                  allocate(input)
                  input%interior_id = model%get_bulk_variable_id(trim(yaml_pair%key))
                  if (fabm_is_variable_used(input%interior_id)) then
                      allocate(input%interior_data _INDEX_LOCATION_)
                      input%interior_data = value
                      call model%link_interior_data(input%interior_id, input%interior_data)
                  else
                      input%horizontal_id = model%get_horizontal_variable_id(trim(yaml_pair%key))
                      if (fabm_is_variable_used(input%horizontal_id)) then
                         allocate(input%horizontal_data _INDEX_HORIZONTAL_LOCATION_)
                         input%horizontal_data = value
                         call model%link_horizontal_data(input%horizontal_id, input%horizontal_data)
                      else
                         input%scalar_id = model%get_scalar_variable_id(trim(yaml_pair%key))
                         if (fabm_is_variable_used(input%scalar_id)) then
                            input%scalar_data = value
                            call model%link_scalar(input%scalar_id, input%scalar_data)
                         else
                            call driver%log_message('WARNING: environment variable '//trim(yaml_pair%key)//' is not used by FABM model and will be ignored.')
                         end if
                      end if
                  end if
              end select
              yaml_pair => yaml_pair%next
          end do
      class default
         call driver%log_message('environment.yaml should contain a dictionary at root level')
         stop 2
      end select
   end subroutine read_environment
Jorn Bruggeman's avatar
Jorn Bruggeman committed
474

Jorn Bruggeman's avatar
Jorn Bruggeman committed
475
   subroutine randomize_mask
476 477 478 479 480 481 482 483 484 485
#if _FABM_BOTTOM_INDEX_==-1
      ! Depth index of bottom varies in the horizontal - pick random numbers between 0 (land) and maximum index
      call random_number(tmp_hz)
      bottom_index = floor(tmp_hz*(1+domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)))
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
     ! Ensure invalid bottom indices [land points] are set such that vertical loops have 0 iterations.
     where (bottom_index == 0) bottom_index = domain_extent(_FABM_DEPTH_DIMENSION_INDEX_) + 1
#  endif
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
486
#ifdef _HAS_MASK_
487 488 489 490 491
#  ifdef _FABM_HORIZONTAL_MASK_
      ! Apply random mask across horizontal domain (half of grid cells masked)
      call random_number(tmp_hz)
      mask_hz = _FABM_UNMASKED_VALUE_
      where (tmp_hz>0.5_rk) mask_hz = _FABM_MASKED_VALUE_
492 493 494 495 496 497
      horizontal_count = count(_IS_UNMASKED_(mask_hz))
#    ifdef _FABM_DEPTH_DIMENSION_INDEX_
      interior_count = horizontal_count * domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)
#    else
      interior_count = horizontal_count
#    endif
498 499 500
#  else
      ! Apply random mask across interior domain (half of grid cells masked)
      call random_number(tmp)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
501
      mask = _FABM_UNMASKED_VALUE_
502 503 504
      where (tmp>0.5_rk) mask = _FABM_MASKED_VALUE_

#    if _FABM_BOTTOM_INDEX_==-1
505
      ! Bottom index varies in the horizontal. Ensure the bottom cell itself is unmasked, and anything deeper is masked.
506
      _BEGIN_GLOBAL_HORIZONTAL_LOOP_
507
         ! Valid bottom index - unmask associated cell, then mask all deeper ones
508
#      ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
509
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) mask _INDEX_GLOBAL_VERTICAL_(bottom_index _INDEX_HORIZONTAL_LOCATION_) = _FABM_UNMASKED_VALUE_
510
         mask _INDEX_GLOBAL_VERTICAL_(:bottom_index _INDEX_HORIZONTAL_LOCATION_ - 1) = _FABM_MASKED_VALUE_
511
#      else
512
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1) mask _INDEX_GLOBAL_VERTICAL_(bottom_index _INDEX_HORIZONTAL_LOCATION_) = _FABM_UNMASKED_VALUE_
513
         mask _INDEX_GLOBAL_VERTICAL_(bottom_index _INDEX_HORIZONTAL_LOCATION_ + 1:) = _FABM_MASKED_VALUE_
514
#      endif
515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
      _END_GLOBAL_HORIZONTAL_LOOP_
#    endif

      ! Infer horizontal mask (mask points that have all column layers masked)
      mask_hz = _FABM_UNMASKED_VALUE_
      where (.not.any(_IS_UNMASKED_(mask),_FABM_DEPTH_DIMENSION_INDEX_)) mask_hz = _FABM_MASKED_VALUE_

      ! For valid points in the horizontal, make sure that index 1 (bottom or surface) is unmasked
      _BEGIN_GLOBAL_HORIZONTAL_LOOP_
         if (_IS_UNMASKED_(mask_hz _INDEX_HORIZONTAL_LOCATION_)) then
            mask _INDEX_GLOBAL_VERTICAL_(1) = _FABM_UNMASKED_VALUE_
#    if _FABM_BOTTOM_INDEX_!=-1
            mask _INDEX_GLOBAL_VERTICAL_(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) = _FABM_UNMASKED_VALUE_
#    endif
         end if
      _END_GLOBAL_HORIZONTAL_LOOP_
531 532
      horizontal_count = count(_IS_UNMASKED_(mask_hz))
      interior_count = count(_IS_UNMASKED_(mask))
533
#  endif
534
#elif _FABM_BOTTOM_INDEX_==-1
535
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
536
      horizontal_count = count(bottom_index <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_))
537 538
      interior_count = sum(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_) - bottom_index + 1)
#  else
539
      horizontal_count = count(bottom_index >= 1)
540 541
      interior_count = sum(bottom_index)
#  endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
542 543 544
#endif
   end subroutine randomize_mask

545 546
   subroutine simulate(n)
      integer, intent(in) :: n
547
      real(rk) :: time_begin, time_end
548 549 550 551 552 553 554 555 556
      integer :: nseed
      integer, allocatable :: seed(:)

      call random_seed(size=nseed)
      allocate(seed(nseed))
      seed(:) = 1
      call random_seed(put=seed)

      call randomize_mask
557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575

      call start_test('fabm_initialize_state')
      _BEGIN_OUTER_INTERIOR_LOOP_
         call fabm_initialize_state(model _ARGUMENTS_INTERIOR_IN_)
      _END_OUTER_INTERIOR_LOOP_
      call report_test_result()

      call start_test('fabm_initialize_bottom_state')
      _BEGIN_OUTER_HORIZONTAL_LOOP_
         call fabm_initialize_bottom_state(model _ARGUMENTS_HORIZONTAL_IN_)
      _END_OUTER_HORIZONTAL_LOOP_
      call report_test_result()

      call start_test('fabm_initialize_surface_state')
      _BEGIN_OUTER_HORIZONTAL_LOOP_
         call fabm_initialize_surface_state(model _ARGUMENTS_HORIZONTAL_IN_)
      _END_OUTER_HORIZONTAL_LOOP_
      call report_test_result()

576
      write (*,'(a,i0,a)') 'Simulating with ', interior_count, ' wet cells...'
577 578 579

      call cpu_time(time_begin)

580
      do i=1,n
581 582
         _BEGIN_GLOBAL_HORIZONTAL_LOOP_
#ifdef _FABM_DEPTH_DIMENSION_INDEX_
583 584 585 586 587
#  if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
            ! No mask but non-constant bottom index. We need to skip everything below bottom
            if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) &
               call fabm_get_light(model,_IRANGE_ _ARG_VERTICAL_FIXED_LOCATION_)
#  else
588
            call fabm_get_light(model,1,domain_extent(_FABM_DEPTH_DIMENSION_INDEX_) _ARG_VERTICAL_FIXED_LOCATION_)
589
#  endif
590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610
#else
            call fabm_get_light(model _ARGUMENTS_HORIZONTAL_IN_)
#endif
         _END_GLOBAL_HORIZONTAL_LOOP_

         _BEGIN_OUTER_HORIZONTAL_LOOP_
            flux = 0
            sms_bt = 0
            call fabm_do_bottom(model _ARGUMENTS_HORIZONTAL_IN_,flux,sms_bt)
         _END_OUTER_HORIZONTAL_LOOP_

         _BEGIN_OUTER_HORIZONTAL_LOOP_
            flux = 0
            sms_sf = 0
            call fabm_do_surface(model _ARGUMENTS_HORIZONTAL_IN_,flux,sms_sf)
         _END_OUTER_HORIZONTAL_LOOP_

         _BEGIN_OUTER_INTERIOR_LOOP_
            dy = 0
            call fabm_do(model _ARGUMENTS_INTERIOR_IN_,dy)
         _END_OUTER_INTERIOR_LOOP_
611 612

         if (mod(i, 100) == 0) write (*,'(i0,a)') int(100*i/real(n, rk)), ' % complete'
613 614 615 616
      end do

      call cpu_time(time_end)

617 618
      write (*,'(a)') 'Simulation complete.'
      write (*,'(a,f8.3,a)') 'Total time: ', time_end - time_begin, ' s'
619 620
   end subroutine

Jorn Bruggeman's avatar
Jorn Bruggeman committed
621
   subroutine test_update
Jorn Bruggeman's avatar
Jorn Bruggeman committed
622 623
      real(rk),pointer _DIMENSION_GLOBAL_ :: pdata
      real(rk),pointer _DIMENSION_GLOBAL_HORIZONTAL_ :: pdata_hz
624
      logical :: valid
Jorn Bruggeman's avatar
Jorn Bruggeman committed
625

Jorn Bruggeman's avatar
Jorn Bruggeman committed
626 627 628 629 630
      call randomize_mask

      ! ======================================================================
      ! Initialize all state variables
      ! ======================================================================
Jorn Bruggeman's avatar
Jorn Bruggeman committed
631

Jorn Bruggeman's avatar
Jorn Bruggeman committed
632
      call start_test('fabm_initialize_state')
633
      _BEGIN_OUTER_INTERIOR_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
634
         call fabm_initialize_state(model _ARGUMENTS_INTERIOR_IN_)
635
      _END_OUTER_INTERIOR_LOOP_
636 637 638
      do ivar=1,size(model%state_variables)
         call check_interior(interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar), model%state_variables(ivar)%missing_value, ivar+interior_state_offset+1._rk)
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
639
      call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
640

Jorn Bruggeman's avatar
Jorn Bruggeman committed
641
      call start_test('fabm_initialize_bottom_state')
642 643 644
      _BEGIN_OUTER_HORIZONTAL_LOOP_
         call fabm_initialize_bottom_state(model _ARGUMENTS_HORIZONTAL_IN_)
      _END_OUTER_HORIZONTAL_LOOP_
645 646 647
      do ivar=1,size(model%bottom_state_variables)
         call check_horizontal(bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar), model%bottom_state_variables(ivar)%missing_value, ivar+bottom_state_offset+1._rk)
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
648
      call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
649

Jorn Bruggeman's avatar
Jorn Bruggeman committed
650
      call start_test('fabm_initialize_surface_state')
651 652 653
      _BEGIN_OUTER_HORIZONTAL_LOOP_
         call fabm_initialize_surface_state(model _ARGUMENTS_HORIZONTAL_IN_)
      _END_OUTER_HORIZONTAL_LOOP_
654 655 656
      do ivar=1,size(model%surface_state_variables)
         call check_horizontal(surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar), model%surface_state_variables(ivar)%missing_value, ivar+surface_state_offset+1._rk)
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
657
      call report_test_result()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
658

Jorn Bruggeman's avatar
Jorn Bruggeman committed
659 660 661 662 663 664 665 666 667 668 669 670
      ! ======================================================================
      ! Initialize environmental dependencies
      ! ======================================================================

      temperature = 1+interior_dependency_offset
      call apply_mask_3d(temperature,-999._rk-interior_dependency_offset)

      wind_speed = 1+horizontal_dependency_offset
      call apply_mask_2d(wind_speed,-999._rk-horizontal_dependency_offset)

#ifdef _FABM_DEPTH_DIMENSION_INDEX_
      ! Model has depth dimension: make sure depth varies from 0 at the surface till 1 at the bottom
671
      _BEGIN_GLOBAL_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
672
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
673 674 675 676 677 678 679
#    if _FABM_BOTTOM_INDEX_==-1
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ == domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) then
            depth _INDEX_LOCATION_ = 2
         else
            depth _INDEX_LOCATION_ = real(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)-_VERTICAL_ITERATOR_,rk)/(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)-bottom_index _INDEX_HORIZONTAL_LOCATION_)
         end if
#    else
680
         depth _INDEX_LOCATION_ = real(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)-_VERTICAL_ITERATOR_,rk)/(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)-1)
681
#    endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
682
#  else
683
#    if _FABM_BOTTOM_INDEX_==-1
684 685 686 687 688
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_==1) then
            depth _INDEX_LOCATION_ = 2
         else
            depth _INDEX_LOCATION_ = real(_VERTICAL_ITERATOR_-1,rk)/(bottom_index _INDEX_HORIZONTAL_LOCATION_-1)
         end if
689
#    else
690
         depth _INDEX_LOCATION_ = real(_VERTICAL_ITERATOR_-1,rk)/(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)-1)
691
#    endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
692
#  endif
693
      _END_GLOBAL_LOOP_
Jorn Bruggeman's avatar
Jorn Bruggeman committed
694 695 696 697 698 699 700 701
#else
      ! No depth dimension
      depth = 2
#endif
      call apply_mask_3d(depth,-999._rk-interior_dependency_offset)

      do ivar=1,size(model%state_variables)
         interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar) = ivar+interior_state_offset
702
         call apply_mask_3d(interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar),model%state_variables(ivar)%missing_value)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
703 704 705
      end do
      do ivar=1,size(model%surface_state_variables)
         surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar) = ivar+surface_state_offset
706
         call apply_mask_2d(surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar),model%surface_state_variables(ivar)%missing_value)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
707 708 709
      end do
      do ivar=1,size(model%bottom_state_variables)
         bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar) = ivar+bottom_state_offset
710
         call apply_mask_2d(bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar),model%bottom_state_variables(ivar)%missing_value)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
711 712 713 714 715 716 717
      end do

      ! ======================================================================
      ! Retrieve source terms of interior state variables.
      ! ======================================================================

      call start_test('fabm_do')
718
      interior_loop_count = 0
719 720
      _BEGIN_OUTER_INTERIOR_LOOP_
         dy = 0
721 722 723 724 725
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_) && _FABM_VECTORIZED_DIMENSION_INDEX_==_FABM_DEPTH_DIMENSION_INDEX_ && defined(_FABM_DEPTH_DIMENSION_INDEX_)
         ! We are looping over depth, but as we have a non-constant bottom index (yet no mask), we need to skip everything below bottom
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) &
            call fabm_do(model,_IMIN_,_IMAX_ _ARG_INTERIOR_FIXED_LOCATION_,dy(_IMIN_:_IMAX_,:))
#else
726
         call fabm_do(model _ARGUMENTS_INTERIOR_IN_,dy)
727
#endif
728 729 730 731
         do ivar=1,size(model%state_variables)
            call check_interior_slice_plus_1(dy,ivar,0.0_rk,-real(ivar+interior_state_offset,rk) _ARGUMENTS_INTERIOR_IN_)
         end do
      _END_OUTER_INTERIOR_LOOP_
732
      call assert(interior_loop_count == interior_count, 'fabm_do', 'call count does not match number of (unmasked) interior points')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
733 734 735 736 737 738 739 740

      do ivar=1,size(model%diagnostic_variables)
         if (model%diagnostic_variables(ivar)%save .and. model%diagnostic_variables(ivar)%target%source==source_do) then
            pdata => fabm_get_interior_diagnostic_data(model, ivar)
            call check_interior(pdata, model%diagnostic_variables(ivar)%missing_value, -model%diagnostic_variables(ivar)%missing_value)
         end if
      end do

Jorn Bruggeman's avatar
Jorn Bruggeman committed
741 742 743 744 745 746
      call report_test_result()

      ! ======================================================================
      ! Retrieve surface fluxes of interior state variables, source terms of surface-associated state variables.
      ! ======================================================================

747 748
#if _FABM_BOTTOM_INDEX_==-1
      _BEGIN_GLOBAL_HORIZONTAL_LOOP_
749 750 751 752 753
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ == domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) depth _INDEX_GLOBAL_VERTICAL_(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) = 0
#  else
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ == 1) depth _INDEX_GLOBAL_VERTICAL_(1) = 0
#  endif
754 755 756
      _END_GLOBAL_HORIZONTAL_LOOP_
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
757
      call start_test('fabm_do_surface')
758
      surface_loop_count = 0
759
      _BEGIN_OUTER_HORIZONTAL_LOOP_
760 761 762
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) then
#endif
763 764 765 766 767 768 769 770 771
         flux = 0
         sms_sf = 0
         call fabm_do_surface(model _ARGUMENTS_HORIZONTAL_IN_,flux,sms_sf)
         do ivar=1,size(model%state_variables)
            call check_horizontal_slice_plus_1(flux,ivar,0.0_rk,-real(ivar+interior_state_offset,rk) _ARGUMENTS_HORIZONTAL_IN_)
         end do
         do ivar=1,size(model%surface_state_variables)
            call check_horizontal_slice_plus_1(sms_sf,ivar,0.0_rk,-real(ivar+surface_state_offset,rk) _ARGUMENTS_HORIZONTAL_IN_)
         end do
772 773 774
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
         endif
#endif
775
      _END_OUTER_HORIZONTAL_LOOP_
776
      call assert(surface_loop_count == horizontal_count, 'fabm_do_surface', 'call count does not match number of (unmasked) horizontal points')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
777

Jorn Bruggeman's avatar
Jorn Bruggeman committed
778 779 780 781 782 783 784 785
      do ivar=1,size(model%horizontal_diagnostic_variables)
         if (model%horizontal_diagnostic_variables(ivar)%save .and. model%horizontal_diagnostic_variables(ivar)%target%source == source_do_surface) then
            pdata_hz => fabm_get_horizontal_diagnostic_data(model, ivar)
            call check_horizontal(pdata_hz, model%horizontal_diagnostic_variables(ivar)%missing_value, -model%horizontal_diagnostic_variables(ivar)%missing_value)
         end if
      end do

      call report_test_result()
786

Jorn Bruggeman's avatar
Jorn Bruggeman committed
787 788 789 790
      ! ======================================================================
      ! Retrieve bottom fluxes of interior state variables, source terms of bottom-associated state variables.
      ! ======================================================================

791 792
#if _FABM_BOTTOM_INDEX_==-1
      _BEGIN_GLOBAL_HORIZONTAL_LOOP_
793 794 795 796 797
#  ifdef _FABM_VERTICAL_BOTTOM_TO_SURFACE_
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ == domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) depth _INDEX_GLOBAL_VERTICAL_(domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) = 1
#  else
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ == 1) depth _INDEX_GLOBAL_VERTICAL_(1) = 1
#  endif
798 799 800
      _END_GLOBAL_HORIZONTAL_LOOP_
#endif

Jorn Bruggeman's avatar
Jorn Bruggeman committed
801
      call start_test('fabm_do_bottom')
802
      bottom_loop_count = 0
803
      _BEGIN_OUTER_HORIZONTAL_LOOP_
804 805 806
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) then
#endif
807 808 809 810 811 812 813 814 815
         flux = 0
         sms_bt = 0
         call fabm_do_bottom(model _ARGUMENTS_HORIZONTAL_IN_,flux,sms_bt)
         do ivar=1,size(model%state_variables)
            call check_horizontal_slice_plus_1(flux,ivar,0.0_rk,-real(ivar+interior_state_offset,rk) _ARGUMENTS_HORIZONTAL_IN_)
         end do
         do ivar=1,size(model%bottom_state_variables)
            call check_horizontal_slice_plus_1(sms_bt,ivar,0.0_rk,-real(ivar+bottom_state_offset,rk) _ARGUMENTS_HORIZONTAL_IN_)
         end do
816 817 818
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
         endif
#endif
819
      _END_OUTER_HORIZONTAL_LOOP_
820
      call assert(bottom_loop_count == horizontal_count, 'fabm_do_surface', 'call count does not match number of (unmasked) horizontal points')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
821 822 823 824 825 826 827 828 829 830 831 832 833 834 835

      do ivar=1,size(model%horizontal_diagnostic_variables)
         if (model%horizontal_diagnostic_variables(ivar)%save .and. model%horizontal_diagnostic_variables(ivar)%target%source == source_do_bottom) then
            pdata_hz => fabm_get_horizontal_diagnostic_data(model, ivar)
            call check_horizontal(pdata_hz, model%horizontal_diagnostic_variables(ivar)%missing_value, -model%horizontal_diagnostic_variables(ivar)%missing_value)
         end if
      end do

      call report_test_result()

      ! ======================================================================
      ! Column-based operators
      ! ======================================================================

      call start_test('fabm_get_light')
836
      column_loop_count = 0
Jorn Bruggeman's avatar
Jorn Bruggeman committed
837 838
      _BEGIN_GLOBAL_HORIZONTAL_LOOP_
#ifdef _FABM_DEPTH_DIMENSION_INDEX_
839 840 841
#  if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_)
         ! No mask but non-constant bottom index. We need to skip everything below bottom
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) &
842
            call fabm_get_light(model,_IMIN_,_IMAX_ _ARG_VERTICAL_FIXED_LOCATION_)
843
#  else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
844
         call fabm_get_light(model,1,domain_extent(_FABM_DEPTH_DIMENSION_INDEX_) _ARG_VERTICAL_FIXED_LOCATION_)
845
#  endif
Jorn Bruggeman's avatar
Jorn Bruggeman committed
846
#else
847
         call fabm_get_light(model _ARGUMENTS_HORIZONTAL_IN_)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
848 849
#endif
      _END_GLOBAL_HORIZONTAL_LOOP_
850
      call assert(column_loop_count == interior_count, 'fabm_get_light', 'call count does not match number of (unmasked) interior points')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
851 852 853 854 855 856 857

      do ivar=1,size(model%diagnostic_variables)
         if (model%diagnostic_variables(ivar)%save .and. model%diagnostic_variables(ivar)%target%source==source_do_column) then
            pdata => fabm_get_interior_diagnostic_data(model, ivar)
            call check_interior(pdata, model%diagnostic_variables(ivar)%missing_value, -model%diagnostic_variables(ivar)%missing_value)
         end if
      end do
858 859 860 861 862 863 864 865

      do ivar=1,size(model%horizontal_diagnostic_variables)
         if (model%horizontal_diagnostic_variables(ivar)%save .and. model%horizontal_diagnostic_variables(ivar)%target%source == source_do_column) then
            pdata_hz => fabm_get_horizontal_diagnostic_data(model, ivar)
            call check_horizontal(pdata_hz, model%horizontal_diagnostic_variables(ivar)%missing_value, -model%horizontal_diagnostic_variables(ivar)%missing_value)
         end if
      end do

Jorn Bruggeman's avatar
Jorn Bruggeman committed
866 867 868 869 870 871 872
      call report_test_result()

      ! ======================================================================
      ! Retrieve vertical velocities (sinking, floating, active movement).
      ! ======================================================================

      call start_test('fabm_get_vertical_movement')
873
      vertical_movement_loop_count = 0
874
      _BEGIN_OUTER_INTERIOR_LOOP_
875 876 877 878 879
#if _FABM_BOTTOM_INDEX_==-1 && !defined(_HAS_MASK_) && _FABM_VECTORIZED_DIMENSION_INDEX_==_FABM_DEPTH_DIMENSION_INDEX_ && defined(_FABM_DEPTH_DIMENSION_INDEX_)
         ! We are looping over depth, but as we have a non-constant bottom index (yet no mask), we need to skip everything below bottom
         if (bottom_index _INDEX_HORIZONTAL_LOCATION_ >= 1 .and. bottom_index _INDEX_HORIZONTAL_LOCATION_ <= domain_extent(_FABM_DEPTH_DIMENSION_INDEX_)) &
            call fabm_get_vertical_movement(model,_IMIN_,_IMAX_ _ARG_INTERIOR_FIXED_LOCATION_,w(_IMIN_:_IMAX_,:))
#else
880
         call fabm_get_vertical_movement(model _ARGUMENTS_INTERIOR_IN_,w)
881
#endif
882
         do ivar=1,size(model%state_variables)
883 884 885 886 887
            if (mod(ivar, 2) == 0) then
               call check_interior_slice_plus_1(w,ivar,0.0_rk,real(ivar+interior_state_offset,rk) _ARGUMENTS_INTERIOR_IN_)
            else
               call check_interior_slice_plus_1(w,ivar,0.0_rk,-real(ivar+interior_state_offset,rk) _ARGUMENTS_INTERIOR_IN_)
            end if
888 889
         end do
      _END_OUTER_INTERIOR_LOOP_
890
      call assert(vertical_movement_loop_count == interior_count, 'fabm_get_vertical_movement', 'call count does not match number of (unmasked) interior points')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
891 892
      call report_test_result()

893 894 895 896 897 898
      ! ======================================================================
      ! Check state with valid state
      ! ======================================================================

      call start_test('fabm_check_state')
      _BEGIN_OUTER_INTERIOR_LOOP_
899 900
         call fabm_check_state(model _ARGUMENTS_INTERIOR_IN_, .true., valid)
         if (.not. valid) call driver%fatal_error('fabm_check_state', 'state is reported as invalid')
901 902 903 904 905
      _END_OUTER_INTERIOR_LOOP_
      call report_test_result()

      call start_test('fabm_check_surface_state')
      _BEGIN_OUTER_HORIZONTAL_LOOP_
906 907
         call fabm_check_surface_state(model _ARGUMENTS_HORIZONTAL_IN_, .true., valid)
         if (.not. valid) call driver%fatal_error('fabm_check_surface_state', 'state is reported as invalid')
908 909 910 911 912
      _END_OUTER_HORIZONTAL_LOOP_
      call report_test_result()

      call start_test('fabm_check_bottom_state')
      _BEGIN_OUTER_HORIZONTAL_LOOP_
913 914
         call fabm_check_bottom_state(model _ARGUMENTS_HORIZONTAL_IN_, .true., valid)
         if (.not. valid) call driver%fatal_error('fabm_check_bottom_state', 'state is reported as invalid')
915 916 917 918 919 920 921 922 923
      _END_OUTER_HORIZONTAL_LOOP_
      call report_test_result()

      ! ======================================================================
      ! Check state with state below minimum
      ! ======================================================================

      ! Now destroy the state by setting all values to below the minimum
      do ivar=1,size(model%state_variables)
924 925
        interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar) = model%state_variables(ivar)%minimum - abs(model%state_variables(ivar)%minimum) - 1
        call apply_mask_3d(interior_state(_PREARG_LOCATION_DIMENSIONS_ ivar), model%state_variables(ivar)%missing_value)
926 927
      end do
      do ivar=1,size(model%bottom_state_variables)
928 929
        bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar) = model%bottom_state_variables(ivar)%minimum - abs(model%bottom_state_variables(ivar)%minimum) - 1
        call apply_mask_2d(bottom_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar), model%bottom_state_variables(ivar)%missing_value)
930 931
      end do
      do ivar=1,size(model%surface_state_variables)
932 933
        surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar) = model%surface_state_variables(ivar)%minimum - abs(model%surface_state_variables(ivar)%minimum) - 1
        call apply_mask_2d(surface_state(_PREARG_HORIZONTAL_LOCATION_DIMENSIONS_ ivar), model%surface_state_variables(ivar)%missing_value)
934 935 936 937
      end do

      call start_test('fabm_check_state < min')
      _BEGIN_OUTER_INTERIOR_LOOP_
938 939 940
         call fabm_check_state(model _ARGUMENTS_INTERIOR_IN_, .true., valid)
#ifdef _HAS_MASK_
#  ifdef _FABM_HORIZONTAL_MASK_
941
         call assert(valid .neqv. any(_IS_UNMASKED_(mask_hz _INDEX_GLOBAL_HORIZONTAL_(_START_:_STOP_))), 'fabm_check_state', 'invalid result')
942
#  else
943
         call assert(valid .neqv. any(_IS_UNMASKED_(mask _INDEX_GLOBAL_INTERIOR_(_START_:_STOP_))), 'fabm_check_state'