python_fabm.F90 33.5 KB
Newer Older
Jorn Bruggeman's avatar
Jorn Bruggeman committed
1
#include "fabm_driver.h"
2 3 4

module fabm_python

5
   use iso_c_binding, only: c_double, c_int, c_char, C_NULL_CHAR, c_f_pointer, c_loc, c_ptr, c_null_ptr
6 7 8

   !DIR$ ATTRIBUTES DLLEXPORT :: STATE_VARIABLE,DIAGNOSTIC_VARIABLE,CONSERVED_QUANTITY

Jorn Bruggeman's avatar
Jorn Bruggeman committed
9
   use fabm, only: type_fabm_model, type_fabm_variable, fabm_get_version, status_start_done, fabm_create_model
10
   use fabm_types, only:rk => rke,attribute_length,type_model_list_node,type_base_model, &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
11
                        factory,type_link,type_link_list,type_internal_variable
12
   use fabm_driver, only: type_base_driver, driver, fatal_error
13
   use fabm_properties
Jorn Bruggeman's avatar
Jorn Bruggeman committed
14
   use fabm_python_helper
15
   use fabm_c_helper
Jorn Bruggeman's avatar
Jorn Bruggeman committed
16 17

   implicit none
18

Jorn Bruggeman's avatar
Jorn Bruggeman committed
19 20
   public

21 22 23 24 25 26
   integer, parameter :: INTERIOR_STATE_VARIABLE        = 1
   integer, parameter :: SURFACE_STATE_VARIABLE         = 2
   integer, parameter :: BOTTOM_STATE_VARIABLE          = 3
   integer, parameter :: INTERIOR_DIAGNOSTIC_VARIABLE   = 4
   integer, parameter :: HORIZONTAL_DIAGNOSTIC_VARIABLE = 5
   integer, parameter :: CONSERVED_QUANTITY             = 6
27

28
   logical, save :: error_occurred = .false.
29
   character(len=:), allocatable, save :: error_message
Jorn Bruggeman's avatar
Jorn Bruggeman committed
30

31 32
   type,extends(type_base_driver) :: type_python_driver
   contains
33 34
      procedure :: fatal_error => python_driver_fatal_error
      procedure :: log_message => python_driver_log_message
35
   end type
Jorn Bruggeman's avatar
Jorn Bruggeman committed
36

37 38 39 40 41 42 43 44 45
   type type_model_wrapper
      class (type_fabm_model), pointer :: p => null()
      character(len=1024), dimension(:), allocatable :: environment_names,environment_units
      integer :: index_column_depth
      type (type_link_list) :: coupling_link_list
      real(c_double),pointer :: column_depth => null()
      type (type_property_dictionary) :: forced_parameters,forced_couplings
   end type

46
contains
Jorn Bruggeman's avatar
Jorn Bruggeman committed
47

48
   subroutine get_version(length, version_string) bind(c)
49
      !DIR$ ATTRIBUTES DLLEXPORT :: get_version
50 51
      integer(c_int), value, intent(in) :: length
      character(kind=c_char)            :: version_string(length)
52

53
      character(len=length-1) :: string
54

55 56
      call fabm_get_version(string)
      call copy_to_c_string(string, version_string)
57 58
   end subroutine get_version

59 60 61 62
   function create_model(path) result(ptr) bind(c)
      !DIR$ ATTRIBUTES DLLEXPORT :: create_model
      character(kind=c_char), target, intent(in)  :: path(*)
      type(c_ptr)                                 :: ptr
63

64 65 66
      type (type_model_wrapper),       pointer :: model
      character(len=attribute_length), pointer :: ppath
      class (type_property),           pointer :: property
67

Jorn Bruggeman's avatar
Jorn Bruggeman committed
68
      ! Initialize driver object used by FABM for logging/error reporting.
69
      if (.not. associated(driver)) allocate(type_python_driver::driver)
70

Jorn Bruggeman's avatar
Jorn Bruggeman committed
71
      ! Build FABM model tree (configuration will be read from file specified as argument).
