-
Notifications
You must be signed in to change notification settings - Fork 93
/
Copy pathPRTGenericMod.F90
1510 lines (1088 loc) · 61.7 KB
/
PRTGenericMod.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
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
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module PRTGenericMod
! ------------------------------------------------------------------------------------
! Plant Allocation and Reactive Transport (PART) +
! Extensible Hypotheses (EH) = PARTEH
!
! Non-Specific (Generic) Classes and Functions
! This contains the base classes for both the variables and the global class
!
! General idea: PARTEH treats its state variables as objects. Each object
! can be mapped to, or associated with:
! 1) an organ
! 2) a spatial position associated with that organ
! 3) a chemical element (ie carbon isotope or nutrient), aka chemical species
!
!
! THIS ROUTINE SHOULD NOT HAVE TO BE MODIFIED TO ACCOMODATE NEW HYPOTHESES
! (in principle ...)
!
! Ryan Knox, April 2018
! ------------------------------------------------------------------------------------
use FatesConstantsMod, only : r8 => fates_r8
use FatesConstantsMod, only : i4 => fates_int
use FatesConstantsMod, only : nearzero
use FatesConstantsMod, only : calloc_abs_error
use FatesConstantsMod, only : years_per_day
use FatesConstantsMod, only : days_per_sec
use FatesGlobals , only : endrun => fates_endrun
use FatesGlobals , only : fates_log
use shr_log_mod , only : errMsg => shr_log_errMsg
use PRTParametersMod , only : prt_params
implicit none
private ! Modules are private by default
character(len=*), parameter, private :: sourcefile = &
__FILE__
integer, parameter, public :: maxlen_varname = 128
integer, parameter, public :: maxlen_varsymbol = 32
integer, parameter, public :: maxlen_varunits = 32
integer, parameter, public :: len_baseunit = 6
! We use this parameter as the value for which we set un-initialized values
real(r8), parameter, public :: un_initialized = -9.9e32_r8
! We use this parameter as the value for which we check un-initialized values
real(r8), parameter, public :: check_initialized = -8.8e32_r8
! -------------------------------------------------------------------------------------
! IMPORTANT!
! All elements in all organs should be expressed in terms of KILOGRAMS
! All rates of change are expressed in terms of kilograms / day
! This assumption cannot be broken!
! -------------------------------------------------------------------------------------
character(len=len_baseunit), parameter, public :: mass_unit = 'kg'
character(len=len_baseunit), parameter, public :: mass_rate_unit = 'kg/day'
! -------------------------------------------------------------------------------------
! Allocation Hypothesis Types
! These should each have their own module
! -------------------------------------------------------------------------------------
integer, parameter, public :: prt_carbon_allom_hyp = 1
integer, parameter, public :: prt_cnp_flex_allom_hyp = 2
! -------------------------------------------------------------------------------------
! Organ types
! These are public indices used to map the organs
! in each hypothesis to organs that acknowledged in the calling model
! -------------------------------------------------------------------------------------
integer, parameter, public :: num_organ_types = 6
integer, parameter, public :: all_organs = 0 ! index for all organs
integer, parameter, public :: leaf_organ = 1 ! index for leaf organs
integer, parameter, public :: fnrt_organ = 2 ! index for fine-root organs
integer, parameter, public :: sapw_organ = 3 ! index for sapwood organs
integer, parameter, public :: store_organ = 4 ! index for storage organs
integer, parameter, public :: repro_organ = 5 ! index for reproductive organs
integer, parameter, public :: struct_organ = 6 ! index for structure (dead) organs
! -------------------------------------------------------------------------------------
! Element types
! These are public indices used to map the elements (chem species) in each hypothesis
! to the element that are acknowledged in the calling model
! -------------------------------------------------------------------------------------
integer, parameter, public :: num_element_types = 6 ! Total number of unique element
! curently recognized by PARTEH
! should be max index in list below
! The following list are the unique indices associated with the
! element used in each hypothesis. Note these are just POTENTIAL
! element. At the time of writing this, we are very far away from
! creating allocation schemes that even use potassium.
integer, parameter, public :: carbon12_element = 1
integer, parameter, public :: carbon13_element = 2
integer, parameter, public :: carbon14_element = 3
integer, parameter, public :: nitrogen_element = 4
integer, parameter, public :: phosphorus_element = 5
integer, parameter, public :: potassium_element = 6
! The following elements are just placeholders. In the future
! if someone wants to develope an allocation hypothesis
! that uses nickel, we can just uncomment it from this list
! integer, parameter, public :: calcium_element = 7
! integer, parameter, public :: magnesium_element = 8
! integer, parameter, public :: sulfur_element = 9
! integer, parameter, public :: chlorine_element = 10
! integer, parameter, public :: iron_element = 11
! integer, parameter, public :: manganese_element = 12
! integer, parameter, public :: zinc_element = 13
! integer, parameter, public :: copper_element = 14
! integer, parameter, public :: boron_element = 15
! integer, parameter, public :: molybdenum_element = 16
! integer, parameter, public :: nickel_element = 17
! We have some lists of elements or lists of organs, such as
! a list of all carbon elements. To keep routines simple
! we set a global to the maximum list size for scratch arrays.
integer, parameter, public :: max_spec_per_group = 3 ! we may query these lists
! the carbon elements are the biggest list
! right now
! List of all carbon elements, the special index "all_carbon_elements"
! implies the following list of carbon organs (NOT USED)
integer, parameter, dimension(3), public :: carbon_elements_list = &
[carbon12_element, carbon13_element, carbon14_element]
! This is the maximum number of leaf age pools allowed on each plant
! (used for allocating scratch space)
integer, parameter, public :: max_nleafage = 4
! This is the minimum allowable L2FR, this is needed so that plants
! in the understory don't shrink their roots down so far that
! they dissappear and cause numerical issues
real(r8), parameter, public :: l2fr_min = 0.01_r8
! -------------------------------------------------------------------------------------
!
! The following is the data structure that holds the state (ie carbon,
! nutrients, etc) for each pool of each plant.
!
! For example, this could be the carbon 12 of the leaf pool; its instantaneous state,
! and its fluxes.
!
! Note also that these are vectors and not scalars, which indicates that there
! may be more than 1 discrete spatial positions. For instance, there might be multiple
! leaf layers or something.
!
! Since there are many variables, as well as boundary conditions, this object is
! NESTED in the prt_vartypes (<---- see the "s" at the end?) structure that follows.
!
! Each object will have a unique index associated with it, it will also be mapped
! to a specific organ and element combination.
!
! It is assumed that over the control period (probably 1 day) that
! changes in the current state (val) relative to the value at the start of the
! control period (val0), are equal to the time integrated flux terms
! (net_alloc, turnover, etc)
!
! -------------------------------------------------------------------------------------
type, public :: prt_vartype
real(r8),pointer :: val(:) ! Instantaneous state variable [kg]
real(r8),allocatable :: val0(:) ! State variable at the beginning
! of the control period [kg]
real(r8),allocatable :: net_alloc(:) ! Net change due to allocation/transport [kg]
! over the control period [kg]
real(r8),allocatable :: turnover(:) ! Losses due to turnover [kg]
! or, any mass destined for litter
! over the control period
real(r8),allocatable :: burned(:) ! Losses due to burn [kg]
real(r8),allocatable :: damaged(:) ! Losses due to damage [kg]
real(r8),allocatable :: herbivory(:) ! Losses due to herbivory [kg]
! Placeholder
! To save on memory, keep this commented out, or simply
! add this only in the extension ... ?
! real(r8),allocatable :: coordinate(:,:)
end type prt_vartype
! -------------------------------------------------------------------------------------
! Input boundary conditions. These will be allocated as an array for each plant.
! This type will also be broken into 3 types of boundary conditions: input only,
! output only, and input-output.
! -------------------------------------------------------------------------------------
type, public :: prt_bctype
real(r8), pointer :: rval
integer, pointer :: ival
end type prt_bctype
! -------------------------------------------------------------------------------------
! The following is the object that is directly attached to each plant.
!
! ie this is the parent object.
! It contains the state variable object: variables
! as well as the boundary condition pointers bc_inout, bc_in and bc_out
!
! This object also contains the bulk of the PRT routines, including
! extended (hypothesis specific routines) and generic routines (eg
! routines that can operate on any hypothesis)
!
! There are procedures that are specialized for each module. And then
! there are procedures that are supposed to be generic and should support
! all the different modules.
! -------------------------------------------------------------------------------------
type, public :: prt_vartypes
type(prt_vartype),allocatable :: variables(:) ! The state variables and fluxes
type(prt_bctype), allocatable :: bc_inout(:) ! These boundaries may be changed
type(prt_bctype), allocatable :: bc_in(:) ! These are protected
type(prt_bctype), allocatable :: bc_out(:) ! These are overwritten
real(r8) :: ode_opt_step
contains
! These are extendable procedures that have specialized
! content in each of the different hypotheses
procedure :: DailyPRT => DailyPRTBase
procedure :: FastPRT => FastPRTBase
procedure :: DamageRecovery => DamageRecoveryBase
procedure :: GetNutrientTarget => GetNutrientTargetBase
! These are generic functions that should work on all hypotheses
procedure, non_overridable :: InitAllocate
procedure, non_overridable :: InitPRTVartype
procedure, non_overridable :: FlushBCs
procedure, non_overridable :: InitializeInitialConditions
procedure, non_overridable :: CheckInitialConditions
procedure, non_overridable :: RegisterBCIn
procedure, non_overridable :: RegisterBCOut
procedure, non_overridable :: RegisterBCInout
procedure, non_overridable :: GetState
procedure, non_overridable :: GetTurnover
procedure, non_overridable :: GetBurned
procedure, non_overridable :: GetHerbivory
procedure, non_overridable :: GetNetAlloc
procedure, non_overridable :: ZeroRates
procedure, non_overridable :: CheckMassConservation
procedure, non_overridable :: DeallocatePRTVartypes
procedure, non_overridable :: WeightedFusePRTVartypes
procedure, non_overridable :: CopyPRTVartypes
procedure :: AgeLeaves ! This routine may be used generically
! but also leaving the door open for over-rides
end type prt_vartypes
! Global identifiers for which elements we are using (apply mostly to litter)
integer, public :: num_elements ! This is the number of elements in this simulation
! e.g. (C,N,P,K, etc)
integer, allocatable, public :: element_list(:) ! This vector holds the list of global element identifiers
! examples are carbon12_element
! nitrogen_element, etc.
integer, public :: element_pos(num_element_types) ! This is the reverse lookup
! for element types. Pick an element
! global index, and it gives you
! the position in the element_list
! -------------------------------------------------------------------------------------
! This next section contains the objects that describe the mapping for each specific
! hypothesis. It is also a way to call the descriptions of variables for any
! arbitrary hypothesis.
! These are things that are globally true, not specific to each plant.
! For instance the map just contains the list of variable names, not the values for
! each plant.
! These are not instanced on every plant, they are just instanced once on every model
! machine or memory space. They should only be initialized once and used
! as read-only from that point on.
! -------------------------------------------------------------------------------------
! -------------------------------------------------------------------------------------
! This type simply packs the names and symbols associated with all
! the variables for any given hypothesis
! -------------------------------------------------------------------------------------
type, public :: state_descriptor_type
character(len=maxlen_varname) :: longname
character(len=maxlen_varsymbol) :: symbol
integer :: organ_id ! global id for organ
integer :: element_id ! global id for element
integer :: num_pos ! number of descrete spatial positions
! Also, will probably need flags to define different types of groups that this variable
! belongs too, which will control things like fusion, normalization, when to zero, etc...
end type state_descriptor_type
! This type will help us loop through all the different variables associated
! with a specific organ type. Since variables are a combination of organ and
! element, the number of unique variables is capped at the number of elements
! per each organ.
type, public :: organ_map_type
integer, dimension(1:num_element_types) :: var_id
integer :: num_vars
end type organ_map_type
! This structure packs both the mapping structure and the variable descriptors
! -------------------------------------------------------------------------------------
! This array should contain the lists of indices to
! the element x organ variable structure that is used to map variables to the outside
! world.
!
!
! | carbon | nitrogen | phosphorus | .... |
! ------------------------------------------
! leaf | | | | |
! fine-root | | | | |
! sapwood | | | | |
! storage | | | | |
! reproduction | | | | |
! structure | | | | |
! .... | | | | |
! ------------------------------------------
!
! -------------------------------------------------------------------------------------
type, public :: prt_global_type
! Note that index 0 is reserved for "all" or "irrelevant"
character(len=maxlen_varname) :: hyp_name
! This is the hypothesis index, used internally for some non-specific
! routines where different methods are applied
integer :: hyp_id
! This will save the specific variable id associated with each organ and element
integer, dimension(0:num_organ_types,0:num_element_types) :: sp_organ_map
! This holds the verbose descriptions of the variables, symbols, names, etc
type(state_descriptor_type), allocatable :: state_descriptor(:)
! This will save the list of variable ids associated with each organ. There
! are multiple of these because we may have multiple element per organ.
type(organ_map_type), dimension(1:num_organ_types) :: organ_map
! The number of input boundary conditions
integer :: num_bc_in
! The number of output boundary conditions
integer :: num_bc_out
! The number of combo input-output boundary conditions
integer :: num_bc_inout
! The number of variables set by each hypothesis
integer :: num_vars
contains
procedure, non_overridable :: ZeroGlobal
procedure, non_overridable :: RegisterVarInGlobal
end type prt_global_type
class(prt_global_type),pointer,public :: prt_global
! Make necessary procedures public
public :: GetCoordVal
public :: SetState
public :: StorageNutrientTarget
contains
! =====================================================================================
! Module Functions and Subroutines
! =====================================================================================
subroutine ZeroGlobal(this)
! This subroutine zero's out the map between variable indexes and the
! elements and organs they are associated with.
! It also sets the counts of the variables and boundary conditions as
! a nonsense number that will trigger a fail if they are specified later.
! This routine must be called
class(prt_global_type) :: this
integer :: io ! Organ loop counter
integer :: is ! Element loop counter
! First zero out the array
do io = 1,num_organ_types
do is = 1,num_element_types
this%sp_organ_map(io,is) = 0
this%organ_map(io)%var_id(is) = 0
end do
this%organ_map(io)%num_vars = 0
end do
! Set the number of boundary conditions as a bogus value
this%num_bc_in = -9
this%num_bc_out = -9
this%num_bc_inout = -9
! Set the number of variables to a bogus value. This should be
! immediately over-written in the routine that is calling this
this%num_vars = -9
return
end subroutine ZeroGlobal
! =====================================================================================
subroutine RegisterVarInGlobal(this, var_id, long_name, symbol, organ_id, element_id, num_pos)
! This subroutine is called for each variable that is defined in each specific hypothesis.
! For instance, this is called six times in the carbon only hypothesis,
! each time providing names, symbols, associated organs and element for each pool.
class(prt_global_type) :: this
integer, intent(in) :: var_id
character(len=*),intent(in) :: long_name
character(len=*),intent(in) :: symbol
integer, intent(in) :: organ_id
integer, intent(in) :: element_id
integer, intent(in) :: num_pos
! Set the descriptions and the associated organs/element in the variable's
! own array
this%state_descriptor(var_id)%longname = long_name
this%state_descriptor(var_id)%symbol = symbol
this%state_descriptor(var_id)%organ_id = organ_id
this%state_descriptor(var_id)%element_id = element_id
this%state_descriptor(var_id)%num_pos = num_pos
! Set the mapping tables for the external model
this%sp_organ_map(organ_id,element_id) = var_id
! Set another map that helps to locate all the relevant pools associated
! with an organ
this%organ_map(organ_id)%num_vars = this%organ_map(organ_id)%num_vars + 1
this%organ_map(organ_id)%var_id(this%organ_map(organ_id)%num_vars) = var_id
return
end subroutine RegisterVarInGlobal
! =====================================================================================
subroutine InitPRTVartype(this)
class(prt_vartypes) :: this
! This subroutine should be the first call whenever a prt_vartype object is
! instantiated.
!
! Most likely, this will occur whenever a new plant or cohort is created.
!
! This routine handles the allocation (extended procedure)
! and then the initializing of states with bogus information, and then
! the flushing of all boundary conditions to null.
call this%InitAllocate() ! Allocate memory spaces
call this%InitializeInitialConditions() ! Set states to a nan-like starter value
call this%FlushBCs() ! Set all boundary condition pointers
! to null
return
end subroutine InitPRTVartype
! =====================================================================================
subroutine InitAllocate(this)
! ----------------------------------------------------------------------------------
! This initialization is called everytime a plant/cohort
! is newly recruited. Like the name implies, we are just allocating space here.
! ----------------------------------------------------------------------------------
class(prt_vartypes) :: this
integer :: i_var ! Variable loop index
integer :: num_pos ! The number of positions for each variable
! Allocate the boundar condition arrays and flush them to no-data flags
! ----------------------------------------------------------------------------------
if(prt_global%num_bc_in > 0) then
allocate(this%bc_in(prt_global%num_bc_in))
end if
if(prt_global%num_bc_inout > 0) then
allocate(this%bc_inout(prt_global%num_bc_inout))
end if
if(prt_global%num_bc_out > 0) then
allocate(this%bc_out(prt_global%num_bc_out))
end if
! Allocate the state variables
allocate(this%variables(prt_global%num_vars))
do i_var = 1, prt_global%num_vars
num_pos = prt_global%state_descriptor(i_var)%num_pos
allocate(this%variables(i_var)%val(num_pos))
allocate(this%variables(i_var)%val0(num_pos))
allocate(this%variables(i_var)%turnover(num_pos))
allocate(this%variables(i_var)%net_alloc(num_pos))
allocate(this%variables(i_var)%burned(num_pos))
allocate(this%variables(i_var)%damaged(num_pos))
allocate(this%variables(i_var)%herbivory(num_pos))
end do
return
end subroutine InitAllocate
! =====================================================================================
subroutine InitializeInitialConditions(this)
! ----------------------------------------------------------------------------------
! This routine sets all PARTEH variables to a nonsense value.
! This ensures that a fail is triggered if a value is not initialized correctly.
! ----------------------------------------------------------------------------------
class(prt_vartypes) :: this
integer :: i_var ! Variable index
do i_var = 1, prt_global%num_vars
this%variables(i_var)%val(:) = un_initialized
this%variables(i_var)%val0(:) = un_initialized
this%variables(i_var)%turnover(:) = un_initialized
this%variables(i_var)%burned(:) = un_initialized
this%variables(i_var)%damaged(:) = un_initialized
this%variables(i_var)%net_alloc(:) = un_initialized
this%variables(i_var)%herbivory(:) = un_initialized
end do
! Initialize the optimum step size as very large.
this%ode_opt_step = 1e6_r8
return
end subroutine InitializeInitialConditions
! =============================================================
subroutine CheckInitialConditions(this)
! This subroutine makes sure that every variable defined
! in the hypothesis has been given an initial value.
!
! This should be called following any blocks where initial
! conditions are set. In fates, these calls already
! exist and when new hypotheses are added, they will
! already be checked if the initial conditions are
! specified in parallel with the other hypotheses.
class(prt_vartypes) :: this
integer :: i_var ! index for iterating variables
integer :: n_cor_ids ! Number of coordinate ids
integer :: i_cor ! index for iterating coordinate dimension
integer :: i_organ ! The global organ id for this variable
integer :: i_element ! The global element id for this variable
do i_var = 1, prt_global%num_vars
n_cor_ids = size(this%variables(i_var)%val,1)
do i_cor = 1, n_cor_ids
if(this%variables(i_var)%val(i_cor) < check_initialized) then
i_organ = prt_global%state_descriptor(i_var)%organ_id
i_element = prt_global%state_descriptor(i_var)%element_id
write(fates_log(),*)'Not all initial conditions for state variables'
write(fates_log(),*)' in PRT hypothesis: ',trim(prt_global%hyp_name)
write(fates_log(),*)' were written out.'
write(fates_log(),*)' i_var: ',i_var
write(fates_log(),*)' i_cor: ',i_cor
write(fates_log(),*)' organ_id:',i_organ
write(fates_log(),*)' element_id',i_element
write(fates_log(),*)'Exiting'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
end do
end do
return
end subroutine CheckInitialConditions
! =====================================================================================
subroutine FlushBCs(this)
! Boundary conditions are pointers to real's and integers in the calling model.
! To flush these, we set all pointers to null
! Arguments
class(prt_vartypes) :: this
! Local
integer :: num_bc_in
integer :: num_bc_out
integer :: num_bc_inout
integer :: i
if(allocated(this%bc_in))then
num_bc_in = size(this%bc_in,1)
do i = 1,num_bc_in
this%bc_in(i)%rval => null()
this%bc_in(i)%ival => null()
end do
end if
if(allocated(this%bc_out))then
num_bc_out = size(this%bc_out,1)
do i = 1,num_bc_out
this%bc_out(i)%rval => null()
this%bc_out(i)%ival => null()
end do
end if
if(allocated(this%bc_inout))then
num_bc_inout = size(this%bc_inout,1)
do i = 1,num_bc_inout
this%bc_inout(i)%rval => null()
this%bc_inout(i)%ival => null()
end do
end if
return
end subroutine FlushBCs
! =====================================================================================
subroutine RegisterBCIn(this,bc_id, bc_rval, bc_ival )
! This routine must be called once for each "input only" boundary condition of each
! hypothesis.
! The group of calls only needs to happen once, following InitPRTVartype.
! Since we use pointers, we don't need to constantly ask for new boundary conditions
!
! The only complication to this would occur, if the boundary condition variable
! that these pointers point to is being disassociated. In that case, one would
! need to re-register that boundary condition variable.
! Input Arguments
class(prt_vartypes) :: this
integer,intent(in) :: bc_id
real(r8),optional, intent(inout), target :: bc_rval
integer, optional, intent(inout), target :: bc_ival
if(present(bc_ival)) then
this%bc_in(bc_id)%ival => bc_ival
end if
if(present(bc_rval)) then
this%bc_in(bc_id)%rval => bc_rval
end if
return
end subroutine RegisterBCIn
! =====================================================================================
subroutine RegisterBCOut(this,bc_id, bc_rval, bc_ival )
! This routine is similar to the routine above RegisterBCIn, except this
! is for registering "output only" boundary conditions.
! Input Arguments
class(prt_vartypes) :: this
integer,intent(in) :: bc_id
real(r8), optional, intent(inout),target :: bc_rval
integer, optional, intent(inout),target :: bc_ival
if(present(bc_ival)) then
this%bc_out(bc_id)%ival => bc_ival
end if
if(present(bc_rval)) then
this%bc_out(bc_id)%rval => bc_rval
end if
return
end subroutine RegisterBCOut
! =====================================================================================
subroutine RegisterBCInOut(this,bc_id, bc_rval, bc_ival )
! This routine is similar to the two routines above, except this
! is for registering "input-output" boundary conditions.
! These are conditions that are passed into PARTEH, and are expected
! to be updated (or not), and passed back to the host (FATES).
! Input Arguments
class(prt_vartypes) :: this
integer,intent(in) :: bc_id
real(r8), optional, intent(inout),target :: bc_rval
integer, optional, intent(inout),target :: bc_ival
if(present(bc_ival)) then
this%bc_inout(bc_id)%ival => bc_ival
end if
if(present(bc_rval)) then
this%bc_inout(bc_id)%rval => bc_rval
end if
return
end subroutine RegisterBCInOut
! =====================================================================================
subroutine CopyPRTVartypes(this, donor_prt_obj)
! Here we copy over all information from a donor_prt_object into the current PRT
! object. It is assumed that the current PRT object
! has already been initialized ( ie. InitAllocate() )
! variable val0 is omitted, because it is ephemeral and used only during the
! allocation process
! Arguments
class(prt_vartypes) :: this
class(prt_vartypes), intent(in), pointer :: donor_prt_obj
! Locals
integer :: i_var ! loop iterator for variable objects
integer :: i_bc ! loop iterator for boundary pointers
integer :: num_bc_in
integer :: num_bc_inout
integer :: num_bc_out
do i_var = 1, prt_global%num_vars
this%variables(i_var)%val(:) = donor_prt_obj%variables(i_var)%val(:)
this%variables(i_var)%val0(:) = donor_prt_obj%variables(i_var)%val0(:)
this%variables(i_var)%net_alloc(:) = donor_prt_obj%variables(i_var)%net_alloc(:)
this%variables(i_var)%turnover(:) = donor_prt_obj%variables(i_var)%turnover(:)
this%variables(i_var)%burned(:) = donor_prt_obj%variables(i_var)%burned(:)
this%variables(i_var)%damaged(:) = donor_prt_obj%variables(i_var)%damaged(:)
this%variables(i_var)%herbivory(:) = donor_prt_obj%variables(i_var)%herbivory(:)
end do
this%ode_opt_step = donor_prt_obj%ode_opt_step
return
end subroutine CopyPRTVartypes
! =====================================================================================
subroutine WeightedFusePRTVartypes(this,donor_prt_obj, recipient_fuse_weight)
! This subroutine fuses two PRT objects together based on a fusion weighting
! assigned for the recipient (the class calling this)
! Arguments
class(prt_vartypes) :: this
class(prt_vartypes), intent(in), pointer :: donor_prt_obj
real(r8),intent(in) :: recipient_fuse_weight ! This is the weighting
! for the recipient
! integer,intent(in),optional :: position_id
! Locals
integer :: i_var ! Loop iterator over variables
integer :: pos_id ! coordinate id (defaults to 1, if not position_id)
! If no argument regarding positions is supplied, then fuse all positions
! sequentially
do i_var = 1, prt_global%num_vars
do pos_id = 1,prt_global%state_descriptor(i_var)%num_pos
this%variables(i_var)%val(pos_id) = recipient_fuse_weight * this%variables(i_var)%val(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val(pos_id)
this%variables(i_var)%val0(pos_id) = recipient_fuse_weight * this%variables(i_var)%val0(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%val0(pos_id)
this%variables(i_var)%net_alloc(pos_id) = recipient_fuse_weight * this%variables(i_var)%net_alloc(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%net_alloc(pos_id)
this%variables(i_var)%turnover(pos_id) = recipient_fuse_weight * this%variables(i_var)%turnover(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%turnover(pos_id)
this%variables(i_var)%burned(pos_id) = recipient_fuse_weight * this%variables(i_var)%burned(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%burned(pos_id)
this%variables(i_var)%damaged(pos_id) = recipient_fuse_weight * this%variables(i_var)%damaged(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%damaged(pos_id)
this%variables(i_var)%herbivory(pos_id) = recipient_fuse_weight * this%variables(i_var)%herbivory(pos_id) + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%variables(i_var)%herbivory(pos_id)
end do
end do
this%ode_opt_step = recipient_fuse_weight * this%ode_opt_step + &
(1.0_r8-recipient_fuse_weight) * donor_prt_obj%ode_opt_step
return
end subroutine WeightedFusePRTVartypes
! =====================================================================================
subroutine DeallocatePRTVartypes(this)
! ---------------------------------------------------------------------------------
! Unfortunately ... all plants must die. It is sad, but when this happens
! we must also deallocate our memory of them. Man, thats really is sad. Why
! must we also forget them... Well, anyway, any time a plant/cohort
! is deallocated, we must also deallocate all this memory bound in the PARTEH
! data structure. But on the bright side, there will always be new recruits,
! a new generation, to allocate as well. Life must go on.
! I suppose since we are recording their life in the history output, in a way
! we are remembering them. I feel better now.
! ---------------------------------------------------------------------------------
class(prt_vartypes) :: this
integer :: i_var, istat
character(len=255) :: smsg
! Check to see if there is any value in these pools?
! SHould not deallocate if there is any carbon left
if(allocated(this%variables)) then
do i_var = 1, prt_global%num_vars
deallocate( &
this%variables(i_var)%val, &
this%variables(i_var)%val0, &
this%variables(i_var)%net_alloc, &
this%variables(i_var)%turnover, &
this%variables(i_var)%burned, &
this%variables(i_var)%herbivory, &
stat=istat, errmsg=smsg )
if (istat/=0) call endrun(msg='DeallocatePRTVartypes 1 stat/=0:'//trim(smsg)//errMsg(sourcefile, __LINE__))
end do
deallocate(this%variables, stat=istat, errmsg=smsg)
if (istat/=0) call endrun(msg='DeallocatePRTVartypes 2 stat/=0:'//trim(smsg)//errMsg(sourcefile, __LINE__))
end if
if(allocated(this%bc_in))then
deallocate(this%bc_in, stat=istat, errmsg=smsg)
if (istat/=0) call endrun(msg='DeallocatePRTVartypes bc_in stat/=0:'//trim(smsg)//errMsg(sourcefile, __LINE__))
end if
if(allocated(this%bc_out))then
deallocate(this%bc_out, stat=istat, errmsg=smsg)
if (istat/=0) call endrun(msg='DeallocatePRTVartypes bc_out stat/=0:'//trim(smsg)//errMsg(sourcefile, __LINE__))
end if
if(allocated(this%bc_inout))then
deallocate(this%bc_inout, stat=istat, errmsg=smsg)
if (istat/=0) call endrun(msg='DeallocatePRTVartypes bc_inout stat/=0:'//trim(smsg)//errMsg(sourcefile, __LINE__))
end if
return
end subroutine DeallocatePRTVartypes
! =====================================================================================
subroutine ZeroRates(this)
! ---------------------------------------------------------------------------------
! This subroutine zeros all of the rates of change for our variables.
! It also sets the initial value to the current state.
! This allows us to make mass conservation checks, where
! val - val0 = net_alloc + turnover
!
! This subroutine is called each day in FATES, which is the control interval
! that we conserve carbon from the allocation and turnover process.
! ---------------------------------------------------------------------------------
class(prt_vartypes) :: this
integer :: i_var ! Variable index
do i_var = 1, prt_global%num_vars
this%variables(i_var)%val0(:) = this%variables(i_var)%val(:)
this%variables(i_var)%net_alloc(:) = 0.0_r8
this%variables(i_var)%turnover(:) = 0.0_r8
this%variables(i_var)%burned(:) = 0.0_r8
this%variables(i_var)%herbivory(:) = 0.0_r8
this%variables(i_var)%damaged(:) = 0.0_r8
end do
end subroutine ZeroRates
! ====================================================================================
subroutine CheckMassConservation(this,ipft,position_id)
! ---------------------------------------------------------------------------------
! At any time, the sum of fluxes should equal the difference between val and val0.
! This routine loops over all variables and ensures this is true.
! The final argument is any uniqely identifying index that can be used
! to differentiate where in the call sequence a failure in conservation occurs.
! ---------------------------------------------------------------------------------
class(prt_vartypes) :: this
integer, intent(in) :: ipft ! functional type of the plant
integer, intent(in) :: position_id ! Helps to know where
! in the call sequence this was called
integer :: i_var ! Variable index
integer :: i_pos ! Position (coordinate) index
real(r8) :: err ! absolute error [kg]
real(r8) :: rel_err ! error relative to the pool's size [kg]
do i_var = 1, prt_global%num_vars
do i_pos = 1, prt_global%state_descriptor(i_var)%num_pos
err = abs((this%variables(i_var)%val(i_pos) - this%variables(i_var)%val0(i_pos)) - &
(this%variables(i_var)%net_alloc(i_pos) &
-this%variables(i_var)%turnover(i_pos) &
-this%variables(i_var)%burned(i_pos) &
-this%variables(i_var)%herbivory(i_pos) &
-this%variables(i_var)%damaged(i_pos)))
if(this%variables(i_var)%val(i_pos) > nearzero ) then
rel_err = err / this%variables(i_var)%val(i_pos)
else
rel_err = 0.0_r8
end if
if( abs(err) > calloc_abs_error ) then
write(fates_log(),*) 'PARTEH mass conservation check failed'
write(fates_log(),*) ' Change in mass over control period should'
write(fates_log(),*) ' always equal the integrated fluxes.'
write(fates_log(),*) ' pft id: ',ipft
write(fates_log(),*) ' position id: ',position_id
write(fates_log(),*) ' organ id: ',prt_global%state_descriptor(i_var)%organ_id
write(fates_log(),*) ' element_id: ',prt_global%state_descriptor(i_var)%element_id