fabm0d.F90 35.1 KB
Newer Older
1
#include "fabm_driver.h"
2
#include "fabm_0d.h"
3 4 5 6 7 8 9 10 11 12 13 14 15
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: 0D independent driver for the Framework for Aquatic Biogeochemical Models (FABM)
!
! !INTERFACE:
   module fabm_0d
!
! !DESCRIPTION:
! TODO
!
! !USES:
   use time
16
   use input
17 18
   use eqstate,only:rho_feistel

19
   use fabm
20 21
   use fabm_driver
   use fabm_expressions
22

23 24 25
   use shared
   use output

26
   implicit none
27 28 29 30 31 32
   private
!
! !PUBLIC MEMBER FUNCTIONS:
   public init_run, time_loop, clean_up
!
! !DEFINED PARAMETERS:
Jorn Bruggeman's avatar
Jorn Bruggeman committed
33 34
   integer, parameter :: namlst=10, yaml_unit=23
   integer, parameter :: CENTER=0,SURFACE=1,BOTTOM=2
35 36 37 38
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
Jorn Bruggeman's avatar
Jorn Bruggeman committed
39 40
!  Run configuration file
   character(len=PATH_MAX) :: run_nml_file='run.nml'
41 42

!  FABM nml configuration file
Jorn Bruggeman's avatar
Jorn Bruggeman committed
43
   character(len=PATH_MAX) :: fabm_nml_file='fabm.nml'
44 45

!  FABM yaml configuration file
Jorn Bruggeman's avatar
Jorn Bruggeman committed
46
   character(len=PATH_MAX) :: fabm_yaml_file='fabm.yaml'
47

48
   ! Bio model info
49
   integer  :: ode_method
50
   logical  :: repair_state
51
   integer  :: swr_method
52
   logical  :: albedo_correction
53 54 55 56
   real(rk) :: cloud
   real(rk) :: par_fraction
   real(rk) :: par_background_extinction
   logical  :: apply_self_shading
57
   integer  :: model_type
58
   real(rk),allocatable :: current_rhs(:)
59

Jorn Bruggeman's avatar
Jorn Bruggeman committed
60 61
   real(rk), pointer :: bio_albedo, bio_extinction

Jorn Bruggeman's avatar
Jorn Bruggeman committed
62 63
   ! Shortcuts to number of state variables (interior, surface, bottom)
   integer :: n_int, n_sf, n_bt
64

65
   ! Environment
66
   real(rk),target :: current_depth,dens,decimal_yearday
67
   real(rk)        :: swr_sf,par_sf,par_bt,par_ct,extinction
68

69
   real(rk),allocatable :: expression_data(:)
70
   real(rk),allocatable :: totals0(:)
71

Jorn Bruggeman's avatar
Jorn Bruggeman committed
72 73 74
   type (type_scalar_input) :: mixing_rate
   type (type_scalar_input) :: mixed_layer_depth
   type (type_scalar_input), allocatable :: cc_deep(:)
75
   real(rk),allocatable,target :: w(:)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
76
   integer(timestepkind), save :: itime
77

Jorn Bruggeman's avatar
Jorn Bruggeman committed
78 79
   type (type_fabm_interior_variable_id), save :: id_dens, id_par
   logical                                     :: compute_density
80

Jorn Bruggeman's avatar
Jorn Bruggeman committed
81
   type,extends(type_base_driver) :: type_fabm0d_driver
82
   contains
Jorn Bruggeman's avatar
Jorn Bruggeman committed
83 84
      procedure :: fatal_error => fabm0d_driver_fatal_error
      procedure :: log_message => fabm0d_driver_log_message
85
   end type
86
!EOP
87 88 89 90
!-----------------------------------------------------------------------

   contains

Knut's avatar
Knut committed
91
#define _ODE_ZEROD_
92 93
#include "ode_solvers_template.F90"

94 95 96 97 98 99
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Parse the command line
!
! !INTERFACE:
100
   subroutine cmdline
101 102 103 104 105 106 107

!   character(len=*), parameter :: version = '1.0'
   character(len=32) :: arg
   integer :: i
!EOP
!-----------------------------------------------------------------------
!BOC
108 109
   i=1
   do while (i <= command_argument_count())
110
      call get_command_argument(i, arg)
Knut's avatar
Knut committed
111

112 113 114
      select case (arg)
#if 0
      case ('-v', '--version')
115
         print '(2a)', 'fabm0d version ', RELEASE
116 117 118 119 120
         stop
#endif
      case ('-h', '--help')
         call print_help()
         stop