72
      allocate(model)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
73
      call c_f_pointer(c_loc(path), ppath)
74
      model%p => fabm_create_model(path=ppath(:index(ppath, C_NULL_CHAR) - 1), parameters=model%forced_parameters)
75

76
      ! Get a list of all parameters that had an explicit value specified.
77
      property => model%p%root%parameters%first
78
      do while (associated(property))
79 80
         if (.not. model%p%root%parameters%missing%contains(property%name)) &
            call model%forced_parameters%set_property(property)
81 82
         property => property%next
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
83

84
      ! Get a list of all active couplings.
85
      property => model%p%root%couplings%first
86
      do while (associated(property))
87 88
         if (.not. model%p%root%couplings%missing%contains(property%name)) &
            call model%forced_couplings%set_property(property)
89 90 91
         property => property%next
      end do

Jorn Bruggeman's avatar
Jorn Bruggeman committed
92
      ! Send information on spatial domain to FABM (this also allocates memory for diagnostics)
93
      call model%p%set_domain(1._rk)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
94

95
      ! Retrieve arrays to hold values for environmental variables and corresponding metadata.
96 97 98
      call get_environment_metadata(model%p, model%environment_names, model%environment_units, model%index_column_depth)

      call get_couplings(model%p, model%coupling_link_list)
99

100 101
      ptr = c_loc(model)
   end function create_model
102

103 104 105 106
   subroutine reinitialize(model)
      type (type_model_wrapper), intent(inout) :: model

      class (type_fabm_model),     pointer :: newmodel
107 108 109 110
      type (type_model_list_node), pointer :: node
      class (type_base_model),     pointer :: childmodel
      class (type_property),       pointer :: property,next

Jorn Bruggeman's avatar
Jorn Bruggeman committed
111
      ! Create new model object.
112 113 114
      allocate(newmodel)

      ! Transfer forced parameters to root of the model.
115 116
      call newmodel%root%parameters%update(model%forced_parameters)
      call newmodel%root%couplings%update(model%forced_couplings)
117 118

      ! Re-create original models
119
      node => model%p%root%children%first
120
      do while (associated(node))
121
         if (node%model%user_created) then
122
            call factory%create(node%model%type_name, childmodel)
123
            childmodel%user_created = .true.
124
            call newmodel%root%add_child(childmodel, node%model%name, node%model%long_name, configunit=-1)
125 126 127 128 129
         end if
         node => node%next
      end do

      ! Clean up old model
130 131
      call finalize(model)
      model%p => newmodel
132 133

      ! Initialize new model
134
      call model%p%initialize()
135 136

      ! Removed unused forced parameters from root model.
137
      property => model%p%root%parameters%first
138
      do while (associated(property))
139
         if (.not. model%p%root%parameters%retrieved%contains(property%name)) then
140
            next => property%next
141
            call model%p%root%parameters%delete(property%name)
142 143 144 145 146 147 148
            property => next
         else
            property => property%next
         end if
      end do

      ! Send information on spatial domain to FABM (this also allocates memory for diagnostics)
149
      call model%p%set_domain(1._rk)
150 151

      ! Retrieve arrays to hold values for environmental variables and corresponding metadata.
152 153
      call get_environment_metadata(model%p, model%environment_names, model%environment_units, model%index_column_depth)
      model%column_depth => null()
154

155
      call get_couplings(model%p, model%coupling_link_list)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
156
   end subroutine reinitialize
157

158 159 160 161 162 163 164 165 166
   subroutine start(pmodel) bind(c)
      !DIR$ ATTRIBUTES DLLEXPORT :: start
      type (c_ptr), intent(in), value :: pmodel

      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      call model%p%start()
   end subroutine start
167

168 169 170 171 172
   integer(c_int) function get_error_state() bind(c)
      !DIR$ ATTRIBUTES DLLEXPORT :: get_error_state
      get_error_state = logical2int(error_occurred)
   end function get_error_state