121 122 123 124 125 126 127 128
      case ('-r', '--run_nml')
         i = i+1
         call get_command_argument(i, run_nml_file)
#if 0
      case ('-n', '--nml')
         i = i+1
         call get_command_argument(i, fabm_nml_file)
#endif
129
      case ('-y', '--yaml')
130 131
         i = i+1
         call get_command_argument(i, fabm_yaml_file)
132 133 134 135 136
      case default
         print '(a,a,/)', 'Unrecognized command-line option: ', arg
         call print_help()
         stop
      end select
137
      i = i+1
138
   end do
139 140 141 142 143
#if 0
   print '(a)', trim(run_nml_file)
   print '(a)', trim(fabm_nml_file)
   print '(a)', trim(fabm_yaml_file)
#endif
144 145 146 147 148 149 150 151 152 153 154

   contains

   subroutine print_help()
      print '(a)', 'usage: fabm0d [OPTIONS]'
      print '(a)', ''
      print '(a)', 'Without further options, fabm0d run using default input filenames.'
      print '(a)', ''
      print '(a)', 'fabm0d options:'
      print '(a)', ''
      print '(a)', '  -h, --help        print usage information and exit'
Jorn Bruggeman's avatar
Jorn Bruggeman committed
155 156
      print '(a)', '  -r, --run_nml     namelist file with simualtion settings - default run.nml' 
      print '(a)', '  -y, --yaml file   yaml-formatted file FABM configuration - default fabm.yaml'
157
      print '(a)', ''
158 159 160 161
   end subroutine print_help

   end subroutine  cmdline

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Initialise the model
!
! !INTERFACE:
   subroutine init_run()
!
! !DESCRIPTION:
!  This internal routine triggers the initialization of the model.
!  The first section reads the namelists of {\tt run.nml} with
!  the user specifications. Then, one by one each of the modules are
!  initialised.
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
!
! !LOCAL VARIABLES:
181
   character(len=PATH_MAX)   :: env_file
Jorn Bruggeman's avatar
Jorn Bruggeman committed
182
   real(rk)                  :: depth, dt
183
   real(rk),parameter        :: invalid_latitude = -100._rk,invalid_longitude = -400.0_rk
184
   logical                   :: file_exists
Jorn Bruggeman's avatar
Jorn Bruggeman committed
185
   integer                   :: ios
186

187
   namelist /model_setup/ title,start,stop,dt,ode_method,repair_state,model_type
188
   namelist /environment/ env_file,swr_method,albedo_correction, &
189 190
                          latitude,longitude,cloud,par_fraction, &
                          depth,par_background_extinction,apply_self_shading
191
!EOP
192 193
!-----------------------------------------------------------------------
!BOC
194 195

   ! Make FABM use our custom logger/error reporter
Jorn Bruggeman's avatar
Jorn Bruggeman committed
196
   allocate(type_fabm0d_driver::driver)
197

198 199
   call cmdline

200
   LEVEL1 'init_run'
201 202 203
   STDERR LINE

   ! Open the namelist file.
204
   LEVEL2 'reading model setup namelists from ',trim(run_nml_file)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
205 206
   open(namlst,file=run_nml_file,status='old',action='read',iostat=ios)
   if (ios/=0) call fatal_error('init_run','I could not open '//trim(run_nml_file)//' for reading.')
207

208
   ! Initialize environment
Jorn Bruggeman's avatar
Jorn Bruggeman committed
209 210 211
   temp%value = 0.0_rk
   salt%value = 0.0_rk
   light%value = 0.0_rk
212
   dens = 0.0_rk
213
   par_sf = 0.0_rk
214 215
   par_bt = 0.0_rk
   par_ct = 0.0_rk
216
   decimal_yearday = 0.0_rk
217
   model_type = 0
218 219 220 221 222 223 224

   ! Read all namelists
   title = ''
   start = ''
   stop = ''
   dt = 0.0_rk
   ode_method = 1
225
   repair_state = .false.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
226 227
   read(namlst,nml=model_setup,iostat=ios)
   if (ios/=0) call fatal_error('init_run','I could not read the "model_setup" namelist from '//trim(run_nml_file)//'.')
Knut's avatar
Knut committed
228

229 230 231
   ! Read environment namelist
   env_file = ''
   swr_method = 0
232
   albedo_correction = .true.
233 234
   latitude = invalid_latitude
   longitude = invalid_longitude
235 236 237 238 239
   cloud = 0.0_rk
   par_fraction = 1.0_rk
   depth = -1.0_rk
   par_background_extinction = 0.0_rk
   apply_self_shading = .true.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
240 241
   read(namlst,nml=environment,iostat=ios)
   if (ios/=0) call fatal_error('init_run','I could not read the "environment" namelist from '//trim(run_nml_file)//'.')
242

243
   compute_conserved_quantities = .false.
244
   call configure_output(namlst)
245

246 247 248
   ! Close the namelist file.
   close (namlst)

Jorn Bruggeman's avatar
Jorn Bruggeman committed
249 250 251 252
   if (start=='')  call fatal_error('init_run',trim(run_nml_file)//': start time "start" must be set in "model_setup" namelist.')
   if (stop=='')   call fatal_error('init_run',trim(run_nml_file)//': stop time "stop" must be set in "model_setup" namelist.')
   if (dt<=0.0_rk) call fatal_error('init_run',trim(run_nml_file)//': time step "dt" must be set to a positive value in "model_setup" namelist.')
   if (env_file=='') call fatal_error('init_run',trim(run_nml_file)//': "env_file" must be set to a valid file path in "environment" namelist.')
253
   if (latitude/=invalid_latitude.and.(latitude<-90._rk.or.latitude>90._rk)) &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
254
      call fatal_error('init_run',trim(run_nml_file)//': latitude must lie between -90 and 90.')
255
   if (longitude/=invalid_longitude.and.(longitude<-360._rk.or.longitude>360._rk)) &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
256
      call fatal_error('init_run',trim(run_nml_file)//': longitude must lie between -360 and 360.')
257

258
   ! Make sure depth has been provided.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
259
   if (depth<=0.0_rk) call fatal_error('init_run',trim(run_nml_file)//': &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
260
      &a positive value for "depth" must be provided in "environment" namelist.')
261
   column_depth = depth ! Provided depth is the column depth. The modelled biogeochemistry will be positioned at half this depth.
262
   call update_depth(CENTER)
Knut's avatar
Knut committed
263

264 265
   ! If longitude and latitude are used, make sure they have been provided and are valid.
   if (swr_method==0) then
Jorn Bruggeman's avatar
Jorn Bruggeman committed
266
      if (latitude==invalid_latitude) call fatal_error('init_run',trim(run_nml_file)//': &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
267
         &a valid value for "latitude" must be provided in "environment" if "swr_method" is 0.')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
268
      if (longitude==invalid_longitude) call fatal_error('init_run',trim(run_nml_file)//': &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
269
         &a valid value for "longitude" must be provided in "environment" if "swr_method" is 0.')
270
   end if
Knut's avatar
Knut committed
271

272 273 274 275 276 277 278
   ! Configure the time module to use actual start and stop dates.
   timefmt = 2

   ! Transfer the time step to the time module.
   timestep = dt

   ! Write information for this run to the console.
279
   LEVEL2 'Simulation: '//trim(title)
280 281 282 283 284 285 286
   select case (swr_method)
      case (0)
         LEVEL2 'Surface photosynthetically active radiation will be calculated from time,'
         LEVEL2 'cloud cover, and the simulated location at (lat,long)'
         LEVEL2 latitude,longitude
         LEVEL2 'Local PAR will be calculated from the surface value,'
         LEVEL2 'depth, and light extinction coefficient.'
287
         LEVEL2 'albedo_correction =',albedo_correction
288 289 290 291 292 293 294
      case (1)
         LEVEL2 'Surface photosynthetically active radiation (PAR) is provided as input.'
         LEVEL2 'Local PAR will be calculated from the surface value,'
         LEVEL2 'depth, and light extinction coefficient.'
      case (2)
         LEVEL2 'Local photosynthetically active radiation is provided as input.'
   end select
295 296 297 298 299 300 301

   LEVEL2 'initializing modules....'

   ! Initialize the time module.
   call init_time(MinN,MaxN)

   ! Open the file with observations of the local environment.
302 303
   LEVEL1 'init environment'
   LEVEL2 'reading local environment data from:'
304
   LEVEL3 trim(env_file)
305
   call init_input()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
306 307 308
   call light%configure(method=2, path=env_file, index=1, name='shortwave radiation')
   call temp%configure(method=2, path=env_file, index=2, name='temperature')
   call salt%configure(method=2, path=env_file, index=3, name='salinity')
309 310 311
   call register_input(light)
   call register_input(temp)
   call register_input(salt)
312

313 314 315 316
   ! Build FABM model tree. Use 'fabm_yaml_file' if available, otherwise fall back to fabm.nml.
   LEVEL1 'initialize FABM'
   LEVEL2 'reading configuration from:'
   inquire(file=trim(fabm_yaml_file),exist=file_exists)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
317 318 319
   if (.not. file_exists) call fatal_error('init_run','can not find '//trim(fabm_yaml_file)//'.')
   LEVEL3 trim(fabm_yaml_file)
   model => fabm_create_model(path=trim(fabm_yaml_file))
320

321
   ! Shortcuts to the number of state variables.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
322
   n_int = size(model%interior_state_variables)
323 324 325 326 327
   n_sf  = size(model%surface_state_variables)
   n_bt  = size(model%bottom_state_variables)

   allocate(cc(n_int+n_bt+n_sf))

328 329 330 331
   if (model_type==1) then
      call driver%log_message('The model type is set to mixed layer model (model_type = 1).')
      call driver%log_message('Therefore, bottom-associated processes will be deactivated.')
      allocate(cc_deep(n_int))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
332 333
      cc_deep(:)%value = 0.0_rk
      mixing_rate%value = 0.0_rk
334 335 336
      allocate(w(n_int))
   end if

Jorn Bruggeman's avatar
Jorn Bruggeman committed
337
   ! Allocate memory to hold totals of conserved quantities
338 339 340 341
   allocate(totals0             (size(model%conserved_quantities)))  ! at initial time (depth-integrated, interior + interfaces)
   allocate(totals              (size(model%conserved_quantities)))  ! at current time (depth-explicit, interior only)
   allocate(int_change_in_totals(size(model%conserved_quantities)))  ! change since start of simulation (depth-integrated, interior + interfaces)

Jorn Bruggeman's avatar
Jorn Bruggeman committed
342
   call register_output_fields()
343

344
   ! Send information on spatial domain to FABM (this also allocates memory for diagnostics)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
345
   call model%set_domain(seconds_per_time_unit=timestep)
346

347 348
   ! Create state variable vector, using the initial values specified by the model,
   ! and link state data to FABM.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
349 350 351
   call model%link_all_interior_state_data(cc(1:n_int))
   call model%link_all_bottom_state_data  (cc(n_int+1:n_int+n_bt))
   call model%link_all_surface_state_data (cc(n_int+n_bt+1:n_int+n_bt+n_sf))
352

Jorn Bruggeman's avatar
Jorn Bruggeman committed
353
   id_dens = model%get_interior_variable_id(fabm_standard_variables%density)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
354
   compute_density = model%variable_needs_values(id_dens)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
355
   if (compute_density) call model%link_interior_data(id_dens,dens)
356

Jorn Bruggeman's avatar
Jorn Bruggeman committed
357
   id_par = model%get_interior_variable_id(fabm_standard_variables%downwelling_photosynthetic_radiative_flux)
358

359
   ! Link environmental data to FABM
Jorn Bruggeman's avatar
Jorn Bruggeman committed
360 361
   call model%link_interior_data(fabm_standard_variables%temperature,temp%value)
   call model%link_interior_data(fabm_standard_variables%practical_salinity,salt%value)
362
   if (model%variable_needs_values(id_par)) call model%link_interior_data(id_par,par)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
363 364 365 366 367 368 369 370 371 372 373 374
   call model%link_interior_data(fabm_standard_variables%pressure,current_depth)
   call model%link_interior_data(fabm_standard_variables%cell_thickness,column_depth)
   call model%link_interior_data(fabm_standard_variables%depth,current_depth)
   call model%link_interior_data(fabm_standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux,extinction)
   call model%link_horizontal_data(fabm_standard_variables%surface_downwelling_photosynthetic_radiative_flux,par_sf)
   call model%link_horizontal_data(fabm_standard_variables%surface_downwelling_shortwave_flux,swr_sf)
   call model%link_horizontal_data(fabm_standard_variables%cloud_area_fraction,cloud)
   call model%link_horizontal_data(fabm_standard_variables%bottom_depth,column_depth)
   call model%link_horizontal_data(fabm_standard_variables%bottom_depth_below_geoid,column_depth)
   if (latitude /=invalid_latitude ) call model%link_horizontal_data(fabm_standard_variables%latitude,latitude)
   if (longitude/=invalid_longitude) call model%link_horizontal_data(fabm_standard_variables%longitude,longitude)
   call model%link_scalar(fabm_standard_variables%number_of_days_since_start_of_the_year,decimal_yearday)
375

376
   ! Read forcing data specified in input.yaml.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
377
   call init_input_from_file('input.yaml')
378

Jorn Bruggeman's avatar
Jorn Bruggeman committed
379 380 381 382
   ! Request computation of contributions by BGC models to surface albedo and light attenuation
   call model%require_interior_data(fabm_standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux)
   call model%require_horizontal_data(fabm_standard_variables%surface_albedo)

383
   ! Check whether all dependencies of biogeochemical models have now been fulfilled.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
384
   call model%start()
385

Jorn Bruggeman's avatar
Jorn Bruggeman committed
386 387 388 389
   ! Get pointers to contributions by BGC models to surface albedo and light attenuation
   bio_extinction => model%get_interior_data(model%get_interior_variable_id(fabm_standard_variables%attenuation_coefficient_of_photosynthetic_radiative_flux))
   bio_albedo => model%get_horizontal_data(model%get_horizontal_variable_id(fabm_standard_variables%surface_albedo))

390
   ! Update time and all time-dependent inputs.
391
   call update_environment(0_timestepkind)
392

393
   ! Perform custom initialization per biogeochemical model
Jorn Bruggeman's avatar
Jorn Bruggeman committed
394 395 396
   call model%initialize_interior_state()
   call model%initialize_surface_state()
   call model%initialize_bottom_state()
397

398 399 400
   ! Let FABM update the light field (requires state variables to be initialized!)
   call update_light()

401
   ! Allow the model to compute all diagnostics, so output for initial time contains sensible values.
402 403
   allocate(current_rhs(size(cc)))
   call get_rhs(.false.,size(cc),cc,current_rhs)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
404 405 406
   call model%link_all_interior_state_data(cc(1:n_int))
   call model%link_all_bottom_state_data  (cc(n_int+1:n_int+n_bt))
   call model%link_all_surface_state_data (cc(n_int+n_bt+1:n_int+n_bt+n_sf))
407

408 409
   call get_conserved_quantities(totals0)
   int_change_in_totals = 0.0_rk
410

411 412
   LEVEL1 'init_output'
   call init_output(start)
413

414
   call do_output(0_timestepkind)
415 416

   STDERR LINE
417

418 419 420
   end subroutine init_run
!EOC

Jorn Bruggeman's avatar
Jorn Bruggeman committed
421
   subroutine init_input_from_file(path)
422 423
      use yaml_types
      use yaml,yaml_parse=>parse,yaml_error_length=>error_length
424

Jorn Bruggeman's avatar
Jorn Bruggeman committed
425 426
      character(len=*),intent(in) :: path

427
      logical                            :: exists
428
      character(len=yaml_error_length)   :: yaml_error
Jorn Bruggeman's avatar
Jorn Bruggeman committed
429
      class (type_node),         pointer :: root
430

Jorn Bruggeman's avatar
Jorn Bruggeman committed
431 432
      ! Determine whether input configuration file exists. If not, return.
      inquire(file=path,exist=exists)
433 434
      if (.not.exists) return

435
      ! Parse YAML.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
436 437
      root => yaml_parse(path,yaml_unit,yaml_error)
      if (yaml_error/='') call driver%fatal_error('init_input_from_file',trim(yaml_error))
438

439
      ! Process root-level dictionary.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
440
      select type (root)
441
         class is (type_dictionary)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
442
         call init_input_from_yaml_node(root)
443
         class default
Jorn Bruggeman's avatar
Jorn Bruggeman committed
444
         call fatal_error('init_input_from_file',trim(path)//' must contain a dictionary with (variable name : information) pairs,&
Jorn Bruggeman's avatar
Jorn Bruggeman committed
445
               & not a single value.')
446
      end select
Jorn Bruggeman's avatar
Jorn Bruggeman committed
447
   end subroutine init_input_from_file
448

Jorn Bruggeman's avatar
Jorn Bruggeman committed
449
   subroutine init_input_from_yaml_node(mapping)
450
      use yaml_types
451 452 453 454 455

      class (type_dictionary),intent(in)  :: mapping

      character(len=64)                  :: variable_name
      type (type_key_value_pair),pointer :: pair
456 457
      integer                            :: i
      logical                            :: found
458

459
      pair => mapping%first
460
      if (associated(pair)) call driver%log_message('Forcing data specified in input.yaml:')
461 462
      do while (associated(pair))
         variable_name = trim(pair%key)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
463
         if (variable_name=='') call driver%fatal_error('init_input_from_yaml_node','Empty variable name specified.')
464 465 466 467
         found = .false.
         if (model_type==1) then
            select case (variable_name)
            case ('mixed_layer_depth')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
468
               call parse_input_variable(pair%key,pair%value,mixed_layer_depth)
469 470 471 472 473 474
               found = .true.
            case ('mixing_rate')
               call parse_input_variable(pair%key,pair%value,mixing_rate)
               found = .true.
            case default
               do i=1,n_int
Jorn Bruggeman's avatar
Jorn Bruggeman committed
475
                  if (variable_name=='deep/'//trim(model%interior_state_variables(i)%path)) then
476 477 478 479 480 481 482
                     call parse_input_variable(pair%key,pair%value,cc_deep(i))
                     found = .true.
                  end if
               end do
            end select
         end if
         if (.not.found) call parse_input_variable(pair%key,pair%value)
483 484
         pair => pair%next
      end do
Jorn Bruggeman's avatar
Jorn Bruggeman committed
485
   end subroutine init_input_from_yaml_node
486

Jorn Bruggeman's avatar
Jorn Bruggeman committed
487
   subroutine parse_input_variable(variable_name,value_node,input_)
488
      use yaml_types
489

490 491
      character(len=*),        intent(in) :: variable_name
      class (type_node),target,intent(in) :: value_node
Jorn Bruggeman's avatar
Jorn Bruggeman committed
492
      type (type_scalar_input), target, optional :: input_
493

494
      class (type_dictionary),   pointer :: mapping
495 496
      type (type_error),         pointer :: config_error
      class (type_node),         pointer :: node
497
      class (type_scalar),       pointer :: constant_value_node, file_node
498 499
      real(rk)                           :: relaxation_time
      logical                            :: is_state_variable
500
      type (type_key_value_pair),pointer :: pair
Jorn Bruggeman's avatar
Jorn Bruggeman committed
501 502 503
      type (type_fabm_interior_variable_id)   :: interior_id
      type (type_fabm_horizontal_variable_id) :: horizontal_id
      type (type_fabm_scalar_variable_id)     :: scalar_id
Jorn Bruggeman's avatar
Jorn Bruggeman committed
504
      type (type_scalar_input), pointer  :: input
505 506 507 508 509 510 511

      select type (value_node)
      class is (type_dictionary)
         mapping => value_node
      class default
         call fatal_error('init_input_from_yaml_node','Contents of '//trim(value_node%path)//' must be a dictionary, not a single value.')
      end select
512 513

      is_state_variable = .false.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
514 515 516 517 518 519
      if (present(input_)) then
         input => input_
      else
         allocate(input)
         call extra_inputs%add(input)

520
         ! Try to locate the forced variable among interior, horizontal, and global variables in the active biogeochemical models.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
521 522
         interior_id = model%get_interior_variable_id(variable_name)
         if (model%is_variable_used(interior_id)) then
523
            is_state_variable = interior_id%variable%state_indices%value/=-1
524
         else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
525
            horizontal_id = model%get_horizontal_variable_id(variable_name)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
526
            if (model%is_variable_used(horizontal_id)) then
527 528
               is_state_variable = horizontal_id%variable%state_indices%value/=-1
            else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
529
               scalar_id = model%get_scalar_variable_id(variable_name)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
530
               if (.not. model%is_variable_used(scalar_id)) call log_message('WARNING! input.yaml: &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
531
                  &Variable "'//trim(variable_name)//'" is not present in any biogeochemical model.')
532
            end if
533 534 535
         end if
      end if

536
      ! Prepend to list of input data.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
537
      input%name = trim(variable_name)
538

539 540 541
      constant_value_node => mapping%get_scalar('constant_value',required=.false.,error=config_error)
      file_node => mapping%get_scalar('file',required=.false.,error=config_error)
      if (associated(constant_value_node)) then
542
         ! Input variable is set to a constant value.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
543
         input%method = 0
Jorn Bruggeman's avatar
Jorn Bruggeman committed
544
         input%constant_value = mapping%get_real('constant_value',error=config_error)
545 546
         if (associated(config_error)) call fatal_error('parse_input_variable',config_error%message)

547
         ! Make sure keys related to time-varying input are not present.
548
         if (associated(file_node)) call fatal_error('parse_input_variable','input.yaml, variable "'//trim(variable_name)//'": &
Jorn Bruggeman's avatar
Jorn Bruggeman committed
549
            &keys "constant_value" and "file" cannot both be present.')
550
         node => mapping%get('column')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
551 552
         if (associated(node)) call fatal_error('parse_input_variable','input.yaml, variable "'//trim(variable_name)//'": &
            &keys "constant_value" and "column" cannot both be present.')
553
         node => mapping%get('scale_factor')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
554 555
         if (associated(node)) call fatal_error('parse_input_variable','input.yaml, variable "'//trim(variable_name)//'": &
            &keys "constant_value" and "scale_factor" cannot both be present.')
556
      elseif (associated(file_node)) then
557
         ! Input variable is set to a time-varying value. Obtain path, column number and scale factor.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
558 559
         input%method = 2
         input%path = mapping%get_string('file',error=config_error)
560
         if (associated(config_error)) call fatal_error('parse_input_variable',config_error%message)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
561
         input%index = mapping%get_integer('column',default=1,error=config_error)
562
         if (associated(config_error)) call fatal_error('parse_input_variable',config_error%message)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
563
         input%scale_factor = mapping%get_real('scale_factor',default=1.0_rk,error=config_error)
564
         if (associated(config_error)) call fatal_error('parse_input_variable',config_error%message)
565 566 567
      else
         call fatal_error('parse_input_variable','input.yaml, variable "'//trim(variable_name)//'": &
            &either key "constant_value" or key "file" must be present.')
568
      end if
Jorn Bruggeman's avatar
Jorn Bruggeman committed
569
      call register_input(input)
570 571

      if (is_state_variable) then
572
         ! This is a state variable. Obtain associated relaxation time.
573 574 575
         relaxation_time = mapping%get_real('relaxation_time',default=1.e15_rk,error=config_error)
         if (associated(config_error)) call fatal_error('parse_input_variable',config_error%message)
      else
576
         ! This is not a state variable. Make sure no relaxation time is specified.
577
         node => mapping%get('relaxation_time')
Jorn Bruggeman's avatar
Jorn Bruggeman committed
578 579
         if (associated(node)) call fatal_error('parse_input_variable','input.yaml, variable "'//trim(variable_name)//'": &
            &key "relaxation_time" is not supported because "'//trim(variable_name)//'" is not a state variable.')
580
      end if
581

582 583 584
      ! Warn about uninterpreted keys.
      pair => mapping%first
      do while (associated(pair))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
585 586
         if (.not.pair%accessed) call fatal_error('parse_input_variable','input.yaml: &
            &Unrecognized option "'//trim(pair%key)//'" found for variable "'//trim(variable_name)//'".')
587 588 589
         pair => pair%next
      end do

590
      ! If a data pointer was provided, this variable for the host, not FABM, so return.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
591
      if (present(input_)) return
592

593
      ! Link forced data to target variable.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
594
      if (model%is_variable_used(interior_id)) then
Jorn Bruggeman's avatar
Jorn Bruggeman committed
595
         call model%link_interior_data(interior_id, input%value, source=data_source_user)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
596
      elseif (model%is_variable_used(horizontal_id)) then
Jorn Bruggeman's avatar
Jorn Bruggeman committed
597
         call model%link_horizontal_data(horizontal_id, input%value, source=data_source_user)
598
      else
Jorn Bruggeman's avatar
Jorn Bruggeman committed
599
         call model%link_scalar(scalar_id, input%value, source=data_source_user)
600
      end if
601

602
   end subroutine parse_input_variable
603

604
   subroutine update_environment(n)
605
      integer(timestepkind),intent(in) :: n
606

607
      ! Update time in time manager
608 609
      call update_time(n)

610
      ! Compute decimal year day (input for some biogeochemical models)
611
      decimal_yearday = yearday-1 + dble(secondsofday)/86400
612

613
      ! Update environment (i.e., read from input files)
614 615
      call do_input(julianday,secondsofday)

Jorn Bruggeman's avatar
Jorn Bruggeman committed
616 617
      if (model_type==1) column_depth = mixed_layer_depth%value

618
      ! Compute density from temperature and salinity, if required by biogeochemistry.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
619
      if (compute_density) dens = rho_feistel(salt%value,temp%value,5._rk*column_depth,.true.)
620 621 622
   end subroutine update_environment

   subroutine update_light()
Jorn Bruggeman's avatar
Jorn Bruggeman committed
623 624 625 626
      real(rk) :: zenith_angle,solar_zenith_angle
      real(rk) :: shortwave_radiation
      real(rk) :: albedo,albedo_water
      real(rk) :: hh
627

628
      ! Calculate photosynthetically active radiation at surface, if it is not provided in the input file.
629 630
      if (swr_method==0) then
         ! Calculate photosynthetically active radiation from geographic location, time, cloud cover.
631
         hh = secondsofday*(1._rk/3600)
632
         zenith_angle = solar_zenith_angle(yearday,hh,longitude,latitude)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
633
         swr_sf = shortwave_radiation(zenith_angle,yearday,longitude,latitude,cloud)
634 635
         if (albedo_correction) then
            albedo = albedo_water(1,zenith_angle,yearday)
636
            swr_sf = swr_sf*(1._rk-albedo-bio_albedo)
637
         end if
Jorn Bruggeman's avatar
Jorn Bruggeman committed
638 639
      else
         swr_sf = light%value
640 641 642
      end if

      ! Multiply by fraction of short-wave radiation that is photosynthetically active.
643
      par_sf = par_fraction*swr_sf
644 645 646

      ! Apply light attentuation with depth, unless local light is provided in the input file.
      if (swr_method/=2) then
647 648
         ! Calculate light extinction
         extinction = 0.0_rk
Jorn Bruggeman's avatar
Jorn Bruggeman committed
649
         if (apply_self_shading) extinction = bio_extinction
650 651
         extinction = extinction + par_background_extinction

652 653 654 655 656 657 658 659 660 661
         ! Either we calculate surface PAR, or surface PAR is provided.
         ! Calculate the local PAR at the given depth from par fraction, extinction coefficient, and depth.
         par_ct = par_sf*exp(-0.5_rk*column_depth*extinction)
         par_bt = par_sf*exp(-column_depth*extinction)
      else
         par_ct = par_sf
         par_bt = par_sf
      end if
      call update_depth(CENTER)

662
   end subroutine update_light
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: Manage global time--stepping \label{timeLoop}
!
! !INTERFACE:
   subroutine time_loop()
!
! !DESCRIPTION:
! This internal routine is the heart of the code. It contains
! the main time-loop inside of which all routines required
! during the time step are called.
!
! !REVISION HISTORY:
!  Original author(s): Jorn Bruggeman
!
! !LOCAL VARIABLES:
681
   logical                   :: valid_state
682
   integer                   :: progress,k
683 684 685 686 687
!EOP
!-----------------------------------------------------------------------
!BOC
   LEVEL1 'time_loop'

688 689
   progress = (MaxN-MinN+1)/10
   k = 0
Jorn Bruggeman's avatar
Jorn Bruggeman committed
690 691
   do itime=MinN,MaxN
      if(mod(itime,progress) == 0 .or. itime == MinN) then
692 693
         LEVEL0 k,'%'
         k = k+10
694 695
      end if

696
      ! Update time and all time-dependent inputs.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
697
      call update_environment(itime)
698
      call update_light()
699

700
      ! Integrate one time step
Jorn Bruggeman's avatar
Jorn Bruggeman committed
701
      call ode_solver(ode_method,size(cc),timestep,cc,get_rhs,get_ppdd)
702
      call get_rhs(.false.,size(cc),cc,current_rhs)
703

Jorn Bruggeman's avatar
Jorn Bruggeman committed
704
      ! ODE solver may have redirected the current state to an array with intermediate values.
705
      ! Reset to global array.
Jorn Bruggeman's avatar
Jorn Bruggeman committed
706 707 708
      call model%link_all_interior_state_data(cc(1:n_int))
      call model%link_all_bottom_state_data  (cc(n_int+1:n_int+n_bt))
      call model%link_all_surface_state_data (cc(n_int+n_bt+1:n_int+n_sf+n_bt))
709

710
      ! Verify whether the model state is still valid (clip if needed and allowed)
Jorn Bruggeman's avatar
Jorn Bruggeman committed
711 712 713
      call model%check_interior_state(repair_state,valid_state)
      if (valid_state .or. repair_state) call model%check_bottom_state(repair_state,valid_state)
      if (valid_state .or. repair_state) call model%check_surface_state(repair_state,valid_state)
714 715 716 717
      if (.not. (valid_state .or. repair_state)) &
         call fatal_error('time_loop','State variable values are invalid and repair is not allowed. &
            &This may be fixed by setting repair_state=.true. (clip state to nearest valid value), &
            &but this should be used with caution. Try and decrease the time step (dt) first - and see if that helps.')
718

719 720 721
      if (compute_conserved_quantities) then
         call get_conserved_quantities(int_change_in_totals)
         int_change_in_totals = int_change_in_totals - totals0
722
      end if
723

724
      ! Do output
Jorn Bruggeman's avatar
Jorn Bruggeman committed
725
      call do_output(itime)
726 727 728 729 730 731
   end do
   STDERR LINE

   end subroutine time_loop
!EOC

732 733 734
   subroutine get_conserved_quantities(depth_int_totals)
      real(rk), intent(inout) :: depth_int_totals(size(model%conserved_quantities))
      real(rk) :: totals_hz(size(model%conserved_quantities))
Jorn Bruggeman's avatar
Jorn Bruggeman committed
735 736
      call model%get_interior_conserved_quantities(totals)
      call model%get_horizontal_conserved_quantities(totals_hz)
737 738
      depth_int_totals = totals*column_depth + totals_hz
   end subroutine