173 174 175 176 177 178 179
   subroutine get_error(length, message) bind(c)
      !DIR$ ATTRIBUTES DLLEXPORT :: get_error
      integer(c_int),        intent(in), value             :: length
      character(kind=c_char),intent(out),dimension(length) :: message
      call copy_to_c_string(error_message, message)
   end subroutine get_error

180
   integer(c_int) function model_count(pmodel) bind(c)
181
      !DIR$ ATTRIBUTES DLLEXPORT :: model_count
182
      type (c_ptr), intent(in), value :: pmodel
183

184
      type (type_model_wrapper),   pointer :: model
185 186
      type (type_model_list_node), pointer :: node

187
      call c_f_pointer(pmodel, model)
188
      model_count = 0
189
      node => model%p%root%children%first
190 191 192
      do while (associated(node))
         model_count = model_count + 1
         node => node%next
Jorn Bruggeman's avatar
Jorn Bruggeman committed
193
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
194
   end function model_count
195

196 197
   subroutine get_counts(pmodel, nstate_interior, nstate_surface, nstate_bottom, ndiagnostic_interior, ndiagnostic_horizontal, &
      nconserved, ndependencies, nparameters, ncouplings) bind(c)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
198
      !DIR$ ATTRIBUTES DLLEXPORT :: get_counts
199 200 201 202 203 204 205 206
      type (c_ptr),   intent(in), value :: pmodel
      integer(c_int), intent(out)       :: nstate_interior, nstate_surface, nstate_bottom
      integer(c_int), intent(out)       :: ndiagnostic_interior, ndiagnostic_horizontal
      integer(c_int), intent(out)       :: nconserved, ndependencies, nparameters, ncouplings

      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
207
      nstate_interior = size(model%p%interior_state_variables)
208 209
      nstate_surface = size(model%p%surface_state_variables)
      nstate_bottom = size(model%p%bottom_state_variables)
210
      ndiagnostic_interior = size(model%p%interior_diagnostic_variables)
211 212 213 214 215
      ndiagnostic_horizontal = size(model%p%horizontal_diagnostic_variables)
      nconserved = size(model%p%conserved_quantities)
      ndependencies = size(model%environment_names)
      nparameters = model%p%root%parameters%size()
      ncouplings = model%coupling_link_list%count()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
216
   end subroutine get_counts
217

218
   subroutine get_variable_metadata(pmodel, category, index, length, name, units, long_name, path) bind(c)
219
      !DIR$ ATTRIBUTES DLLEXPORT :: get_variable_metadata
220
      type (c_ptr),           intent(in), value              :: pmodel
221 222
      integer(c_int),         intent(in), value              :: category, index, length
      character(kind=c_char), intent(out), dimension(length) :: name, units, long_name, path
223

224 225
      type (type_model_wrapper),  pointer :: model
      class (type_fabm_variable), pointer :: variable
226 227

      call c_f_pointer(pmodel, model)
228 229 230

      ! Get a pointer to the target variable
      select case (category)
231
      case (INTERIOR_STATE_VARIABLE)
232
         variable => model%p%interior_state_variables(index)
233
      case (SURFACE_STATE_VARIABLE)
234
         variable => model%p%surface_state_variables(index)
235
      case (BOTTOM_STATE_VARIABLE)
236
         variable => model%p%bottom_state_variables(index)
237
      case (INTERIOR_DIAGNOSTIC_VARIABLE)
238
         variable => model%p%interior_diagnostic_variables(index)
239
      case (HORIZONTAL_DIAGNOSTIC_VARIABLE)
240
         variable => model%p%horizontal_diagnostic_variables(index)
241
      case (CONSERVED_QUANTITY)
242
         variable => model%p%conserved_quantities(index)
243
      end select
244 245 246 247
      call copy_to_c_string(variable%name,            name)
      call copy_to_c_string(variable%units,           units)
      call copy_to_c_string(variable%local_long_name, long_name)
      call copy_to_c_string(variable%path,            path)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
248
   end subroutine get_variable_metadata
Jorn Bruggeman's avatar
Jorn Bruggeman committed
249

250
   function get_variable(pmodel, category, index) bind(c) result(pvariable)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
251
      !DIR$ ATTRIBUTES DLLEXPORT :: get_variable
252
      type (c_ptr),   intent(in), value :: pmodel
253 254
      integer(c_int), intent(in), value :: category, index
      type (c_ptr)                      :: pvariable
255

256 257 258 259
      type (type_model_wrapper),     pointer :: model
      type (type_internal_variable), pointer :: variable

      call c_f_pointer(pmodel, model)
260 261 262

      ! Get a pointer to the target variable
      select case (category)
263
      case (INTERIOR_STATE_VARIABLE)
264
         variable => model%p%interior_state_variables(index)%target
265
      case (SURFACE_STATE_VARIABLE)
266
         variable => model%p%surface_state_variables(index)%target
267
      case (BOTTOM_STATE_VARIABLE)
268
         variable => model%p%bottom_state_variables(index)%target
269
      case (INTERIOR_DIAGNOSTIC_VARIABLE)
270
         variable => model%p%interior_diagnostic_variables(index)%target
271
      case (HORIZONTAL_DIAGNOSTIC_VARIABLE)
272
         variable => model%p%horizontal_diagnostic_variables(index)%target
273
      case (CONSERVED_QUANTITY)
274
         variable => model%p%conserved_quantities(index)%target
275 276
      end select
      pvariable = c_loc(variable)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
277
   end function get_variable
278

279
   subroutine get_parameter_metadata(pmodel, index, length, name, units, long_name, typecode, has_default) bind(c)
280
      !DIR$ ATTRIBUTES DLLEXPORT :: get_parameter_metadata
281
      type (c_ptr),           intent(in), value              :: pmodel
282 283 284
      integer(c_int),         intent(in), value              :: index, length
      character(kind=c_char), intent(out), dimension(length) :: name, units, long_name
      integer(c_int),         intent(out)                    :: typecode, has_default
285

286 287 288 289 290
      type (type_model_wrapper), pointer :: model
      integer                            :: i
      class (type_property),     pointer :: property

      call c_f_pointer(pmodel, model)
291 292

      i = 1
293
      property => model%p%root%parameters%first
294
      do while (associated(property))
295
         if (index == i) exit
296 297
         property => property%next
         i = i + 1
Jorn Bruggeman's avatar
Jorn Bruggeman committed
298 299
      end do

300 301 302
      call copy_to_c_string(property%name, name)
      call copy_to_c_string(property%units, units)
      call copy_to_c_string(property%long_name, long_name)
303
      typecode = property%typecode()
304
      has_default = logical2int(property%has_default)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
305
   end subroutine get_parameter_metadata
Jorn Bruggeman's avatar
Jorn Bruggeman committed
306

307
   subroutine get_dependency_metadata(pmodel, index, length, name, units) bind(c)
308
      !DIR$ ATTRIBUTES DLLEXPORT :: get_dependency_metadata
309
      type (c_ptr),           intent(in), value              :: pmodel
310 311
      integer(c_int),         intent(in), value              :: index, length
      character(kind=c_char), intent(out), dimension(length) :: name, units
Jorn Bruggeman's avatar
Jorn Bruggeman committed
312

313 314 315 316 317
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      call copy_to_c_string(model%environment_names(index), name)
      call copy_to_c_string(model%environment_units(index), units)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
318
   end subroutine get_dependency_metadata
319

320
   subroutine get_coupling(pmodel, index, slave, master) bind(c)
321
      !DIR$ ATTRIBUTES DLLEXPORT :: get_coupling
322
      type (c_ptr),   intent(in), value :: pmodel
323 324
      integer(c_int), intent(in), value :: index
      type (c_ptr),   intent(out)       :: slave, master
325

326
      type (type_model_wrapper), pointer :: model
327 328
      type (type_link),pointer :: link_slave
      integer                  :: i
Knut's avatar
Knut committed
329

330 331
      call c_f_pointer(pmodel, model)
      link_slave => model%coupling_link_list%first
332 333 334 335 336
      do i=2,index
         link_slave => link_slave%next
      end do
      slave = c_loc(link_slave%original)
      master = c_loc(link_slave%target)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
337
   end subroutine get_coupling
338

339
   function variable_get_suitable_masters(pmodel, pvariable) result(plist) bind(c)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
340
      !DIR$ ATTRIBUTES DLLEXPORT :: variable_get_suitable_masters
341
      type (c_ptr), intent(in), value :: pmodel
342 343
      type (c_ptr), intent(in), value :: pvariable
      type (c_ptr)                    :: plist
344

345 346 347
      type (type_model_wrapper),     pointer :: model
      type (type_internal_variable), pointer :: variable
      type (type_link_list),         pointer :: list
348

349
      call c_f_pointer(pmodel, model)
350
      call c_f_pointer(pvariable, variable)
351
      list => get_suitable_masters(model%p, variable)
352
      plist = c_loc(list)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
353
   end function variable_get_suitable_masters
354

355
   subroutine get_model_metadata(pmodel, name, length, long_name, user_created) bind(c)
356
      !DIR$ ATTRIBUTES DLLEXPORT :: get_model_metadata
357
      type (c_ptr),           intent(in), value  :: pmodel
358 359 360 361
      character(kind=c_char), intent(in), target :: name(*)
      integer(c_int),         intent(in), value  :: length
      character(kind=c_char), intent(out)        :: long_name(length)
      integer(c_int),         intent(out)        :: user_created
362

363 364 365
      type (type_model_wrapper),       pointer :: model
      character(len=attribute_length), pointer :: pname
      class (type_base_model),         pointer :: found_model
366

367
      call c_f_pointer(pmodel, model)
368
      call c_f_pointer(c_loc(name), pname)
369
      found_model => model%p%root%find_model(pname(:index(pname, C_NULL_CHAR) - 1))
370
      if (.not.associated(found_model)) call fatal_error('get_model_metadata', &
371 372
         'model "'//pname(:index(pname, C_NULL_CHAR) - 1) // '" not found.')
      call copy_to_c_string(found_model%long_name, long_name)
373
      user_created = logical2int(found_model%user_created)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
374
   end subroutine get_model_metadata
375

376
   subroutine link_dependency_data(pmodel, index, value) bind(c)
377
      !DIR$ ATTRIBUTES DLLEXPORT :: link_dependency_data
378
      type (c_ptr),   intent(in), value  :: pmodel
379 380
      integer(c_int), intent(in), value  :: index
      real(c_double), intent(in), target :: value
381

382 383 384 385 386 387 388
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      call model%p%link_interior_data(model%environment_names(index), value)
      call model%p%link_horizontal_data(model%environment_names(index), value)
      call model%p%link_scalar(model%environment_names(index), value)
      if (index == model%index_column_depth) model%column_depth => value
Jorn Bruggeman's avatar
Jorn Bruggeman committed
389
   end subroutine link_dependency_data
390

391
   subroutine link_interior_state_data(pmodel, index, value) bind(c)
392
      !DIR$ ATTRIBUTES DLLEXPORT :: link_interior_state_data
393
      type (c_ptr),   intent(in), value     :: pmodel
394 395
      integer(c_int), intent(in), value     :: index
      real(c_double), intent(inout), target :: value
Jorn Bruggeman's avatar
Jorn Bruggeman committed
396

397 398 399
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
400
      value = model%p%interior_state_variables(index)%initial_value
401
      call model%p%link_interior_state_data(index, value)
402
   end subroutine link_interior_state_data
Jorn Bruggeman's avatar
Jorn Bruggeman committed
403

404
   subroutine link_surface_state_data(pmodel, index, value) bind(c)
405
      !DIR$ ATTRIBUTES DLLEXPORT :: link_surface_state_data
406
      type (c_ptr),   intent(in), value     :: pmodel
407 408
      integer(c_int), intent(in), value     :: index
      real(c_double), intent(inout), target :: value
Jorn Bruggeman's avatar
Jorn Bruggeman committed
409

410 411 412 413 414
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      value = model%p%surface_state_variables(index)%initial_value
      call model%p%link_surface_state_data(index, value)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
415
   end subroutine link_surface_state_data
Knut's avatar
Knut committed
416

417
   subroutine link_bottom_state_data(pmodel, index, value) bind(c)
418
      !DIR$ ATTRIBUTES DLLEXPORT :: link_bottom_state_data
419
      type (c_ptr),   intent(in), value     :: pmodel
420 421
      integer(c_int), intent(in), value     :: index
      real(c_double), intent(inout), target :: value
422

423 424 425 426 427
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      value = model%p%bottom_state_variables(index)%initial_value
      call model%p%link_bottom_state_data(index, value)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
428
   end subroutine link_bottom_state_data
429

430
   subroutine get_rates(pmodel, t, dy, do_surface, do_bottom) bind(c)
431
      !DIR$ ATTRIBUTES DLLEXPORT :: get_rates
432
      type (c_ptr),           intent(in), value :: pmodel
433 434 435
      real(rk), value,        intent(in) :: t
      real(c_double), target, intent(in) :: dy(*)
      integer(c_int), value,  intent(in) :: do_surface, do_bottom
436

437
      type (type_model_wrapper), pointer :: model
438
      real(c_double), pointer :: dy_(:)
439

440
      call c_f_pointer(pmodel, model)
441
      if (model%p%status < status_start_done) then
442 443 444 445
         call fatal_error('get_rates', 'start has not been called yet.')
         return
      end if

Jorn Bruggeman's avatar
Jorn Bruggeman committed
446 447
      call c_f_pointer(c_loc(dy), dy_, (/size(model%p%interior_state_variables) &
         + size(model%p%surface_state_variables) + size(model%p%bottom_state_variables)/))
448

449
      if (t < 0) then
450
         call model%p%prepare_inputs()
451
      else
452
         call model%p%prepare_inputs(t)
453 454
      end if
      dy_ = 0.0_rk
455
      if (int2logical(do_surface)) call model%p%get_surface_sources(dy_(1:size(model%p%interior_state_variables)), &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
456 457
         dy_(size(model%p%interior_state_variables)+1:size(model%p%interior_state_variables) &
         + size(model%p%surface_state_variables)))
458 459
      if (int2logical(do_bottom)) call model%p%get_bottom_sources(dy_(1:size(model%p%interior_state_variables)), &
         dy_(size(model%p%interior_state_variables) + size(model%p%surface_state_variables) + 1:))
460
      if (int2logical(do_surface) .or. int2logical(do_bottom)) then
461 462
         if (.not.associated(model%column_depth)) call fatal_error('get_rates', &
            'Value for environmental dependency ' // trim(model%environment_names(model%index_column_depth)) // &
463
            ' must be provided if get_rates is called with the do_surface and/or do_bottom flags.')
464
         dy_(1:size(model%p%interior_state_variables)) = dy_(1:size(model%p%interior_state_variables)) / model%column_depth
465
      end if
466
      call model%p%get_interior_sources(dy_(1:size(model%p%interior_state_variables)))
467
      call model%p%finalize_outputs()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
468 469

      ! Compute rate of change in conserved quantities
470
      !call fabm_state_to_conserved_quantities(model,pelagic_rates,conserved_rates)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
471 472

      ! Normalize rate of change in conserved quantities to sum of absolute rates of change.
473 474
      !call fabm_state_to_conserved_quantities(model,abs(pelagic_rates),abs_conserved_rates)
      !where (abs_conserved_rates>0.0_rk) conserved_rates = conserved_rates/abs_conserved_rates
Jorn Bruggeman's avatar
Jorn Bruggeman committed
475
   end subroutine get_rates
476

477
   function check_state(pmodel, repair_) bind(c) result(valid_)
478
      !DIR$ ATTRIBUTES DLLEXPORT :: check_state
479 480 481
      type (c_ptr),         intent(in), value :: pmodel
      integer(c_int),value, intent(in)        :: repair_
      integer(c_int)                          :: valid_
482

483
      type (type_model_wrapper), pointer :: model
484 485
      logical :: repair, interior_valid, surface_valid, bottom_valid

486
      call c_f_pointer(pmodel, model)
487
      if (model%p%status < status_start_done) then
488 489 490 491
         call fatal_error('check_state', 'start has not been called yet.')
         return
      end if

492
      repair = int2logical(repair_)
493 494 495
      call model%p%check_interior_state(repair, interior_valid)
      call model%p%check_surface_state(repair, surface_valid)
      call model%p%check_bottom_state(repair, bottom_valid)
496 497 498
      valid_ = logical2int(interior_valid .and. surface_valid .and. bottom_valid)
   end function check_state

499
   subroutine integrate(pmodel, nt, ny, t_, y_ini_, y_, dt, do_surface, do_bottom) bind(c)
500
      !DIR$ ATTRIBUTES DLLEXPORT :: integrate
501
      type (c_ptr),  value, intent(in) :: pmodel
502 503 504 505 506
      integer(c_int),value, intent(in) :: nt, ny
      real(c_double),target,intent(in) :: t_(*), y_ini_(*), y_(*)
      real(c_double),value, intent(in) :: dt
      integer(c_int),value, intent(in) :: do_surface, do_bottom

507
      type (type_model_wrapper), pointer :: model
508 509 510 511
      real(c_double),pointer :: t(:), y_ini(:), y(:,:)
      integer                :: it
      real(rk)               :: t_cur
      real(rk), target       :: y_cur(ny)
512
      real(rk)               :: dy(ny)
513 514
      logical                :: surface, bottom

515
      call c_f_pointer(pmodel, model)
516
      if (model%p%status < status_start_done) then
517 518 519
         call fatal_error('integrate', 'start has not been called yet.')
         return
      end if
Jorn Bruggeman's avatar
Jorn Bruggeman committed
520 521
      if (ny /= size(model%p%interior_state_variables) + size(model%p%surface_state_variables) &
         + size(model%p%bottom_state_variables)) then
522 523 524
         call fatal_error('integrate', 'ny is wrong length')
         return
      end if
525 526 527 528 529 530 531

      call c_f_pointer(c_loc(t_), t, (/nt/))
      call c_f_pointer(c_loc(y_ini_), y_ini, (/ny/))
      call c_f_pointer(c_loc(y_), y, (/ny, nt/))

      surface = int2logical(do_surface)
      bottom = int2logical(do_bottom)
532 533
      if ((surface .or. bottom) .and. .not. associated(model%column_depth)) then
          call fatal_error('get_rates', &
534
            'Value for environmental dependency ' // trim(model%environment_names(model%index_column_depth)) // &
535
            ' must be provided if integrate is called with the do_surface and/or do_bottom flags.')
536
          return
537
      end if
538 539 540
      call model%p%link_all_interior_state_data(y_cur(1:size(model%p%interior_state_variables)))
      call model%p%link_all_surface_state_data(y_cur(size(model%p%interior_state_variables) + 1: &
         size(model%p%interior_state_variables) + size(model%p%surface_state_variables)))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
541 542
      call model%p%link_all_bottom_state_data(y_cur(size(model%p%interior_state_variables) &
         + size(model%p%surface_state_variables) + 1:))
543 544 545 546 547 548 549 550 551 552

      it = 1
      t_cur = t(1)
      y_cur = y_ini
      do while (it <= nt)
          if (t_cur >= t(it)) then
              y(:, it) = y_cur
              it = it + 1
          end if

553
          call model%p%prepare_inputs(t_cur)
554
          dy = 0.0_rk
555
          if (surface) call model%p%get_surface_sources(dy(1:size(model%p%interior_state_variables)), &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
556 557
             dy(size(model%p%interior_state_variables) + 1:size(model%p%interior_state_variables) &
             + size(model%p%surface_state_variables)))
558 559
          if (bottom) call model%p%get_bottom_sources(dy(1:size(model%p%interior_state_variables)), &
             dy(size(model%p%interior_state_variables) + size(model%p%surface_state_variables) + 1:))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
560 561
          if (surface .or. bottom) dy(1:size(model%p%interior_state_variables)) = dy(1:size(model%p%interior_state_variables)) &
             / model%column_depth
562
          call model%p%get_interior_sources(dy(1:size(model%p%interior_state_variables)))
563
          y_cur = y_cur + dt * dy * 86400
564 565 566 567
          t_cur = t_cur + dt
      end do
   end subroutine integrate

568
   subroutine get_interior_diagnostic_data(pmodel, index, ptr) bind(c)
569
      !DIR$ ATTRIBUTES DLLEXPORT :: get_interior_diagnostic_data
570
      type (c_ptr),   intent(in), value :: pmodel
571 572
      integer(c_int), intent(in), value :: index
      type(c_ptr),    intent(out)       :: ptr
573 574
      real(rk), pointer :: pvalue

575 576 577
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
578
      ptr = c_null_ptr
579
      pvalue => model%p%get_interior_diagnostic_data(index)
580
      if (associated(pvalue)) ptr = c_loc(pvalue)
581
   end subroutine get_interior_diagnostic_data
582

583
   subroutine get_horizontal_diagnostic_data(pmodel, index, ptr) bind(c)
584
      !DIR$ ATTRIBUTES DLLEXPORT :: get_horizontal_diagnostic_data
585
      type (c_ptr),   intent(in), value :: pmodel
586 587
      integer(c_int), intent(in), value :: index
      type(c_ptr),    intent(out)       :: ptr
588 589
      real(rk), pointer :: pvalue

590 591 592
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
593
      ptr = c_null_ptr
594
      pvalue => model%p%get_horizontal_diagnostic_data(index)
595
      if (associated(pvalue)) ptr = c_loc(pvalue)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
596
   end subroutine get_horizontal_diagnostic_data
Jorn Bruggeman's avatar
Jorn Bruggeman committed
597

598 599 600 601 602 603 604 605 606
   subroutine finalize(model)
      type (type_model_wrapper), intent(inout) :: model

      call model%p%finalize()
      if (allocated(model%environment_names)) deallocate(model%environment_names)
      if (allocated(model%environment_units)) deallocate(model%environment_units)
      call model%forced_parameters%finalize()
      call model%forced_couplings%finalize()
      deallocate(model%p)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
607
   end subroutine finalize
Jorn Bruggeman's avatar
Jorn Bruggeman committed
608

609
   subroutine reset_parameter(pmodel, index) bind(c)
610
      !DIR$ ATTRIBUTES DLLEXPORT :: reset_parameter
611
      type (c_ptr), value,   intent(in) :: pmodel
612 613
      integer(c_int), value, intent(in) :: index
      class (type_property), pointer    :: property
614

615 616 617 618
      type (type_model_wrapper), pointer :: model

      call c_f_pointer(pmodel, model)
      property => model%p%root%parameters%get_property(index)
619
      if (.not. associated(property)) return
620
      call model%forced_parameters%delete(property%name)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
621 622

      ! Re-initialize the model using updated parameter values
623
      call reinitialize(model)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
624
   end subroutine reset_parameter
Jorn Bruggeman's avatar
Jorn Bruggeman committed
625

626
   subroutine set_real_parameter(pmodel, name, value) bind(c)
627
      !DIR$ ATTRIBUTES DLLEXPORT :: set_real_parameter
628
      type (c_ptr),   value,          intent(in) :: pmodel
629
      character(kind=c_char), target, intent(in) :: name(*)
630
      real(c_double), value,          intent(in) :: value
631

632 633
      type (type_model_wrapper),       pointer :: model
      character(len=attribute_length), pointer :: pname
634

635
      call c_f_pointer(pmodel, model)
636
      call c_f_pointer(c_loc(name), pname)
Jorn Bruggeman's avatar