Skip to content

solution

aquacrop.solution

aquacrop.solution.adjust_CCx

adjust_CCx(cc_prev, CCo, CCx, CGC, CDC, dt, tSum, Crop_CanopyDevEnd, Crop_CCx)

Function to adjust CCx value for changes in CGC due to water stress during the growing season

Reference Manual: canopy_cover stress response (pg. 27-33)

Arguments:

cc_prev (float): Canopy Cover at previous timestep.

CCo (float): Fractional canopy cover size at emergence

CCx (float): Maximum canopy cover (fraction of soil cover)

CGC (float): Canopy growth coefficient (fraction per gdd)

CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

dt (float): Time delta of canopy growth (1 calander day or ... gdd)

tSum (float): time since germination (CD or gdd)

Crop_CanopyDevEnd (float): time that Canopy developement ends

Crop_CCx (float): Maximum canopy cover (fraction of soil cover)

Returns:

CCxAdj (float): Adjusted CCx
Source code in aquacrop/solution/adjust_CCx.py 
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
def adjust_CCx(
    cc_prev: float,
    CCo: float,
    CCx: float,
    CGC: float,
    CDC: float,
    dt: float,
    tSum: float,
    Crop_CanopyDevEnd: float,
    Crop_CCx: float
    ) -> float:
    """
    Function to adjust CCx value for changes in CGC due to water stress during the growing season

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=36" target="_blank">Reference Manual: canopy_cover stress response</a> (pg. 27-33)


    Arguments:

        cc_prev (float): Canopy Cover at previous timestep.

        CCo (float): Fractional canopy cover size at emergence

        CCx (float): Maximum canopy cover (fraction of soil cover)

        CGC (float): Canopy growth coefficient (fraction per gdd)

        CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

        dt (float): Time delta of canopy growth (1 calander day or ... gdd)

        tSum (float): time since germination (CD or gdd)

        Crop_CanopyDevEnd (float): time that Canopy developement ends

        Crop_CCx (float): Maximum canopy cover (fraction of soil cover)

    Returns:

        CCxAdj (float): Adjusted CCx



    """

    ## Get time required to reach canopy_cover on previous day ##
    tCCtmp = cc_required_time(cc_prev, CCo, CCx, CGC, CDC, "CGC")

    ## Determine CCx adjusted ##
    if tCCtmp > 0:
        tCCtmp = tCCtmp + (Crop_CanopyDevEnd - tSum) + dt
        CCxAdj = cc_development(CCo, CCx, CGC, CDC, tCCtmp, "Growth", Crop_CCx)
    else:
        CCxAdj = 0

    return CCxAdj

aquacrop.solution.aeration_stress

aeration_stress(NewCond_AerDays, Crop_LagAer, thRZ)

Function to calculate aeration stress coefficient

Reference Manual: aeration stress (pg. 89-90)

Arguments:

NewCond_AerDays (int): number aeration stress days

Crop_LagAer (int): lag days before aeration stress

thRZ (NamedTuple): object that contains information on the total water in the root zone

Returns:

Ksa_Aer (float): aeration stress coefficient

NewCond_AerDays (float): updated aer days
Source code in aquacrop/solution/aeration_stress.py 
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
@cc.export("aeration_stress", (f8,f8,thRZNT_type_sig))
def aeration_stress(
    NewCond_AerDays: float,
    Crop_LagAer: float,
    thRZ: NamedTuple,
    ) -> Tuple[float, float]:
    """
    Function to calculate aeration stress coefficient

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=90" target="_blank">Reference Manual: aeration stress</a> (pg. 89-90)


    Arguments:

        NewCond_AerDays (int): number aeration stress days

        Crop_LagAer (int): lag days before aeration stress

        thRZ (NamedTuple): object that contains information on the total water in the root zone


    Returns:

        Ksa_Aer (float): aeration stress coefficient

        NewCond_AerDays (float): updated aer days



    """

    ## Determine aeration stress (root zone) ##
    if thRZ.Act > thRZ.Aer:
        # Calculate aeration stress coefficient
        if NewCond_AerDays < Crop_LagAer:
            stress = 1 - ((thRZ.S - thRZ.Act) / (thRZ.S - thRZ.Aer))
            Ksa_Aer = 1 - ((NewCond_AerDays / 3) * stress)
        elif NewCond_AerDays >= Crop_LagAer:
            Ksa_Aer = (thRZ.S - thRZ.Act) / (thRZ.S - thRZ.Aer)

        # Increment aeration days counter
        NewCond_AerDays = NewCond_AerDays + 1
        if NewCond_AerDays > Crop_LagAer:
            NewCond_AerDays = Crop_LagAer

    else:
        # Set aeration stress coefficient to one (no stress value)
        Ksa_Aer = 1
        # Reset aeration days counter
        NewCond_AerDays = 0

    return Ksa_Aer, NewCond_AerDays

aquacrop.solution.biomass_accumulation

biomass_accumulation(Crop, NewCond_DAP, NewCond_DelayedCDs, NewCond_HIref, NewCond_PctLagPhase, NewCond_B, NewCond_B_NS, Tr, TrPot, et0, growing_season)

Function to calculate biomass accumulation

Reference Manual: biomass accumulaiton (pg. 98-108)

Arguments:

Crop (NamedTuple): Crop object

NewCond_DAP (int): days since planting

NewCond_DelayedCDs (int): Delayed calendar days

NewCond_HIref (float): reference harvest index

NewCond_PctLagPhase (float): percentage of way through early HI development stage

NewCond_B (float): Current biomass growth

NewCond_B_NS (float): current no stress biomass growth

TrPot (float): Daily crop transpiration

TrPot (float): Daily potential transpiration

et0 (float): Daily reference evapotranspiration

growing_season (bool): is Growing season? (True, False)

Returns:

NewCond_B (float): new biomass growth

NewCond_B_NS (float): new (No stress) biomass growth
Source code in aquacrop/solution/biomass_accumulation.py 
 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
@cc.export("biomass_accumulation", (CropStructNT_type_sig,i8,i8,f8,f8,f8,f8,f8,f8,f8,b1))
def biomass_accumulation(
    Crop: NamedTuple,
    NewCond_DAP: int,
    NewCond_DelayedCDs: int,
    NewCond_HIref: float,
    NewCond_PctLagPhase: float,
    NewCond_B: float,
    NewCond_B_NS: float,
    Tr: float,
    TrPot: float,
    et0: float,
    growing_season: bool,
    ) -> Tuple[float, float]:
    """
    Function to calculate biomass accumulation

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=107" target="_blank">Reference Manual: biomass accumulaiton</a> (pg. 98-108)


    Arguments:

        Crop (NamedTuple): Crop object

        NewCond_DAP (int): days since planting

        NewCond_DelayedCDs (int): Delayed calendar days

        NewCond_HIref (float): reference harvest index

        NewCond_PctLagPhase (float): percentage of way through early HI development stage

        NewCond_B (float): Current biomass growth

        NewCond_B_NS (float): current no stress biomass growth

        TrPot (float): Daily crop transpiration

        TrPot (float): Daily potential transpiration

        et0 (float): Daily reference evapotranspiration

        growing_season (bool): is Growing season? (True, False)

    Returns:

        NewCond_B (float): new biomass growth

        NewCond_B_NS (float): new (No stress) biomass growth


    """

    ## Store initial conditions in a new structure for updating ##
    # NewCond = InitCond

    ## Calculate biomass accumulation (if in growing season) ##
    if growing_season == True:
        # Get time for harvest index build-up
        HIt = NewCond_DAP - NewCond_DelayedCDs - Crop.HIstartCD - 1

        if ((Crop.CropType == 2) or (Crop.CropType == 3)) and (NewCond_HIref > 0):
            # Adjust WP for reproductive stage
            if Crop.Determinant == 1:
                fswitch = NewCond_PctLagPhase / 100
            else:
                if HIt < (Crop.YldFormCD / 3):
                    fswitch = HIt / (Crop.YldFormCD / 3)
                else:
                    fswitch = 1

            WPadj = Crop.WP * (1 - (1 - Crop.WPy / 100) * fswitch)
        else:
            WPadj = Crop.WP

        # Adjust WP for CO2 effects
        WPadj = WPadj * Crop.fCO2

        # Calculate biomass accumulation on current day
        # No water stress
        dB_NS = WPadj * (TrPot / et0)
        # With water stress
        dB = WPadj * (Tr / et0)
        if np.isnan(dB) == True:
            dB = 0

        # Update biomass accumulation
        NewCond_B = NewCond_B + dB
        NewCond_B_NS = NewCond_B_NS + dB_NS
    else:
        # No biomass accumulation outside of growing season
        NewCond_B = 0
        NewCond_B_NS = 0

    return (NewCond_B,
            NewCond_B_NS)

aquacrop.solution.canopy_cover

canopy_cover(Crop, prof, Soil_zTop, InitCond, gdd, et0, growing_season)

Function to simulate canopy growth/decline

Reference Manual: canopy_cover equations (pg. 21-33)

Arguments:

Crop (CropStructNT): NamedTuple of Crop object

prof (SoilProfileNT): NamedTuple of SoilProfile object

Soil_zTop (float): top soil depth

InitCond (InitialCondition): InitCond object

gdd (float): Growing Degree Days

et0 (float): reference evapotranspiration

growing_season (bool): is it currently within the growing season (True, Flase)

Returns:

NewCond (InitialCondition): updated InitCond object
Source code in aquacrop/solution/canopy_cover.py 
 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
def canopy_cover(
    Crop: "CropStructNT",
    prof: "SoilProfileNT",
    Soil_zTop: float,
    InitCond: "InitialCondition",
    gdd: float,
    et0: float,
    growing_season: bool,
    ):
    # def CCCrop,Soil_Profile,Soil_zTop,InitCond,gdd,et0,growing_season):

    """
    Function to simulate canopy growth/decline

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=30" target="_blank">Reference Manual: canopy_cover equations</a> (pg. 21-33)


    Arguments:

        Crop (CropStructNT): NamedTuple of Crop object

        prof (SoilProfileNT): NamedTuple of SoilProfile object

        Soil_zTop (float): top soil depth

        InitCond (InitialCondition): InitCond object

        gdd (float): Growing Degree Days

        et0 (float): reference evapotranspiration

        growing_season (bool): is it currently within the growing season (True, Flase)

    Returns:

        NewCond (InitialCondition): updated InitCond object


    """

    # Function to simulate canopy growth/decline

    InitCond_CC_NS = InitCond.canopy_cover_ns
    InitCond_CC = InitCond.canopy_cover
    InitCond_ProtectedSeed = InitCond.protected_seed
    InitCond_CCxAct = InitCond.ccx_act
    InitCond_CropDead = InitCond.crop_dead
    InitCond_tEarlySen = InitCond.t_early_sen
    InitCond_CCxW = InitCond.ccx_w

    ## Store initial conditions in a new structure for updating ##
    NewCond = InitCond
    NewCond.cc_prev = InitCond.canopy_cover

    ## Calculate canopy development (if in growing season) ##
    if growing_season == True:
        # Calculate root zone water content
        taw = TAW()
        root_zone_depletion = Dr()
        # thRZ = RootZoneWater()
        _, root_zone_depletion.Zt, root_zone_depletion.Rz, taw.Zt, taw.Rz, _,_,_,_,_,_ = root_zone_water(
            prof,
            float(NewCond.z_root),
            NewCond.th,
            Soil_zTop,
            float(Crop.Zmin),
            Crop.Aer,
        )

        # _,root_zone_depletion,taw,_ = root_zone_water(Soil_Profile,float(NewCond.z_root),NewCond.th,Soil_zTop,float(Crop.Zmin),Crop.Aer)
        # Check whether to use root zone or top soil depletions for calculating
        # water stress
        if (root_zone_depletion.Rz / taw.Rz) <= (root_zone_depletion.Zt / taw.Zt):
            # Root zone is wetter than top soil, so use root zone value
            root_zone_depletion = root_zone_depletion.Rz
            taw = taw.Rz
        else:
            # Top soil is wetter than root zone, so use top soil values
            root_zone_depletion = root_zone_depletion.Zt
            taw = taw.Zt

        # Determine if water stress is occurring
        beta = True
        water_stress_coef = Ksw()
        water_stress_coef.exp, water_stress_coef.sto, water_stress_coef.sen, water_stress_coef.pol, water_stress_coef.sto_lin = water_stress(
            Crop.p_up,
            Crop.p_lo,
            Crop.ETadj,
            Crop.beta,
            Crop.fshape_w,
            NewCond.t_early_sen,
            root_zone_depletion,
            taw,
            et0,
            beta,
        )

        # water_stress(Crop, NewCond, root_zone_depletion, taw, et0, beta)

        # Get canopy cover growth time
        if Crop.CalendarType == 1:
            dtCC = 1
            tCCadj = NewCond.dap - NewCond.delayed_cds
        elif Crop.CalendarType == 2:
            dtCC = gdd
            tCCadj = NewCond.gdd_cum - NewCond.delayed_gdds

        ## Canopy development (potential) ##
        if (tCCadj < Crop.Emergence) or (round(tCCadj) > Crop.Maturity):
            # No canopy development before emergence/germination or after
            # maturity
            NewCond.canopy_cover_ns = 0
        elif tCCadj < Crop.CanopyDevEnd:
            # Canopy growth can occur
            if InitCond_CC_NS <= Crop.CC0:
                # Very small initial canopy_cover.
                NewCond.canopy_cover_ns = Crop.CC0 * np.exp(Crop.CGC * dtCC)
                # print(Crop.CC0,np.exp(Crop.CGC*dtCC))
            else:
                # Canopy growing
                tmp_tCC = tCCadj - Crop.Emergence
                NewCond.canopy_cover_ns = cc_development(
                    Crop.CC0, 0.98 * Crop.CCx, Crop.CGC, Crop.CDC, tmp_tCC, "Growth", Crop.CCx
                )

            # Update maximum canopy cover size in growing season
            NewCond.ccx_act_ns = NewCond.canopy_cover_ns
        elif tCCadj > Crop.CanopyDevEnd:
            # No more canopy growth is possible or canopy in decline
            # Set CCx for calculation of withered canopy effects
            NewCond.ccx_w_ns = NewCond.ccx_act_ns
            if tCCadj < Crop.Senescence:
                # Mid-season stage - no canopy growth
                NewCond.canopy_cover_ns = InitCond_CC_NS
                # Update maximum canopy cover size in growing season
                NewCond.ccx_act_ns = NewCond.canopy_cover_ns
            else:
                # Late-season stage - canopy decline
                tmp_tCC = tCCadj - Crop.Senescence
                NewCond.canopy_cover_ns = cc_development(
                    Crop.CC0,
                    NewCond.ccx_act_ns,
                    Crop.CGC,
                    Crop.CDC,
                    tmp_tCC,
                    "Decline",
                    NewCond.ccx_act_ns,
                )

        ## Canopy development (actual) ##
        if (tCCadj < Crop.Emergence) or (round(tCCadj) > Crop.Maturity):
            # No canopy development before emergence/germination or after
            # maturity
            NewCond.canopy_cover = 0
            NewCond.cc0_adj = Crop.CC0
        elif tCCadj < Crop.CanopyDevEnd:
            # Canopy growth can occur
            if InitCond_CC <= NewCond.cc0_adj or (
                (InitCond_ProtectedSeed == True) and (InitCond_CC <= (1.25 * NewCond.cc0_adj))
            ):
                # Very small initial canopy_cover or seedling in protected phase of
                # growth. In this case, assume no leaf water expansion stress
                if InitCond_ProtectedSeed == True:
                    tmp_tCC = tCCadj - Crop.Emergence
                    NewCond.canopy_cover = cc_development(
                        Crop.CC0, Crop.CCx, Crop.CGC, Crop.CDC, tmp_tCC, "Growth", Crop.CCx
                    )
                    # Check if seed protection should be turned off
                    if NewCond.canopy_cover > (1.25 * NewCond.cc0_adj):
                        # Turn off seed protection - lead expansion stress can
                        # occur on future time steps.
                        NewCond.protected_seed = False

                else:
                    NewCond.canopy_cover = NewCond.cc0_adj * np.exp(Crop.CGC * dtCC)

            else:
                # Canopy growing

                if InitCond_CC < (0.9799 * Crop.CCx):
                    # Adjust canopy growth coefficient for leaf expansion water
                    # stress effects
                    CGCadj = Crop.CGC * water_stress_coef.exp
                    if CGCadj > 0:

                        # Adjust CCx for change in CGC
                        CCXadj = adjust_CCx(
                            InitCond_CC,
                            NewCond.cc0_adj,
                            Crop.CCx,
                            CGCadj,
                            Crop.CDC,
                            dtCC,
                            tCCadj,
                            Crop.CanopyDevEnd,
                            Crop.CCx,
                        )
                        if CCXadj < 0:

                            NewCond.canopy_cover = InitCond_CC
                        elif abs(InitCond_CC - (0.9799 * Crop.CCx)) < 0.001:

                            # Approaching maximum canopy cover size
                            tmp_tCC = tCCadj - Crop.Emergence
                            NewCond.canopy_cover = cc_development(
                                Crop.CC0, Crop.CCx, Crop.CGC, Crop.CDC, tmp_tCC, "Growth", Crop.CCx
                            )
                        else:

                            # Determine time required to reach canopy_cover on previous,
                            # day, given CGCAdj value
                            tReq = cc_required_time(
                                InitCond_CC, NewCond.cc0_adj, CCXadj, CGCadj, Crop.CDC, "CGC"
                            )
                            if tReq > 0:

                                # Calclate gdd's for canopy growth
                                tmp_tCC = tReq + dtCC
                                # Determine new canopy size
                                NewCond.canopy_cover = cc_development(
                                    NewCond.cc0_adj,
                                    CCXadj,
                                    CGCadj,
                                    Crop.CDC,
                                    tmp_tCC,
                                    "Growth",
                                    Crop.CCx,
                                )
                                # print(NewCond.dap,CCXadj,tReq)

                            else:
                                # No canopy growth
                                NewCond.canopy_cover = InitCond_CC

                    else:

                        # No canopy growth
                        NewCond.canopy_cover = InitCond_CC
                        # Update CC0
                        if NewCond.canopy_cover > NewCond.cc0_adj:
                            NewCond.cc0_adj = Crop.CC0
                        else:
                            NewCond.cc0_adj = NewCond.canopy_cover

                else:
                    # Canopy approaching maximum size
                    tmp_tCC = tCCadj - Crop.Emergence
                    NewCond.canopy_cover = cc_development(
                        Crop.CC0, Crop.CCx, Crop.CGC, Crop.CDC, tmp_tCC, "Growth", Crop.CCx
                    )
                    NewCond.cc0_adj = Crop.CC0

            if NewCond.canopy_cover > InitCond_CCxAct:
                # Update actual maximum canopy cover size during growing season
                NewCond.ccx_act = NewCond.canopy_cover

        elif tCCadj > Crop.CanopyDevEnd:

            # No more canopy growth is possible or canopy is in decline
            if tCCadj < Crop.Senescence:
                # Mid-season stage - no canopy growth
                NewCond.canopy_cover = InitCond_CC
                if NewCond.canopy_cover > InitCond_CCxAct:
                    # Update actual maximum canopy cover size during growing
                    # season
                    NewCond.ccx_act = NewCond.canopy_cover

            else:
                # Late-season stage - canopy decline
                # Adjust canopy decline coefficient for difference between actual
                # and potential CCx
                CDCadj = Crop.CDC * ((NewCond.ccx_act + 2.29) / (Crop.CCx + 2.29))
                # Determine new canopy size
                tmp_tCC = tCCadj - Crop.Senescence
                NewCond.canopy_cover = cc_development(
                    NewCond.cc0_adj,
                    NewCond.ccx_act,
                    Crop.CGC,
                    CDCadj,
                    tmp_tCC,
                    "Decline",
                    NewCond.ccx_act,
                )

            # Check for crop growth termination
            if (NewCond.canopy_cover < 0.001) and (InitCond_CropDead == False):
                # Crop has died
                NewCond.canopy_cover = 0
                NewCond.crop_dead = True

        ## Canopy senescence due to water stress (actual) ##
        if tCCadj >= Crop.Emergence:
            if (tCCadj < Crop.Senescence) or (InitCond_tEarlySen > 0):
                # Check for early canopy senescence  due to severe water
                # stress.
                if (water_stress_coef.sen < 1) and (InitCond_ProtectedSeed == False):

                    # Early canopy senescence
                    NewCond.premat_senes = True
                    if InitCond_tEarlySen == 0:
                        # No prior early senescence
                        NewCond.ccx_early_sen = InitCond_CC

                    # Increment early senescence gdd counter
                    NewCond.t_early_sen = InitCond_tEarlySen + dtCC
                    # Adjust canopy decline coefficient for water stress
                    beta = False

                    water_stress_coef = Ksw()
                    water_stress_coef.exp, water_stress_coef.sto, water_stress_coef.sen, water_stress_coef.pol, water_stress_coef.sto_lin = water_stress(
                        Crop.p_up,
                        Crop.p_lo,
                        Crop.ETadj,
                        Crop.beta,
                        Crop.fshape_w,
                        NewCond.t_early_sen,
                        root_zone_depletion,
                        taw,
                        et0,
                        beta,
                    )

                    # water_stress_coef = water_stress(Crop, NewCond, root_zone_depletion, taw, et0, beta)
                    if water_stress_coef.sen > 0.99999:
                        CDCadj = 0.0001
                    else:
                        CDCadj = (1 - (water_stress_coef.sen ** 8)) * Crop.CDC

                    # Get new canpy cover size after senescence
                    if NewCond.ccx_early_sen < 0.001:
                        CCsen = 0
                    else:
                        # Get time required to reach canopy_cover at end of previous day, given
                        # CDCadj
                        tReq = (np.log(1 + (1 - InitCond_CC / NewCond.ccx_early_sen) / 0.05)) / (
                            (CDCadj * 3.33) / (NewCond.ccx_early_sen + 2.29)
                        )
                        # Calculate gdd's for canopy decline
                        tmp_tCC = tReq + dtCC
                        # Determine new canopy size
                        CCsen = NewCond.ccx_early_sen * (
                            1
                            - 0.05
                            * (
                                np.exp(tmp_tCC * ((CDCadj * 3.33) / (NewCond.ccx_early_sen + 2.29)))
                                - 1
                            )
                        )
                        if CCsen < 0:
                            CCsen = 0

                    # Update canopy cover size
                    if tCCadj < Crop.Senescence:
                        # Limit canopy_cover to CCx
                        if CCsen > Crop.CCx:
                            CCsen = Crop.CCx

                        # canopy_cover cannot be greater than value on previous day
                        NewCond.canopy_cover = CCsen
                        if NewCond.canopy_cover > InitCond_CC:
                            NewCond.canopy_cover = InitCond_CC

                        # Update maximum canopy cover size during growing
                        # season
                        NewCond.ccx_act = NewCond.canopy_cover
                        # Update CC0 if current canopy_cover is less than initial canopy
                        # cover size at planting
                        if NewCond.canopy_cover < Crop.CC0:
                            NewCond.cc0_adj = NewCond.canopy_cover
                        else:
                            NewCond.cc0_adj = Crop.CC0

                    else:
                        # Update canopy_cover to account for canopy cover senescence due
                        # to water stress
                        if CCsen < NewCond.canopy_cover:
                            NewCond.canopy_cover = CCsen

                    # Check for crop growth termination
                    if (NewCond.canopy_cover < 0.001) and (InitCond_CropDead == False):
                        # Crop has died
                        NewCond.canopy_cover = 0
                        NewCond.crop_dead = True

                else:
                    # No water stress
                    NewCond.premat_senes = False
                    if (tCCadj > Crop.Senescence) and (InitCond_tEarlySen > 0):
                        # Rewatering of canopy in late season
                        # Get new values for CCx and CDC
                        tmp_tCC = tCCadj - dtCC - Crop.Senescence
                        CCXadj, CDCadj = update_CCx_CDC(InitCond_CC, Crop.CDC, Crop.CCx, tmp_tCC)
                        NewCond.ccx_act = CCXadj
                        # Get new canopy_cover value for end of current day
                        tmp_tCC = tCCadj - Crop.Senescence
                        NewCond.canopy_cover = cc_development(
                            NewCond.cc0_adj, CCXadj, Crop.CGC, CDCadj, tmp_tCC, "Decline", CCXadj
                        )
                        # Check for crop growth termination
                        if (NewCond.canopy_cover < 0.001) and (InitCond_CropDead == False):
                            NewCond.canopy_cover = 0
                            NewCond.crop_dead = True

                    # Reset early senescence counter
                    NewCond.t_early_sen = 0

                # Adjust CCx for effects of withered canopy
                if NewCond.canopy_cover > InitCond_CCxW:
                    NewCond.ccx_w = NewCond.canopy_cover

        ## Calculate canopy size adjusted for micro-advective effects ##
        # Check to ensure potential canopy_cover is not slightly lower than actual
        if NewCond.canopy_cover_ns < NewCond.canopy_cover:
            NewCond.canopy_cover_ns = NewCond.canopy_cover
            if tCCadj < Crop.CanopyDevEnd:
                NewCond.ccx_act_ns = NewCond.canopy_cover_ns

        # Actual (with water stress)
        NewCond.canopy_cover_adj = (1.72 * NewCond.canopy_cover) - (NewCond.canopy_cover ** 2) + (0.3 * (NewCond.canopy_cover ** 3))
        # Potential (without water stress)
        NewCond.canopy_cover_adj_ns = (
            (1.72 * NewCond.canopy_cover_ns) - (NewCond.canopy_cover_ns ** 2) + (0.3 * (NewCond.canopy_cover_ns ** 3))
        )

    else:
        # No canopy outside growing season - set various values to zero
        NewCond.canopy_cover = 0
        NewCond.canopy_cover_adj = 0
        NewCond.canopy_cover_ns = 0
        NewCond.canopy_cover_adj_ns = 0
        NewCond.ccx_w = 0
        NewCond.ccx_act = 0
        NewCond.ccx_w_ns = 0
        NewCond.ccx_act_ns = 0

    return NewCond

aquacrop.solution.capillary_rise

capillary_rise(prof, Soil_nLayer, Soil_fshape_cr, NewCond, FluxOut, water_table_presence)

Function to calculate capillary rise from a shallow groundwater table

Reference Manual: capillary rise calculations (pg. 52-61)

Arguments:

prof (SoilProfileNT): Soil profile named tuple

Soil_nLayer (int): number of soil layers

Soil_fshape_cr (float): Capillary rise shape factor

NewCond (InitialCondition): InitialCondition object containing model paramaters

FluxOut (numpy.array): Flux of water out of each soil compartment

water_table_presence (int): water_table present (1:yes, 0:no)

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters

CrTot (float): Total Capillary rise
Source code in aquacrop/solution/capillary_rise.py 
 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
def capillary_rise(
    prof: "SoilProfileNT",
    Soil_nLayer: int,
    Soil_fshape_cr: float,
    NewCond: "InitialCondition",
    FluxOut: "ndarray",
    water_table_presence: int,
    ) -> Tuple["InitialCondition", float]:
    """
    Function to calculate capillary rise from a shallow groundwater table


    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=61" target="_blank">Reference Manual: capillary rise calculations</a> (pg. 52-61)


    Arguments:


        prof (SoilProfileNT): Soil profile named tuple

        Soil_nLayer (int): number of soil layers

        Soil_fshape_cr (float): Capillary rise shape factor

        NewCond (InitialCondition): InitialCondition object containing model paramaters

        FluxOut (numpy.array): Flux of water out of each soil compartment

        water_table_presence (int): water_table present (1:yes, 0:no)


    Returns:

        NewCond (InitialCondition): InitCond object containing updated model paramaters

        CrTot (float): Total Capillary rise





    """

    ## Get groundwater table elevation on current day ##
    z_gw = NewCond.z_gw

    ## Calculate capillary rise ##
    if water_table_presence == 0:  # No water table present
        # Capillary rise is zero
        CrTot = 0
    elif water_table_presence == 1:  # Water table present
        # Get maximum capillary rise for bottom compartment
        zBot = prof.dzsum[-1]
        zBotMid = prof.zMid[-1]
        prof = prof
        if (prof.Ksat[-1] > 0) and (z_gw > 0) and ((z_gw - zBotMid) < 4):
            if zBotMid >= z_gw:
                MaxCR = 99
            else:
                MaxCR = np.exp((np.log(z_gw - zBotMid) - prof.bCR[-1]) / prof.aCR[-1])
                if MaxCR > 99:
                    MaxCR = 99

        else:
            MaxCR = 0

        ######################### this needs fixing, will currently break####################

        # # Find top of next soil layer that is not within modelled soil profile
        zTopLayer = 0
        for layeri in np.sort(np.unique(prof.Layer)):
            # Calculate layer thickness
            l_idx = np.argwhere(prof.Layer==layeri).flatten()

            LayThk = prof.dz[l_idx].sum()
            zTopLayer = zTopLayer+LayThk

        # Check for restrictions on upward flow caused by properties of
        # compartments that are not modelled in the soil water balance
        layeri = prof.Layer[-1]

        assert layeri == Soil_nLayer

        while (zTopLayer < z_gw) and (layeri < Soil_nLayer):
            # this needs fixing, will currently break


            compdf = prof.Layer[layeri]
            if (compdf.Ksat > 0) and (z_gw > 0) and ((z_gw-zTopLayer) < 4):
                if zTopLayer >= z_gw:
                    LimCR = 99
                else:
                    LimCR = np.exp((np.log(z_gw-zTopLayer)-compdf.bCR)/compdf.aCR)
                    if LimCR > 99:
                        LimCR = 99

            else:
                LimCR = 0

            if MaxCR > LimCR:
                MaxCR = LimCR

            zTopLayer = zTopLayer+compdf.dz

            layeri = layeri+1 # could be that the increment should be at the end of the loop rather than at the start (bad indexing otherwise)

        #####################################################################################

        # Calculate capillary rise
        compi = len(prof.Comp) - 1  # Start at bottom of root zone
        WCr = 0  # Capillary rise counter
        while (round(MaxCR * 1000) > 0) and (compi > -1) and (round(FluxOut[compi] * 1000) == 0):
            # Proceed upwards until maximum capillary rise occurs, soil surface
            # is reached, or encounter a compartment where downward
            # drainage/infiltration has already occurred on current day
            # Find layer of current compartment
            # Calculate driving force
            if (NewCond.th[compi] >= prof.th_wp[compi]) and (Soil_fshape_cr > 0):
                Df = 1 - (
                    (
                        (NewCond.th[compi] - prof.th_wp[compi])
                        / (NewCond.th_fc_Adj[compi] - prof.th_wp[compi])
                    )
                    ** Soil_fshape_cr
                )
                if Df > 1:
                    Df = 1
                elif Df < 0:
                    Df = 0

            else:
                Df = 1

            # Calculate relative hydraulic conductivity
            thThr = (prof.th_wp[compi] + prof.th_fc[compi]) / 2
            if NewCond.th[compi] < thThr:
                if (NewCond.th[compi] <= prof.th_wp[compi]) or (thThr <= prof.th_wp[compi]):
                    Krel = 0
                else:
                    Krel = (NewCond.th[compi] - prof.th_wp[compi]) / (thThr - prof.th_wp[compi])

            else:
                Krel = 1

            # Check if room is available to store water from capillary rise
            dth = round(NewCond.th_fc_Adj[compi] - NewCond.th[compi],4)

            # Store water if room is available
            if (dth > 0) and ((zBot - prof.dz[compi] / 2) < z_gw):
                dthMax = Krel * Df * MaxCR / (1000 * prof.dz[compi])
                if dth >= dthMax:
                    NewCond.th[compi] = NewCond.th[compi] + dthMax
                    CRcomp = dthMax * 1000 * prof.dz[compi]
                    MaxCR = 0
                else:
                    NewCond.th[compi] = NewCond.th_fc_Adj[compi]
                    CRcomp = dth * 1000 * prof.dz[compi]
                    MaxCR = (Krel * MaxCR) - CRcomp

                WCr = WCr + CRcomp

            # Update bottom elevation of compartment
            zBot = zBot - prof.dz[compi]
            # Update compartment and layer counters
            compi = compi - 1
            # Update restriction on maximum capillary rise
            if compi > -1:

                zBotMid = zBot - (prof.dz[compi] / 2)
                if (prof.Ksat[compi] > 0) and (z_gw > 0) and ((z_gw - zBotMid) < 4):
                    if zBotMid >= z_gw:
                        LimCR = 99
                    else:
                        LimCR = np.exp((np.log(z_gw - zBotMid) - prof.bCR[compi]) / prof.aCR[compi])
                        if LimCR > 99:
                            LimCR = 99

                else:
                    LimCR = 0

                if MaxCR > LimCR:
                    MaxCR = LimCR

        # Store total depth of capillary rise
        CrTot = WCr

    return NewCond, CrTot

aquacrop.solution.cc_development

cc_development(CCo, CCx, CGC, CDC, dt, Mode, CCx0)

Function to calculate canopy cover development by end of the current simulation day

Reference Manual: canopy_cover devlopment (pg. 21-24)

Arguments:

CCo (float): Fractional canopy cover size at emergence

CCx (float): Maximum canopy cover (fraction of soil cover)

CGC (float): Canopy growth coefficient (fraction per gdd)

CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

dt (float): Time delta of canopy growth (1 calander day or ... gdd)

Mode (str): stage of Canopy developement (Growth or Decline)

CCx0 (float): Maximum canopy cover (fraction of soil cover)

Returns:

canopy_cover (float): Canopy Cover
Source code in aquacrop/solution/cc_development.py 
 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
@cc.export("cc_development", "f8(f8,f8,f8,f8,f8,unicode_type,f8)")
def cc_development(
    CCo: float,
    CCx: float,
    CGC: float,
    CDC: float,
    dt: float,
    Mode: str,
    CCx0: float,
    ) -> float:

    """
    Function to calculate canopy cover development by end of the current simulation day

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=30" target="_blank">Reference Manual: canopy_cover devlopment</a> (pg. 21-24)


    Arguments:


        CCo (float): Fractional canopy cover size at emergence

        CCx (float): Maximum canopy cover (fraction of soil cover)

        CGC (float): Canopy growth coefficient (fraction per gdd)

        CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

        dt (float): Time delta of canopy growth (1 calander day or ... gdd)

        Mode (str): stage of Canopy developement (Growth or Decline)

        CCx0 (float): Maximum canopy cover (fraction of soil cover)

    Returns:

        canopy_cover (float): Canopy Cover



    """

    ## Initialise output ##

    ## Calculate new canopy cover ##
    if Mode == "Growth":
        # Calculate canopy growth
        # Exponential growth stage
        canopy_cover = CCo * np.exp(CGC * dt)
        if canopy_cover > (CCx / 2):
            # Exponential decay stage
            canopy_cover = CCx - 0.25 * (CCx / CCo) * CCx * np.exp(-CGC * dt)

        # Limit canopy_cover to CCx
        if canopy_cover > CCx:
            canopy_cover = CCx

    elif Mode == "Decline":
        # Calculate canopy decline
        if CCx < 0.001:
            canopy_cover = 0
        else:
            canopy_cover = CCx * (
                1
                - 0.05
                * (np.exp(dt * CDC * 3.33 * ((CCx + 2.29) / (CCx0 + 2.29)) / (CCx + 2.29)) - 1)
            )

    ## Limit canopy cover to between 0 and 1 ##
    if canopy_cover > 1:
        canopy_cover = 1
    elif canopy_cover < 0:
        canopy_cover = 0

    return canopy_cover

aquacrop.solution.cc_required_time

cc_required_time(cc_prev, CCo, CCx, CGC, CDC, Mode)

Function to find time required to reach canopy_cover at end of previous day, given current CGC or CDC

Reference Manual: canopy_cover devlopment (pg. 21-24)

Arguments:

cc_prev (float): Canopy Cover at previous timestep.

CCo (float): Fractional canopy cover size at emergence

CCx (float): Maximum canopy cover (fraction of soil cover)

CGC (float): Canopy growth coefficient (fraction per gdd)

CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

Mode (str): Canopy growth/decline coefficient (fraction per gdd/calendar day)

Returns:

tReq (float): time required to reach canopy_cover at end of previous day
Source code in aquacrop/solution/cc_required_time.py 
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
@cc.export("cc_required_time", "f8(f8,f8,f8,f8,f8,unicode_type)")
def cc_required_time(
    cc_prev: float,
    CCo: float,
    CCx: float,
    CGC: float,
    CDC: float,
    Mode: str,
    ) -> float:
    """
    Function to find time required to reach canopy_cover at end of previous day, given current CGC or CDC

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=30" target="_blank">Reference Manual: canopy_cover devlopment</a> (pg. 21-24)


    Arguments:

        cc_prev (float): Canopy Cover at previous timestep.

        CCo (float): Fractional canopy cover size at emergence

        CCx (float): Maximum canopy cover (fraction of soil cover)

        CGC (float): Canopy growth coefficient (fraction per gdd)

        CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

        Mode (str): Canopy growth/decline coefficient (fraction per gdd/calendar day)


    Returns:

        tReq (float): time required to reach canopy_cover at end of previous day





    """

    ## Get CGC and/or time (gdd or CD) required to reach canopy_cover on previous day ##
    if Mode == "CGC":
        if cc_prev <= (CCx / 2):

            # print(cc_prev,CCo,(tSum-dt),tSum,dt)
            CGCx = np.log(cc_prev / CCo)
            # print(np.log(cc_prev/CCo),(tSum-dt),CGCx)
        else:
            # print(CCx,CCo,cc_prev)
            CGCx = np.log((0.25 * CCx * CCx / CCo) / (CCx - cc_prev))

        tReq = CGCx / CGC

    elif Mode == "CDC":
        tReq = (np.log(1 + (1 - cc_prev / CCx) / 0.05)) / (CDC / CCx)

    return tReq

aquacrop.solution.check_groundwater_table

check_groundwater_table(prof, NewCond_zGW, NewCond_th, NewCond_th_fc_Adj, water_table_presence, z_gw)

Function to check for presence of a groundwater table, and, if present, to adjust compartment water contents and field capacities where necessary

Reference manual: water table adjustment equations (pg. 52-57)

Arguments:

prof (SoilProfileNT): soil profile paramaters

NewCond_zGW (float): groundwater depth

NewCond_th (ndarray): water content in each soil commpartment

NewCond_th_fc_Adj (ndarray): adjusted water content at field capacity

water_table_presence (int): indicates if water table is present or not

z_gw (float): groundwater depth

Returns:

NewCond_th_fc_Adj (ndarray): adjusted water content at field capacity

thfcAdj (ndarray): adjusted water content at field capacity
Source code in aquacrop/solution/check_groundwater_table.py 
 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
@cc.export("check_groundwater_table", (SoilProfileNT_typ_sig,f8,f8[:],f8[:],i8,f8))
def check_groundwater_table(
    prof: "SoilProfileNT",
    NewCond_zGW: float,
    NewCond_th: "ndarray",
    NewCond_th_fc_Adj: "ndarray",
    water_table_presence: int,
    z_gw: float,
) -> "ndarray":
    """
    Function to check for presence of a groundwater table, and, if present,
    to adjust compartment water contents and field capacities where necessary

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=61" target="_blank">Reference manual: water table adjustment equations</a> (pg. 52-57)


    Arguments:

        prof (SoilProfileNT): soil profile paramaters

        NewCond_zGW (float): groundwater depth

        NewCond_th (ndarray): water content in each soil commpartment

        NewCond_th_fc_Adj (ndarray): adjusted water content at field capacity

        water_table_presence (int): indicates if water table is present or not

        z_gw (float): groundwater depth

    Returns:

        NewCond_th_fc_Adj (ndarray): adjusted water content at field capacity

        thfcAdj (ndarray): adjusted water content at field capacity



    """

    ## Perform calculations (if variable water table is present) ##
    if water_table_presence == 1:

        # Update groundwater conditions for current day
        NewCond_zGW = z_gw

        # Find compartment mid-points
        zMid = prof.zMid

        # Check if water table is within modelled soil profile
        if NewCond_zGW >= 0:
            if len(zMid[zMid >= NewCond_zGW]) == 0:
                NewCond_WTinSoil = False
            else:
                NewCond_WTinSoil = True

        # If water table is in soil profile, adjust water contents
        # if NewCond_WTinSoil == True:
        #     idx = np.argwhere(zMid >= NewCond_zGW).flatten()[0]
        #     for ii in range(idx, len(prof.Comp)):
        #         NewCond_th[ii] = prof.th_s[ii]

        # Adjust compartment field capacity
        compi = len(prof.Comp) - 1
        thfcAdj = np.zeros(compi + 1)
        # Find thFCadj for all compartments
        while compi >= 0:
            if prof.th_fc[compi] <= 0.1:
                Xmax = 1
            else:
                if prof.th_fc[compi] >= 0.3:
                    Xmax = 2
                else:
                    pF = 2 + 0.3 * (prof.th_fc[compi] - 0.1) / 0.2
                    Xmax = (np.exp(pF * np.log(10))) / 100

            if (NewCond_zGW < 0) or ((NewCond_zGW - zMid[compi]) >= Xmax):
                for ii in range(compi + 1):

                    thfcAdj[ii] = prof.th_fc[ii]

                compi = -1
            else:
                if prof.th_fc[compi] >= prof.th_s[compi]:
                    thfcAdj[compi] = prof.th_fc[compi]
                else:
                    if zMid[compi] >= NewCond_zGW:
                        thfcAdj[compi] = prof.th_s[compi]
                    else:
                        dV = prof.th_s[compi] - prof.th_fc[compi]
                        dFC = (dV / (Xmax * Xmax)) * ((zMid[compi] - (NewCond_zGW - Xmax)) ** 2)
                        thfcAdj[compi] = prof.th_fc[compi] + dFC

                compi = compi - 1

        # Store adjusted field capacity values
        NewCond_th_fc_Adj = thfcAdj
        # prof.th_fc_Adj = thfcAdj
        return (NewCond_th_fc_Adj, NewCond_WTinSoil, NewCond_zGW)

    return (NewCond_th_fc_Adj, None, None)

aquacrop.solution.drainage

drainage(prof, th_init, th_fc_Adj_init)

Function to redistribute stored soil water

Reference Manual: drainage calculations (pg. 42-65)

Arguments:

prof (SoilProfile): jit class object object containing soil paramaters

th_init (numpy.array): initial water content

th_fc_Adj_init (numpy.array): adjusted water content at field capacity

Returns:

thnew (numpy.array): updated water content in each compartment

DeepPerc (float): Total Deep Percolation

FluxOut (numpy.array): flux of water out of each compartment
Source code in aquacrop/solution/drainage.py 
 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
@cc.export("drainage", (SoilProfileNT_typ_sig, f8[:], f8[:]))
def drainage(
    prof: "SoilProfileNT",
    th_init: "ndarray",
    th_fc_Adj_init: "ndarray",
    ) -> Tuple["ndarray", float, float]:
    """
    Function to redistribute stored soil water


    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=51" target="_blank">Reference Manual: drainage calculations</a> (pg. 42-65)


    Arguments:

        prof (SoilProfile): jit class object object containing soil paramaters

        th_init (numpy.array): initial water content

        th_fc_Adj_init (numpy.array): adjusted water content at field capacity


    Returns:

        thnew (numpy.array): updated water content in each compartment

        DeepPerc (float): Total Deep Percolation

        FluxOut (numpy.array): flux of water out of each compartment




    """

    # Store initial conditions in new structure for updating %%
    #     NewCond = InitCond

    #     th_init = InitCond.th
    #     th_fc_Adj_init = InitCond.th_fc_Adj

    #  Preallocate arrays %%
    thnew = np.zeros(th_init.shape[0])
    FluxOut = np.zeros(th_init.shape[0])

    # Initialise counters and states %%
    drainsum = 0

    # Calculate drainage and updated water contents %%
    for ii in range(th_init.shape[0]):
        # Specify layer for compartment
        cth_fc = prof.th_fc[ii]
        cth_s = prof.th_s[ii]
        ctau = prof.tau[ii]
        cdz = prof.dz[ii]
        cdzsum = prof.dzsum[ii]
        cKsat = prof.Ksat[ii]

        # Calculate drainage ability of compartment ii
        if th_init[ii] <= th_fc_Adj_init[ii]:
            dthdt = 0

        elif th_init[ii] >= cth_s:
            dthdt = ctau * (cth_s - cth_fc)

            if (th_init[ii] - dthdt) < th_fc_Adj_init[ii]:
                dthdt = th_init[ii] - th_fc_Adj_init[ii]

        else:
            dthdt = (
                ctau
                * (cth_s - cth_fc)
                * ((np.exp(th_init[ii] - cth_fc) - 1) / (np.exp(cth_s - cth_fc) - 1))
            )

            if (th_init[ii] - dthdt) < th_fc_Adj_init[ii]:
                dthdt = th_init[ii] - th_fc_Adj_init[ii]

        # Drainage from compartment ii (mm)
        draincomp = dthdt * cdz * 1000

        # Check drainage ability of compartment ii against cumulative drainage
        # from compartments above
        excess = 0
        prethick = cdzsum - cdz
        drainmax = dthdt * 1000 * prethick
        if drainsum <= drainmax:
            drainability = True
        else:
            drainability = False

        # Drain compartment ii
        if drainability == True:
            # No storage needed. Update water content in compartment ii
            thnew[ii] = th_init[ii] - dthdt

            # Update cumulative drainage (mm)
            drainsum = drainsum + draincomp

            # Restrict cumulative drainage to saturated hydraulic
            # conductivity and adjust excess drainage flow
            if drainsum > cKsat:
                excess = excess + drainsum - cKsat
                drainsum = cKsat

        elif drainability == False:
            # Storage is needed
            dthdt = drainsum / (1000 * prethick)

            # Calculate value of theta (thX) needed to provide a
            # drainage ability equal to cumulative drainage
            if dthdt <= 0:
                thX = th_fc_Adj_init[ii]
            elif ctau > 0:
                A = 1 + (
                    (dthdt * (np.exp(cth_s - cth_fc) - 1)) / (ctau * (cth_s - cth_fc))
                )
                thX = cth_fc + np.log(A)
                if thX < th_fc_Adj_init[ii]:
                    thX = th_fc_Adj_init[ii]

            else:
                thX = cth_s + 0.01

            # Check thX against hydraulic properties of current soil layer
            if thX <= cth_s:
                # Increase compartment ii water content with cumulative
                # drainage
                thnew[ii] = th_init[ii] + (drainsum / (1000 * cdz))
                # Check updated water content against thX
                if thnew[ii] > thX:
                    # Cumulative drainage is the drainage difference
                    # between theta_x and new theta plus drainage ability
                    # at theta_x.
                    drainsum = (thnew[ii] - thX) * 1000 * cdz
                    # Calculate drainage ability for thX
                    if thX <= th_fc_Adj_init[ii]:
                        dthdt = 0
                    elif thX >= cth_s:
                        dthdt = ctau * (cth_s - cth_fc)
                        if (thX - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thX - th_fc_Adj_init[ii]

                    else:
                        dthdt = (
                            ctau
                            * (cth_s - cth_fc)
                            * (
                                (np.exp(thX - cth_fc) - 1)
                                / (np.exp(cth_s - cth_fc) - 1)
                            )
                        )

                        if (thX - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thX - th_fc_Adj_init[ii]

                    # Update drainage total
                    drainsum = drainsum + (dthdt * 1000 * cdz)
                    # Restrict cumulative drainage to saturated hydraulic
                    # conductivity and adjust excess drainage flow
                    if drainsum > cKsat:
                        excess = excess + drainsum - cKsat
                        drainsum = cKsat

                    # Update water content
                    thnew[ii] = thX - dthdt

                elif thnew[ii] > th_fc_Adj_init[ii]:
                    # Calculate drainage ability for updated water content
                    if thnew[ii] <= th_fc_Adj_init[ii]:
                        dthdt = 0
                    elif thnew[ii] >= cth_s:
                        dthdt = ctau * (cth_s - cth_fc)
                        if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thnew[ii] - th_fc_Adj_init[ii]

                    else:
                        dthdt = (
                            ctau
                            * (cth_s - cth_fc)
                            * (
                                (np.exp(thnew[ii] - cth_fc) - 1)
                                / (np.exp(cth_s - cth_fc) - 1)
                            )
                        )
                        if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thnew[ii] - th_fc_Adj_init[ii]

                    # Update water content in compartment ii
                    thnew[ii] = thnew[ii] - dthdt
                    # Update cumulative drainage
                    drainsum = dthdt * 1000 * cdz
                    # Restrict cumulative drainage to saturated hydraulic
                    # conductivity and adjust excess drainage flow
                    if drainsum > cKsat:
                        excess = excess + drainsum - cKsat
                        drainsum = cKsat

                else:
                    # Drainage and cumulative drainage are zero as water
                    # content has not risen above field capacity in
                    # compartment ii.
                    drainsum = 0

            elif thX > cth_s:
                # Increase water content in compartment ii with cumulative
                # drainage from above
                thnew[ii] = th_init[ii] + (drainsum / (1000 * cdz))
                # Check new water content against hydraulic properties of soil
                # layer
                if thnew[ii] <= cth_s:
                    if thnew[ii] > th_fc_Adj_init[ii]:
                        # Calculate new drainage ability
                        if thnew[ii] <= th_fc_Adj_init[ii]:
                            dthdt = 0
                        elif thnew[ii] >= cth_s:
                            dthdt = ctau * (cth_s - cth_fc)
                            if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                                dthdt = thnew[ii] - th_fc_Adj_init[ii]

                        else:
                            dthdt = (
                                ctau
                                * (cth_s - cth_fc)
                                * (
                                    (np.exp(thnew[ii] - cth_fc) - 1)
                                    / (np.exp(cth_s - cth_fc) - 1)
                                )
                            )
                            if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                                dthdt = thnew[ii] - th_fc_Adj_init[ii]

                        # Update water content in compartment ii
                        thnew[ii] = thnew[ii] - dthdt
                        # Update cumulative drainage
                        drainsum = dthdt * 1000 * cdz
                        # Restrict cumulative drainage to saturated hydraulic
                        # conductivity and adjust excess drainage flow
                        if drainsum > cKsat:
                            excess = excess + drainsum - cKsat
                            drainsum = cKsat

                    else:
                        drainsum = 0

                elif thnew[ii] > cth_s:
                    # Calculate excess drainage above saturation
                    excess = (thnew[ii] - cth_s) * 1000 * cdz
                    # Calculate drainage ability for updated water content
                    if thnew[ii] <= th_fc_Adj_init[ii]:
                        dthdt = 0
                    elif thnew[ii] >= cth_s:
                        dthdt = ctau * (cth_s - cth_fc)
                        if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thnew[ii] - th_fc_Adj_init[ii]

                    else:
                        dthdt = (
                            ctau
                            * (cth_s - cth_fc)
                            * (
                                (np.exp(thnew[ii] - cth_fc) - 1)
                                / (np.exp(cth_s - cth_fc) - 1)
                            )
                        )
                        if (thnew[ii] - dthdt) < th_fc_Adj_init[ii]:
                            dthdt = thnew[ii] - th_fc_Adj_init[ii]

                    # Update water content in compartment ii
                    thnew[ii] = cth_s - dthdt

                    # Update drainage from compartment ii
                    draincomp = dthdt * 1000 * cdz
                    # Update maximum drainage
                    drainmax = dthdt * 1000 * prethick

                    # Update excess drainage
                    if drainmax > excess:
                        drainmax = excess

                    excess = excess - drainmax
                    # Update drainsum and restrict to saturated hydraulic
                    # conductivity of soil layer
                    drainsum = draincomp + drainmax
                    if drainsum > cKsat:
                        excess = excess + drainsum - cKsat
                        drainsum = cKsat

        # Store output flux from compartment ii
        FluxOut[ii] = drainsum

        # Redistribute excess in compartment above
        if excess > 0:
            precomp = ii + 1
            while (excess > 0) and (precomp != 0):
                # Update compartment counter
                precomp = precomp - 1
                # Update layer counter
                # precompdf = Soil.Profile.Comp[precomp]

                # Update flux from compartment
                if precomp < ii:
                    FluxOut[precomp] = FluxOut[precomp] - excess

                # Increase water content to store excess
                thnew[precomp] = thnew[precomp] + (excess / (1000 * prof.dz[precomp]))

                # Limit water content to saturation and adjust excess counter
                if thnew[precomp] > prof.th_s[precomp]:
                    excess = (
                        (thnew[precomp] - prof.th_s[precomp]) * 1000 * prof.dz[precomp]
                    )
                    thnew[precomp] = prof.th_s[precomp]
                else:
                    excess = 0

    ## Update conditions and outputs ##
    # Total deep percolation (mm)
    DeepPerc = drainsum
    # Water contents
    # NewCond.th = thnew

    return thnew, DeepPerc, FluxOut

aquacrop.solution.evap_layer_water_content

evap_layer_water_content(InitCond_th, InitCond_EvapZ, prof)

Function to get water contents in the evaporation layer

Reference Manual: evaporation equations (pg. 73-81)

Arguments:

InitCond_th (numpy.array): Initial water content

InitCond_EvapZ (float): evaporation depth

prof (SoilProfileNT): Soil object containing soil paramaters

Returns:

Wevap_Sat (float): Water storage in evaporation layer at saturation (mm)

Wevap_Fc (float): Water storage in evaporation layer at field capacity (mm)

Wevap_Wp (float): Water storage in evaporation layer at permanent wilting point (mm)

Wevap_Dry (float): Water storage in evaporation layer at air dry (mm)

Wevap_Act (float): Actual water storage in evaporation layer (mm)
Source code in aquacrop/solution/evap_layer_water_content.py 
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
@njit
@cc.export("evap_layer_water_content", (f8[:],f8,SoilProfileNT_typ_sig))
def evap_layer_water_content(
    InitCond_th: "ndarray",
    InitCond_EvapZ: float,
    prof: "SoilProfileNT",
) -> Tuple[float, float, float, float, float]:
    """
    Function to get water contents in the evaporation layer

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=82" target="_blank">Reference Manual: evaporation equations</a> (pg. 73-81)


    Arguments:

        InitCond_th (numpy.array): Initial water content

        InitCond_EvapZ (float): evaporation depth

        prof (SoilProfileNT): Soil object containing soil paramaters


    Returns:


        Wevap_Sat (float): Water storage in evaporation layer at saturation (mm)

        Wevap_Fc (float): Water storage in evaporation layer at field capacity (mm)

        Wevap_Wp (float): Water storage in evaporation layer at permanent wilting point (mm)

        Wevap_Dry (float): Water storage in evaporation layer at air dry (mm)

        Wevap_Act (float): Actual water storage in evaporation layer (mm)



    """

    # Find soil compartments covered by evaporation layer
    comp_sto = np.sum(prof.dzsum < InitCond_EvapZ) + 1

    Wevap_Sat = 0
    Wevap_Fc = 0
    Wevap_Wp = 0
    Wevap_Dry = 0
    Wevap_Act = 0

    for ii in range(int(comp_sto)):

        # Determine fraction of soil compartment covered by evaporation layer
        if prof.dzsum[ii] > InitCond_EvapZ:
            factor = 1 - ((prof.dzsum[ii] - InitCond_EvapZ) / prof.dz[ii])
        else:
            factor = 1

        # Actual water storage in evaporation layer (mm)
        Wevap_Act += factor * 1000 * InitCond_th[ii] * prof.dz[ii]
        # Water storage in evaporation layer at saturation (mm)
        Wevap_Sat += factor * 1000 * prof.th_s[ii] * prof.dz[ii]
        # Water storage in evaporation layer at field capacity (mm)
        Wevap_Fc += factor * 1000 * prof.th_fc[ii] * prof.dz[ii]
        # Water storage in evaporation layer at permanent wilting point (mm)
        Wevap_Wp += factor * 1000 * prof.th_wp[ii] * prof.dz[ii]
        # Water storage in evaporation layer at air dry (mm)
        Wevap_Dry += factor * 1000 * prof.th_dry[ii] * prof.dz[ii]

    if Wevap_Act < 0:
        Wevap_Act = 0

    return Wevap_Sat, Wevap_Fc, Wevap_Wp, Wevap_Dry, Wevap_Act

aquacrop.solution.germination

germination(InitCond, Soil_zGerm, prof, Crop_GermThr, Crop_PlantMethod, gdd, growing_season)

Function to check if crop has germinated

Reference Manual: germination condition (pg. 23)

Arguments:

InitCond (InitialCondition): InitCond object containing model paramaters

Soil_zGerm (float): Soil depth affecting germination

prof (SoilProfileNT): Soil object containing soil paramaters

Crop_GermThr (float): Crop germination threshold

Crop_PlantMethod (bool): sown as seedling True or False

gdd (float): Number of Growing Degree Days on current day

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters
Source code in aquacrop/solution/germination.py 
 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
def germination(
    InitCond: "InitialCondition",
    Soil_zGerm: float,
    prof: "SoilProfileNT",
    Crop_GermThr: float,
    Crop_PlantMethod: bool,
    gdd: float,
    growing_season: bool,
    ) -> "InitialCondition":
    """
    Function to check if crop has germinated


    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=32" target="_blank">Reference Manual: germination condition</a> (pg. 23)


    Arguments:


        InitCond (InitialCondition): InitCond object containing model paramaters

        Soil_zGerm (float): Soil depth affecting germination

        prof (SoilProfileNT): Soil object containing soil paramaters

        Crop_GermThr (float): Crop germination threshold

        Crop_PlantMethod (bool): sown as seedling True or False

        gdd (float): Number of Growing Degree Days on current day

        growing_season (bool): is growing season (True or Flase)


    Returns:


        NewCond (InitialCondition): InitCond object containing updated model paramaters







    """

    ## Store initial conditions in new structure for updating ##
    NewCond = InitCond

    ## Check for germination (if in growing season) ##
    if growing_season == True:

        if (NewCond.germination == False):
            # Find compartments covered by top soil layer affecting germination
            comp_sto = np.argwhere(prof.dzsum >= Soil_zGerm).flatten()[0]
            # Calculate water content in top soil layer
            Wr = 0
            WrFC = 0
            WrWP = 0
            for ii in range(comp_sto + 1):
                # Get soil layer
                # Determine fraction of compartment covered by top soil layer
                if prof.dzsum[ii] > Soil_zGerm:
                    factor = 1 - ((prof.dzsum[ii] - Soil_zGerm) / prof.dz[ii])
                else:
                    factor = 1

                # Increment actual water storage (mm)
                Wr = Wr + round(factor * 1000 * InitCond.th[ii] * prof.dz[ii], 3)
                # Increment water storage at field capacity (mm)
                WrFC = WrFC + round(factor * 1000 * prof.th_fc[ii] * prof.dz[ii], 3)
                # Increment water storage at permanent wilting point (mm)
                WrWP = WrWP + round(factor * 1000 * prof.th_wp[ii] * prof.dz[ii], 3)

            # Limit actual water storage to not be less than zero
            if Wr < 0:
                Wr = 0

            # Calculate proportional water content
            WcProp = 1 - ((WrFC - Wr) / (WrFC - WrWP))

            # Check if water content is above germination threshold
            if (WcProp >= Crop_GermThr):
                # Crop has germinated
                NewCond.germination = True
                # If crop sown as seedling, turn on seedling protection
                if Crop_PlantMethod == True:
                    NewCond.protected_seed = True
                else:
                    # Crop is transplanted so no protection
                    NewCond.protected_seed = False

            # Increment delayed growth time counters if germination is yet to
            # occur, and also set seed protection to False if yet to germinate
            else:
                NewCond.delayed_cds = InitCond.delayed_cds + 1
                NewCond.delayed_gdds = InitCond.delayed_gdds + gdd
                NewCond.protected_seed = False

    else:
        # Not in growing season so no germination calculation is performed.
        NewCond.germination = False
        NewCond.protected_seed = False
        NewCond.delayed_cds = 0
        NewCond.delayed_gdds = 0

    return NewCond

aquacrop.solution.groundwater_inflow

groundwater_inflow(prof, NewCond)

Function to calculate capillary rise in the presence of a shallow groundwater table

Reference Manual: capillary rise calculations (pg. 52-61)

Arguments:

prof (SoilProfileNT): Soil profile parameters

NewCond (InitialCondition): model parameters

Returns:

NewCond (InitialCondition): InitCond object containing updated model parameters

GwIn (float): Groundwater inflow
Source code in aquacrop/solution/groundwater_inflow.py 
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
def groundwater_inflow(
    prof: "SoilProfileNT",
    NewCond: "InitialCondition"
    ) -> Tuple["InitialCondition",float]:
    """
    Function to calculate capillary rise in the presence of a shallow groundwater table

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=61" target="_blank">Reference Manual: capillary rise calculations</a> (pg. 52-61)


    Arguments:

        prof (SoilProfileNT): Soil profile parameters

        NewCond (InitialCondition): model parameters


    Returns:


        NewCond (InitialCondition): InitCond object containing updated model parameters

        GwIn (float): Groundwater inflow


    """

    ## Store initial conditions for updating ##
    GwIn = 0

    ## Perform calculations ##
    if NewCond.wt_in_soil == True:
        # Water table in soil profile. Calculate horizontal inflow.
        # Get groundwater table elevation on current day
        # print(f'Setting z_gw in groundwater_inflow.py as: {NewCond.z_gw}')
        z_gw = NewCond.z_gw

        # Find compartment mid-points
        zMid = prof.zMid

        idx = np.argwhere(zMid >= z_gw).flatten()[0]

        for ii in range(idx, len(prof.Comp)):
            # Get soil layer
            if NewCond.th[ii] < prof.th_s[ii]:
                # Update water content
                dth = prof.th_s[ii] - NewCond.th[ii]
                NewCond.th[ii] = prof.th_s[ii]
                # Update groundwater inflow
                GwIn = GwIn + (dth * 1000 * prof.dz[ii])

    return NewCond, GwIn

aquacrop.solution.growing_degree_day

growing_degree_day(GDDmethod, Tupp, Tbase, temp_max, temp_min)

Function to calculate number of growing degree days on current day

Reference manual: growing degree day calculations (pg. 19-20)

Arguments:

GDDmethod (int): gdd calculation method

Tupp (float): Upper temperature (degC) above which crop development no longer increases

Tbase (float): Base temperature (degC) below which growth does not progress

temp_max (float): Maximum tempature on current day (celcius)

temp_min (float): Minimum tempature on current day (celcius)

Returns:

gdd (float): Growing degree days for current day
Source code in aquacrop/solution/growing_degree_day.py 
 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
@cc.export("growing_degree_day", "f8(i4,f8,f8,f8,f8)")
def growing_degree_day(
    GDDmethod: int,
    Tupp: float,
    Tbase: float,
    temp_max: float,
    temp_min: float,
    ):
    """
    Function to calculate number of growing degree days on current day

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=28" target="_blank">Reference manual: growing degree day calculations</a> (pg. 19-20)


    Arguments:

        GDDmethod (int): gdd calculation method

        Tupp (float): Upper temperature (degC) above which crop development no longer increases

        Tbase (float): Base temperature (degC) below which growth does not progress

        temp_max (float): Maximum tempature on current day (celcius)

        temp_min (float): Minimum tempature on current day (celcius)


    Returns:

        gdd (float): Growing degree days for current day



    """

    ## Calculate GDDs ##
    if GDDmethod == 1:
        # method 1
        Tmean = (temp_max + temp_min) / 2
        Tmean = min(Tmean, Tupp)
        Tmean = max(Tmean, Tbase)
        gdd = Tmean - Tbase
    elif GDDmethod == 2:
        # method 2
        temp_max = min(temp_max, Tupp)
        temp_max = max(temp_max, Tbase)

        temp_min = min(temp_min, Tupp)
        temp_min = max(temp_min, Tbase)

        Tmean = (temp_max + temp_min) / 2
        gdd = Tmean - Tbase
    elif GDDmethod == 3:
        # method 3
        temp_max = min(temp_max, Tupp)
        temp_max = max(temp_max, Tbase)

        temp_min = min(temp_min, Tupp)
        Tmean = (temp_max + temp_min) / 2
        Tmean = max(Tmean, Tbase)
        gdd = Tmean - Tbase

    return gdd

aquacrop.solution.growth_stage

growth_stage(Crop, InitCond, growing_season)

Function to determine current growth stage of crop

(used only for irrigation soil moisture thresholds)

Arguments:

Crop (Crop): Crop object containing Crop paramaters

InitCond (InitialCondition): InitCond object containing model paramaters

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters
Source code in aquacrop/solution/growth_stage.py 
 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
def growth_stage(
    Crop: "Crop",
    InitCond: "InitialCondition",
    growing_season: bool,
    ) -> "InitialCondition":
    """
    Function to determine current growth stage of crop

    (used only for irrigation soil moisture thresholds)

    Arguments:

        Crop (Crop): Crop object containing Crop paramaters

        InitCond (InitialCondition): InitCond object containing model paramaters

        growing_season (bool): is growing season (True or Flase)


    Returns:

        NewCond (InitialCondition): InitCond object containing updated model paramaters



    """

    ## Store initial conditions in new structure for updating ##
    NewCond = InitCond

    ## Get growth stage (if in growing season) ##
    if growing_season == True:
        # Adjust time for any delayed growth
        if Crop.CalendarType == 1:
            tAdj = NewCond.dap - NewCond.delayed_cds
        elif Crop.CalendarType == 2:
            tAdj = NewCond.gdd_cum - NewCond.delayed_gdds

        # Update growth stage
        if tAdj <= Crop.Canopy10Pct:
            NewCond.growth_stage = 1
        elif tAdj <= Crop.MaxCanopy:
            NewCond.growth_stage = 2
        elif tAdj <= Crop.Senescence:
            NewCond.growth_stage = 3
        elif tAdj > Crop.Senescence:
            NewCond.growth_stage = 4

    else:
        # Not in growing season so growth stage is set to dummy value
        NewCond.growth_stage = 0

    return NewCond

aquacrop.solution.harvest_index

harvest_index(prof, Soil_zTop, Crop, InitCond, et0, temp_max, temp_min, growing_season)

Function to simulate build up of harvest index

Reference Manual: harvest index calculations (pg. 110-126)

Arguments:

prof (SoilProfileNT): Soil profile paramaters

Soil_zTop (float): topsoil depth

Crop (CropStructNT): Crop parameters

InitCond (InitialCondition): InitCond object containing model paramaters

et0 (float): reference evapotranspiration on current day

temp_max (float): maximum tempature on current day (celcius)

temp_min (float): minimum tempature on current day (celcius)

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters
Source code in aquacrop/solution/harvest_index.py 
 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
def harvest_index(
    prof: "SoilProfileNT",
    Soil_zTop: float,
    Crop: "CropStructNT",
    InitCond: "InitialCondition",
    et0: float,
    temp_max: float,
    temp_min: float,
    growing_season: bool,
    ) -> "InitialCondition":

    """
    Function to simulate build up of harvest index


     <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual: harvest index calculations</a> (pg. 110-126)

    Arguments:


        prof (SoilProfileNT): Soil profile paramaters

        Soil_zTop (float): topsoil depth

        Crop (CropStructNT): Crop parameters

        InitCond (InitialCondition): InitCond object containing model paramaters

        et0 (float): reference evapotranspiration on current day

        temp_max (float): maximum tempature on current day (celcius)

        temp_min (float): minimum tempature on current day (celcius)

        growing_season (bool): is growing season (True or Flase)


    Returns:


        NewCond (InitialCondition): InitCond object containing updated model paramaters



    """

    ## Store initial conditions for updating ##
    NewCond = InitCond

    InitCond_HI = InitCond.harvest_index
    InitCond_HIadj = InitCond.harvest_index_adj
    InitCond_PreAdj = InitCond.pre_adj

    ## Calculate harvest index build up (if in growing season) ##
    if growing_season == True:
        # Calculate root zone water content

        taw = TAW()
        water_root_depletion = Dr()
        # thRZ = RootZoneWater()
        _, water_root_depletion.Zt, water_root_depletion.Rz, taw.Zt, taw.Rz, _,_,_,_,_,_, = root_zone_water(
            prof,
            float(NewCond.z_root),
            NewCond.th,
            Soil_zTop,
            float(Crop.Zmin),
            Crop.Aer,
        )

        # _,water_root_depletion,taw,_ = root_zone_water(Soil_Profile,float(NewCond.z_root),NewCond.th,Soil_zTop,float(Crop.Zmin),Crop.Aer)
        # Check whether to use root zone or top soil depletions for calculating
        # water stress
        if (water_root_depletion.Rz / taw.Rz) <= (water_root_depletion.Zt / taw.Zt):
            # Root zone is wetter than top soil, so use root zone value
            water_root_depletion = water_root_depletion.Rz
            taw = taw.Rz
        else:
            # Top soil is wetter than root zone, so use top soil values
            water_root_depletion = water_root_depletion.Zt
            taw = taw.Zt

        # Calculate water stress
        beta = True
        # Ksw = water_stress(Crop, NewCond, water_root_depletion, taw, et0, beta)
        # Ksw = Ksw()
        Ksw_Exp, Ksw_Sto, Ksw_Sen, Ksw_Pol, Ksw_StoLin = water_stress(
            Crop.p_up,
            Crop.p_lo,
            Crop.ETadj,
            Crop.beta,
            Crop.fshape_w,
            NewCond.t_early_sen,
            water_root_depletion,
            taw,
            et0,
            beta,
        )
        Ksw = KswNT(exp=Ksw_Exp, sto=Ksw_Sto, sen=Ksw_Sen, pol=Ksw_Pol, sto_lin=Ksw_StoLin )
        # Calculate temperature stress
        (Kst_PolH,Kst_PolC) = temperature_stress(Crop, temp_max, temp_min)
        Kst = KstNT(PolH=Kst_PolH,PolC=Kst_PolC)
        # Get reference harvest index on current day
        HIi = NewCond.hi_ref

        # Get time for harvest index build-up
        HIt = NewCond.dap - NewCond.delayed_cds - Crop.HIstartCD - 1

        # Calculate harvest index
        if (NewCond.yield_form == True) and (HIt >= 0):
            # print(NewCond.dap)
            # Root/tuber or fruit/grain crops
            if (Crop.CropType == 2) or (Crop.CropType == 3):
                # Detemine adjustment for water stress before anthesis
                if InitCond_PreAdj == False:
                    InitCond.pre_adj = True
                    NewCond.f_pre = HIadj_pre_anthesis(NewCond.biomass,
                                                NewCond.biomass_ns,
                                                NewCond.canopy_cover,
                                                Crop.dHI_pre)

                # Determine adjustment for crop pollination failure
                if Crop.CropType == 3:  # Adjustment only for fruit/grain crops
                    if (HIt > 0) and (HIt <= Crop.FloweringCD):

                        NewCond.f_pol = HIadj_pollination(
                            NewCond.canopy_cover,
                            NewCond.f_pol,
                            Crop.FloweringCD,
                            Crop.CCmin,
                            Crop.exc,
                            Ksw,
                            Kst,
                            HIt,
                        )

                    HImax = NewCond.f_pol * Crop.HI0
                else:
                    # No pollination adjustment for root/tuber crops
                    HImax = Crop.HI0

                # Determine adjustments for post-anthesis water stress
                if HIt > 0:
                    (NewCond.s_cor1,
                    NewCond.s_cor2,
                    NewCond.fpost_upp,
                    NewCond.fpost_dwn,
                    NewCond.f_post) = HIadj_post_anthesis(NewCond.delayed_cds,
                                                        NewCond.s_cor1,
                                                        NewCond.s_cor2,
                                                        NewCond.dap,
                                                        NewCond.f_pre,
                                                        NewCond.canopy_cover,
                                                        NewCond.fpost_upp,
                                                        NewCond.fpost_dwn,
                                                        Crop, 
                                                        Ksw)

                # Limit harvest_index to maximum allowable increase due to pre- and
                # post-anthesis water stress combinations
                HImult = NewCond.f_pre * NewCond.f_post
                if HImult > 1 + (Crop.dHI0 / 100):
                    HImult = 1 + (Crop.dHI0 / 100)

                # Determine harvest index on current day, adjusted for stress
                # effects
                if HImax >= HIi:
                    harvest_index_adj = HImult * HIi
                else:
                    harvest_index_adj = HImult * HImax

            elif Crop.CropType == 1:
                # Leafy vegetable crops - no adjustment, harvest index equal to
                # reference value for current day
                harvest_index_adj = HIi

        else:

            # No build-up of harvest index if outside yield_ formation period
            HIi = InitCond_HI
            harvest_index_adj = InitCond_HIadj

        # Store final values for current time step
        NewCond.harvest_index = HIi
        NewCond.harvest_index_adj = harvest_index_adj

    else:
        # No harvestable crop outside of a growing season
        NewCond.harvest_index = 0
        NewCond.harvest_index_adj = 0

    # print([NewCond.dap , Crop.YldFormCD])
    return NewCond

aquacrop.solution.HIadj_pollination

HIadj_pollination(NewCond_CC, NewCond_Fpol, Crop_FloweringCD, Crop_CCmin, Crop_exc, Ksw, Kst, HIt)

Function to calculate adjustment to harvest index for failure of pollination due to water or temperature stress

Reference Manual: harvest index calculations (pg. 110-126)

Arguments:

NewCond_CC (float): InitCond object containing model paramaters

NewCond_Fpol (float): InitCond object containing model paramaters

Crop_FloweringCD (float): Length of flowering stage

Crop_CCmin (float): minimum canopy cover

Crop_exc (float):

Ksw (KswNT): Ksw object containing water stress paramaters

Kst (KstNT): Kst object containing tempature stress paramaters

HIt (float): time for harvest index build-up (calander days)

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters
Source code in aquacrop/solution/HIadj_pollination.py 
 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
@cc.export("HIadj_pollination", (f8,f8,f8,f8,f8,KswNT_type_sig,KstNT_type_sig,f8))
def HIadj_pollination(
    NewCond_CC: float,
    NewCond_Fpol: float,
    Crop_FloweringCD: float, 
    Crop_CCmin: float, 
    Crop_exc: float, 
    Ksw: "KswNT", 
    Kst: "KstNT", 
    HIt: float,
) -> float:
    """
    Function to calculate adjustment to harvest index for failure of
    pollination due to water or temperature stress

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual: harvest index calculations</a> (pg. 110-126)


    Arguments:


        NewCond_CC (float): InitCond object containing model paramaters

        NewCond_Fpol (float): InitCond object containing model paramaters

        Crop_FloweringCD (float): Length of flowering stage

        Crop_CCmin (float): minimum canopy cover

        Crop_exc (float): 

        Ksw (KswNT): Ksw object containing water stress paramaters

        Kst (KstNT): Kst object containing tempature stress paramaters

        HIt (float): time for harvest index build-up (calander days)


    Returns:


        NewCond (InitialCondition): InitCond object containing updated model paramaters



    """

    ## Caclulate harvest index adjustment for pollination ##
    # Get fractional flowering
    if HIt == 0:
        # No flowering yet
        FracFlow = 0
    elif HIt > 0:
        # Fractional flowering on previous day
        t1 = HIt - 1
        if t1 == 0:
            F1 = 0
        else:
            t1Pct = 100 * (t1 / Crop_FloweringCD)
            if t1Pct > 100:
                t1Pct = 100

            F1 = 0.00558 * np.exp(0.63 * np.log(t1Pct)) - (0.000969 * t1Pct) - 0.00383

        if F1 < 0:
            F1 = 0

        # Fractional flowering on current day
        t2 = HIt
        if t2 == 0:
            F2 = 0
        else:
            t2Pct = 100 * (t2 / Crop_FloweringCD)
            if t2Pct > 100:
                t2Pct = 100

            F2 = 0.00558 * np.exp(0.63 * np.log(t2Pct)) - (0.000969 * t2Pct) - 0.00383

        if F2 < 0:
            F2 = 0

        # Weight values
        if abs(F1 - F2) < 0.0000001:
            F = 0
        else:
            F = 100 * ((F1 + F2) / 2) / Crop_FloweringCD

        FracFlow = F

    # Calculate pollination adjustment for current day
    if NewCond_CC < Crop_CCmin:
        # No pollination can occur as canopy cover is smaller than minimum
        # threshold
        dFpol = 0
    else:
        Ks = min([Ksw.pol, Kst.PolC, Kst.PolH])
        dFpol = Ks * FracFlow * (1 + (Crop_exc / 100))

    # Calculate pollination adjustment to date
    NewCond_Fpol = NewCond_Fpol + dFpol
    if NewCond_Fpol > 1:
        # Crop has fully pollinated
        NewCond_Fpol = 1

    return NewCond_Fpol

aquacrop.solution.HIadj_post_anthesis

HIadj_post_anthesis(NewCond_DelayedCDs, NewCond_sCor1, NewCond_sCor2, NewCond_DAP, NewCond_Fpre, NewCond_CC, NewCond_fpost_upp, NewCond_fpost_dwn, Crop, Ksw)

Function to calculate adjustment to harvest index for post-anthesis water stress

Reference Manual: harvest index calculations (pg. 110-126)

Arguments:

NewCond_DelayedCDs (int): delayed calendar days

NewCond_sCor1 (float): canopy exapnsion

NewCond_sCor2 (float): stomatal closure

NewCond_DAP (int): days since planting

NewCond_Fpre (float): delayed calendar days

NewCond_CC (float): current canopy cover

NewCond_fpost_upp (float): delayed calendar days

NewCond_fpost_dwn (float): delayed calendar days

Crop (CropStructNT): Crop paramaters

Ksw (KswNT): water stress paramaters

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters
Source code in aquacrop/solution/HIadj_post_anthesis.py 
 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
@cc.export("HIadj_post_anthesis", (i8,f8,f8,i8,f8,f8,f8,f8,CropStructNT_type_sig,KswNT_type_sig,))
def HIadj_post_anthesis(
    NewCond_DelayedCDs: int,
    NewCond_sCor1: float,
    NewCond_sCor2: float,
    NewCond_DAP: int,
    NewCond_Fpre: float,
    NewCond_CC: float,
    NewCond_fpost_upp: float,
    NewCond_fpost_dwn: float,
    Crop: "CropStructNT", 
    Ksw: "KswNT",
    ) -> Tuple[float, float, float, float, float]:
    """
    Function to calculate adjustment to harvest index for post-anthesis water
    stress

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual: harvest index calculations</a> (pg. 110-126)


    Arguments:


        NewCond_DelayedCDs (int): delayed calendar days

        NewCond_sCor1 (float): canopy exapnsion

        NewCond_sCor2 (float): stomatal closure

        NewCond_DAP (int): days since planting

        NewCond_Fpre (float): delayed calendar days

        NewCond_CC (float): current canopy cover

        NewCond_fpost_upp (float): delayed calendar days

        NewCond_fpost_dwn (float): delayed calendar days

        Crop (CropStructNT): Crop paramaters

        Ksw (KswNT): water stress paramaters

    Returns:


        NewCond (InitialCondition): InitCond object containing updated model paramaters


    """

    ## Store initial conditions in a structure for updating ##
    # NewCond = InitCond

    InitCond_DelayedCDs = NewCond_DelayedCDs*1
    InitCond_sCor1 = NewCond_sCor1*1
    InitCond_sCor2 = NewCond_sCor2*1

    ## Calculate harvest index adjustment ##
    # 1. Adjustment for leaf expansion
    tmax1 = Crop.CanopyDevEndCD - Crop.HIstartCD
    dap = NewCond_DAP - InitCond_DelayedCDs
    if (
        (dap <= (Crop.CanopyDevEndCD + 1))
        and (tmax1 > 0)
        and (NewCond_Fpre > 0.99)
        and (NewCond_CC > 0.001)
        and (Crop.a_HI > 0)
    ):
        dCor = 1 + (1 - Ksw.exp) / Crop.a_HI
        NewCond_sCor1 = InitCond_sCor1 + (dCor / tmax1)
        DayCor = dap - 1 - Crop.HIstartCD
        NewCond_fpost_upp = (tmax1 / DayCor) * NewCond_sCor1

    # 2. Adjustment for stomatal closure
    tmax2 = Crop.YldFormCD
    dap = NewCond_DAP - InitCond_DelayedCDs
    if (
        (dap <= (Crop.HIendCD + 1))
        and (tmax2 > 0)
        and (NewCond_Fpre > 0.99)
        and (NewCond_CC > 0.001)
        and (Crop.b_HI > 0)
    ):
        # print(Ksw.sto)
        dCor = np.power(Ksw.sto, 0.1) * (1 - (1 - Ksw.sto) / Crop.b_HI)
        NewCond_sCor2 = InitCond_sCor2 + (dCor / tmax2)
        DayCor = dap - 1 - Crop.HIstartCD
        NewCond_fpost_dwn = (tmax2 / DayCor) * NewCond_sCor2

    # Determine total multiplier
    if (tmax1 == 0) and (tmax2 == 0):
        NewCond_Fpost = 1
    else:
        if tmax2 == 0:
            NewCond_Fpost = NewCond_fpost_upp
        else:
            if tmax1 == 0:
                NewCond_Fpost = NewCond_fpost_dwn
            elif tmax1 <= tmax2:
                NewCond_Fpost = NewCond_fpost_dwn * (
                    ((tmax1 * NewCond_fpost_upp) + (tmax2 - tmax1)) / tmax2
                )
            else:
                NewCond_Fpost = NewCond_fpost_upp * (
                    ((tmax2 * NewCond_fpost_dwn) + (tmax1 - tmax2)) / tmax1
                )

    return (
            NewCond_sCor1,
            NewCond_sCor2,
            NewCond_fpost_upp,
            NewCond_fpost_dwn,
            NewCond_Fpost)

aquacrop.solution.HIadj_pre_anthesis

HIadj_pre_anthesis(NewCond_B, NewCond_B_NS, NewCond_CC, Crop_dHI_pre)

Function to calculate adjustment to harvest index for pre-anthesis water stress

Reference Manual: harvest index calculations (pg. 110-126)

Arguments:

NewCond_B (float): biomass growth

NewCond_B_NS (float): biomass growth (no stress)

NewCond_CC (float): canopy cover

Crop_dHI_pre (float): Crop_dHI_pre

Returns:

NewCond_Fpre (float): adjustment to harvest index for pre-anthesis water stress
Source code in aquacrop/solution/HIadj_pre_anthesis.py 
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
@cc.export("HIadj_pre_anthesis", (f8,f8,f8,f8))
def HIadj_pre_anthesis(
    NewCond_B: float,
    NewCond_B_NS: float,
    NewCond_CC: float,
    Crop_dHI_pre: float
    ) -> float:
    """
    Function to calculate adjustment to harvest index for pre-anthesis water
    stress

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual: harvest index calculations</a> (pg. 110-126)


    Arguments:

        NewCond_B (float): biomass growth

        NewCond_B_NS (float): biomass growth (no stress)

        NewCond_CC (float): canopy cover

        Crop_dHI_pre (float): Crop_dHI_pre


    Returns:


        NewCond_Fpre (float): adjustment to harvest index for pre-anthesis water stress


    """

    ## Store initial conditions in structure for updating ##
    # NewCond = InitCond

    # check that there is an adjustment to be made
    if Crop_dHI_pre > 0:
        ## Calculate adjustment ##
        # Get parameters
        Br = NewCond_B / NewCond_B_NS
        Br_range = np.log(Crop_dHI_pre) / 5.62
        Br_upp = 1
        Br_low = 1 - Br_range
        Br_top = Br_upp - (Br_range / 3)

        # Get biomass ratios
        ratio_low = (Br - Br_low) / (Br_top - Br_low)
        ratio_upp = (Br - Br_top) / (Br_upp - Br_top)

        # Calculate adjustment factor
        if (Br >= Br_low) and (Br < Br_top):
            NewCond_Fpre = 1 + (
                ((1 + np.sin((1.5 - ratio_low) * np.pi)) / 2) * (Crop_dHI_pre / 100)
            )
        elif (Br > Br_top) and (Br <= Br_upp):
            NewCond_Fpre = 1 + (
                ((1 + np.sin((0.5 + ratio_upp) * np.pi)) / 2) * (Crop_dHI_pre / 100)
            )
        else:
            NewCond_Fpre = 1
    else:
        NewCond_Fpre = 1

    if NewCond_CC <= 0.01:
        # No green canopy cover left at start of flowering so no harvestable
        # crop will develop
        NewCond_Fpre = 0

    return NewCond_Fpre

aquacrop.solution.HIref_current_day

HIref_current_day(NewCond_HIref, NewCond_HIfinal, NewCond_DAP, NewCond_DelayedCDs, NewCond_YieldForm, NewCond_PctLagPhase, NewCond_CC, NewCond_CC_prev, NewCond_CCxW, Crop, growing_season)

Function to calculate reference (no adjustment for stress effects) harvest index on current day

Reference Manual: harvest index calculations (pg. 110-126)

Arguments:

NewCond_HIref (float): reference harvest index

NewCond_HIfinal (float): final harvest index to track effects of early canopy decline

NewCond_DAP (int): days after planting

NewCond_DelayedCDs (int): delayed calendar days

NewCond_YieldForm (bool): yield formation stage

NewCond_PctLagPhase (float): percent through eraly development phase

NewCond_CC (float): canopy cover on current day

NewCond_CCxW (float): max canopy cover during season accounting for any early senescence

Crop (CropStructNT): Crop paramaters

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond (NewCond_HIref): reference harvest index

NewCond (NewCond_YieldForm): yield formation stage

NewCond (NewCond_PctLagPhase): percent through early development phase

NewCond (NewCond_HIfinal): final harvest index to track effects of early canopy decline
Source code in aquacrop/solution/HIref_current_day.py 
 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
@cc.export("HIref_current_day", (f8,f8,i8,i8,b1,f8,f8,f8,f8,CropStructNT_type_sig,b1))
def HIref_current_day(
    NewCond_HIref: float,
    NewCond_HIfinal: float,
    NewCond_DAP: int,
    NewCond_DelayedCDs: int,
    NewCond_YieldForm: bool,
    NewCond_PctLagPhase: float,
    NewCond_CC: float,
    NewCond_CC_prev: float,
    NewCond_CCxW: float,
    Crop: "CropStructNT",
    growing_season: bool,
    ) -> Tuple[float, bool, float]: #, float
    """
    Function to calculate reference (no adjustment for stress effects)
    harvest index on current day

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=119" target="_blank">Reference Manual: harvest index calculations</a> (pg. 110-126)



    Arguments:


        NewCond_HIref (float): reference harvest index

        NewCond_HIfinal (float): final harvest index to track effects of early canopy decline

        NewCond_DAP (int): days after planting

        NewCond_DelayedCDs (int): delayed calendar days

        NewCond_YieldForm (bool): yield formation stage

        NewCond_PctLagPhase (float): percent through eraly development phase

        NewCond_CC (float): canopy cover on current day

        NewCond_CCxW (float): max canopy cover during season accounting for any early senescence

        Crop (CropStructNT): Crop paramaters

        growing_season (bool): is growing season (True or Flase)


    Returns:


        NewCond (NewCond_HIref): reference harvest index

        NewCond (NewCond_YieldForm): yield formation stage

        NewCond (NewCond_PctLagPhase): percent through early development phase

        NewCond (NewCond_HIfinal): final harvest index to track effects of early canopy decline


    """

    ## Store initial conditions for updating ##
    # NewCond = InitCond

    InitCond_HIref = NewCond_HIref*1

    # NewCond.hi_ref = 0.

    ## Calculate reference harvest index (if in growing season) ##
    if growing_season == True:
        # Check if in yield_ formation period
        tAdj = NewCond_DAP - NewCond_DelayedCDs
        if tAdj > Crop.HIstartCD:

            NewCond_YieldForm = True
        else:
            NewCond_YieldForm = False

        # Get time for harvest index calculation
        HIt = NewCond_DAP - NewCond_DelayedCDs - Crop.HIstartCD - 1

        if HIt <= 0:
            # Yet to reach time for harvest_index build-up
            NewCond_HIref = 0
            NewCond_PctLagPhase = 0
        else:
            # Check crop type
            if (Crop.CropType == 1) or (Crop.CropType == 2):
                # If crop type is leafy vegetable or root/tuber, then proceed with
                # logistic growth (i.e. no linear switch)
                NewCond_PctLagPhase = 100  # No lag phase
                # Calculate reference harvest index for current day
                NewCond_HIref = (Crop.HIini * Crop.HI0) / (
                    Crop.HIini + (Crop.HI0 - Crop.HIini) * np.exp(-Crop.HIGC * HIt))
                # Harvest index apprAOSP_hing maximum limit
                if NewCond_HIref >= (0.9799 * Crop.HI0):
                    NewCond_HIref = Crop.HI0

            elif Crop.CropType == 3:
                # If crop type is fruit/grain producing, check for linear switch
                if HIt < Crop.tLinSwitch:
                    # Not yet reached linear switch point, therefore proceed with
                    # logistic build-up
                    NewCond_PctLagPhase = 100 * (HIt / Crop.tLinSwitch)
                    # Calculate reference harvest index for current day
                    # (logistic build-up)
                    NewCond_HIref = (Crop.HIini * Crop.HI0) / (
                        Crop.HIini + (Crop.HI0 - Crop.HIini) * np.exp(-Crop.HIGC * HIt))
                else:
                    # Linear switch point has been reached
                    NewCond_PctLagPhase = 100
                    # Calculate reference harvest index for current day
                    # (logistic portion)
                    NewCond_HIref = (Crop.HIini * Crop.HI0) / (Crop.HIini
                        + (Crop.HI0 - Crop.HIini) * np.exp(-Crop.HIGC * Crop.tLinSwitch))
                    # Calculate reference harvest index for current day
                    # (total - logistic portion + linear portion)
                    NewCond_HIref = NewCond_HIref + (Crop.dHILinear * (HIt - Crop.tLinSwitch))

            # Limit hi_ref and round off computed value
            if NewCond_HIref > Crop.HI0:
                NewCond_HIref = Crop.HI0
            elif NewCond_HIref <= (Crop.HIini + 0.004):
                NewCond_HIref = 0
            elif (Crop.HI0 - NewCond_HIref) < 0.004:
                NewCond_HIref = Crop.HI0

            # Adjust hi_ref for inadequate photosynthesis
            if (NewCond_HIfinal == Crop.HI0) and (HIt <= Crop.YldFormCD) and (NewCond_CC <= 0.05) and (
                NewCond_CCxW > 0) and (NewCond_CC < NewCond_CCxW) and (Crop.CropType==2 or Crop.CropType==3):

                NewCond_HIfinal = NewCond_HIref

            if NewCond_HIref > NewCond_HIfinal:
                NewCond_HIref = NewCond_HIfinal

    else:
        # Reference harvest index is zero outside of growing season
        NewCond_HIref = 0

    return (NewCond_HIref,
            NewCond_YieldForm,
            NewCond_PctLagPhase,
            #NewCond_HIfinal,
            )

aquacrop.solution.infiltration

infiltration(prof, NewCond_SurfaceStorage, NewCond_th_fc_Adj, NewCond_th, Infl, Irr, IrrMngt_AppEff, FieldMngt_Bunds, FieldMngt_zBund, FluxOut, DeepPerc0, Runoff0, growing_season)

Function to infiltrate incoming water (rainfall and irrigation)

Reference Manual: drainage calculations (pg. 42-65)

Arguments:

prof (SoilProfile): Soil object containing soil paramaters

NewCond_SurfaceStorage (float): surface storage

NewCond_th_fc_Adj (ndarray): water content at field capacity

NewCond_th (ndarray): soil water content

Infl (float): Infiltration so far

Irr (float): Irrigation on current day

IrrMngt_AppEff (float`: irrigation application efficiency

FieldMngt (FieldMngtStruct): field management params

FluxOut (np.array): flux of water out of each compartment

DeepPerc0 (float): initial Deep Percolation

Runoff0 (float): initial Surface Runoff

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond_th (nunpy.darray): updated soil water content

NewCond_SurfaceStorage (float): updated surface storage

DeepPerc (float): Total Deep Percolation

RunoffTot (float): Total surface Runoff

Infl (float): Infiltration on current day

FluxOut (numpy.array): flux of water out of each compartment
Source code in aquacrop/solution/infiltration.py 
 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
@cc.export("infiltration", (SoilProfileNT_typ_sig,f8,f8[:],f8[:],f8,f8,f8,b1,f8,f8[:],f8,f8,b1))
def infiltration(
     prof: "SoilProfileNT",
     NewCond_SurfaceStorage: float, 
     NewCond_th_fc_Adj: "ndarray", 
     NewCond_th: "ndarray", 
     Infl: float, 
     Irr: float, 
     IrrMngt_AppEff: float, 
     FieldMngt_Bunds: bool,
     FieldMngt_zBund: float,
     FluxOut: "ndarray", 
     DeepPerc0: float, 
     Runoff0: float, 
     growing_season: bool,
) -> Tuple["ndarray", float, float, float, float, "ndarray"]:
    """
    Function to infiltrate incoming water (rainfall and irrigation)

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=51" target="_blank">Reference Manual: drainage calculations</a> (pg. 42-65)



    Arguments:



        prof (SoilProfile): Soil object containing soil paramaters

        NewCond_SurfaceStorage (float): surface storage

        NewCond_th_fc_Adj (ndarray): water content at field capacity

        NewCond_th (ndarray): soil water content

        Infl (float): Infiltration so far

        Irr (float): Irrigation on current day

        IrrMngt_AppEff (float`: irrigation application efficiency

        FieldMngt (FieldMngtStruct): field management params

        FluxOut (np.array): flux of water out of each compartment

        DeepPerc0 (float): initial Deep Percolation

        Runoff0 (float): initial Surface Runoff

        growing_season (bool): is growing season (True or Flase)


    Returns:

        NewCond_th (nunpy.darray): updated soil water content

        NewCond_SurfaceStorage (float): updated surface storage

        DeepPerc (float): Total Deep Percolation

        RunoffTot (float): Total surface Runoff

        Infl (float): Infiltration on current day

        FluxOut (numpy.array): flux of water out of each compartment




    """
    ## Store initial conditions in new structure for updating ##
    # NewCond = InitCond

    InitCond_SurfaceStorage = NewCond_SurfaceStorage*1
    InitCond_th_fc_Adj = NewCond_th_fc_Adj*1
    InitCond_th = NewCond_th*1

    thnew = NewCond_th*1.

    Soil_nComp = thnew.shape[0]

    Infl = max(Infl,0.)

    ## Update infiltration rate for irrigation ##
    # Note: irrigation amount adjusted for specified application efficiency
    if growing_season == True:
        Infl = Infl + (Irr * (IrrMngt_AppEff / 100))

    assert Infl >= 0

    ## Determine surface storage (if bunds are present) ##
    if FieldMngt_Bunds:
        # bunds on field
        if FieldMngt_zBund > 0.001:
            # Bund height too small to be considered
            InflTot = Infl + NewCond_SurfaceStorage
            if InflTot > 0:
                # Update surface storage and infiltration storage
                if InflTot > prof.Ksat[0]:
                    # Infiltration limited by saturated hydraulic conductivity
                    # of surface soil layer
                    ToStore = prof.Ksat[0]
                    # Additional water ponds on surface
                    NewCond_SurfaceStorage = InflTot - prof.Ksat[0]
                else:
                    # All water infiltrates
                    ToStore = InflTot
                    # Reset surface storage depth to zero
                    NewCond_SurfaceStorage = 0

                # Calculate additional runoff
                if NewCond_SurfaceStorage > FieldMngt_zBund:
                    # Water overtops bunds and runs off
                    RunoffIni = NewCond_SurfaceStorage - FieldMngt_zBund
                    # Surface storage equal to bund height
                    NewCond_SurfaceStorage = FieldMngt_zBund * 1
                else:
                    # No overtopping of bunds
                    RunoffIni = 0

            else:
                # No storage or runoff
                ToStore = 0
                RunoffIni = 0

    elif FieldMngt_Bunds == False:
        # No bunds on field
        if Infl > prof.Ksat[0]:
            # Infiltration limited by saturated hydraulic conductivity of top
            # soil layer
            ToStore = prof.Ksat[0]
            # Additional water runs off
            RunoffIni = Infl - prof.Ksat[0]
        else:
            # All water infiltrates
            ToStore = Infl
            RunoffIni = 0

        # Update surface storage
        NewCond_SurfaceStorage = 0
        # Add any water remaining behind bunds to surface runoff (needed for
        # days when bunds are removed to maintain water balance)
        RunoffIni = RunoffIni + InitCond_SurfaceStorage

    ## Initialise counters
    ii = -1
    Runoff = 0
    ## Infiltrate incoming water ##
    if ToStore > 0:
        while (ToStore > 0) and (ii < Soil_nComp - 1):
            # Update compartment counter
            ii = ii + 1
            # Get soil layer

            # Calculate saturated drainage ability
            dthdtS = prof.tau[ii] * (prof.th_s[ii] - prof.th_fc[ii])
            # Calculate drainage factor
            factor = prof.Ksat[ii] / (dthdtS * 1000 * prof.dz[ii])

            # Calculate drainage ability required
            dthdt0 = ToStore / (1000 * prof.dz[ii])

            # Check drainage ability
            if dthdt0 < dthdtS:
                # Calculate water content, thX, needed to meet drainage dthdt0
                if dthdt0 <= 0:
                    theta0 = InitCond_th_fc_Adj[ii]
                else:
                    A = 1 + (
                        (dthdt0 * (np.exp(prof.th_s[ii] - prof.th_fc[ii]) - 1))
                        / (prof.tau[ii] * (prof.th_s[ii] - prof.th_fc[ii]))
                    )

                    theta0 = prof.th_fc[ii] + np.log(A)

                # Limit thX to between saturation and field capacity
                if theta0 > prof.th_s[ii]:
                    theta0 = prof.th_s[ii]
                elif theta0 <= InitCond_th_fc_Adj[ii]:
                    theta0 = InitCond_th_fc_Adj[ii]
                    dthdt0 = 0

            else:
                # Limit water content and drainage to saturation
                theta0 = prof.th_s[ii]
                dthdt0 = dthdtS

            # Calculate maximum water flow through compartment ii
            drainmax = factor * dthdt0 * 1000 * prof.dz[ii]
            # Calculate total drainage from compartment ii
            drainage = drainmax + FluxOut[ii]
            # Limit drainage to saturated hydraulic conductivity
            if drainage > prof.Ksat[ii]:
                drainmax = prof.Ksat[ii] - FluxOut[ii]

            # Calculate difference between threshold and current water contents
            diff = theta0 - InitCond_th[ii]

            if diff > 0:
                # Increase water content of compartment ii
                thnew[ii] = thnew[ii] + (ToStore / (1000 * prof.dz[ii]))
                if thnew[ii] > theta0:
                    # Water remaining that can infiltrate to compartments below
                    ToStore = (thnew[ii] - theta0) * 1000 * prof.dz[ii]
                    thnew[ii] = theta0
                else:
                    # All infiltrating water has been stored
                    ToStore = 0

            # Update outflow from current compartment (drainage + infiltration
            # flows)
            FluxOut[ii] = FluxOut[ii] + ToStore

            # Calculate back-up of water into compartments above
            excess = ToStore - drainmax
            if excess < 0:
                excess = 0

            # Update water to store
            ToStore = ToStore - excess

            # Redistribute excess to compartments above
            if excess > 0:
                precomp = ii + 1
                while (excess > 0) and (precomp != 0):
                    # Keep storing in compartments above until soil surface is
                    # reached
                    # Update compartment counter
                    precomp = precomp - 1
                    # Update layer number

                    # Update outflow from compartment
                    FluxOut[precomp] = FluxOut[precomp] - excess
                    # Update water content
                    thnew[precomp] = thnew[precomp] + (excess / (prof.dz[precomp] * 1000))
                    # Limit water content to saturation
                    if thnew[precomp] > prof.th_s[precomp]:
                        # Update excess to store
                        excess = (thnew[precomp] - prof.th_s[precomp]) * 1000 * prof.dz[precomp]
                        # Set water content to saturation
                        thnew[precomp] = prof.th_s[precomp]
                    else:
                        # All excess stored
                        excess = 0

                if excess > 0:
                    # Any leftover water not stored becomes runoff
                    Runoff = Runoff + excess

        # Infiltration left to store after bottom compartment becomes deep
        # percolation (mm)
        DeepPerc = ToStore
    else:
        # No infiltration
        DeepPerc = 0
        Runoff = 0

    ## Update total runoff ##
    Runoff = Runoff + RunoffIni

    ## Update surface storage (if bunds are present) ##
    if Runoff > RunoffIni:
        if FieldMngt_Bunds:
            if FieldMngt_zBund > 0.001:
                # Increase surface storage
                NewCond_SurfaceStorage = NewCond_SurfaceStorage + (Runoff - RunoffIni)
                # Limit surface storage to bund height
                if NewCond_SurfaceStorage > FieldMngt_zBund:
                    # Additonal water above top of bunds becomes runoff
                    Runoff = RunoffIni + (NewCond_SurfaceStorage - FieldMngt_zBund)
                    # Set surface storage to bund height
                    NewCond_SurfaceStorage = FieldMngt_zBund
                else:
                    # No additional overtopping of bunds
                    Runoff = RunoffIni

    ## Store updated water contents ##
    NewCond_th = thnew

    ## Update deep percolation, surface runoff, and infiltration values ##
    DeepPerc = DeepPerc + DeepPerc0
    Infl = Infl - Runoff
    RunoffTot = Runoff + Runoff0

    return NewCond_th,NewCond_SurfaceStorage, DeepPerc, RunoffTot, Infl, FluxOut

aquacrop.solution.irrigation

irrigation(IrrMngt_IrrMethod, IrrMngt_SMT, IrrMngt_AppEff, IrrMngt_MaxIrr, IrrMngt_IrrInterval, IrrMngt_Schedule, IrrMngt_depth, IrrMngt_MaxIrrSeason, NewCond_GrowthStage, NewCond_IrrCum, NewCond_Epot, NewCond_Tpot, NewCond_Zroot, NewCond_th, NewCond_DAP, NewCond_TimeStepCounter, Crop, prof, Soil_zTop, growing_season, Rain, Runoff)

Function to get irrigation depth for current day

Reference Manual: irrigation description (pg. 31-32)

Arguments:

IrrMngt_IrrMethod (int): irrigation method

IrrMngt_SMT (numpy.array): soil-moisture thresholds

IrrMngt_AppEff (float): application efficiency

IrrMngt_MaxIrr (float): max irrigation depth per event

IrrMngt_IrrInterval (int): irrigation event interval (days)

IrrMngt_Schedule (numpy.array): irrigation depth schedule

IrrMngt_depth (float): depth to apply next day

IrrMngt_MaxIrrSeason (float): max irrigation for the season

NewCond_GrowthStage (float): crop growth stage

NewCond_IrrCum (float): irrigation applied so far this season

NewCond_Epot (float): potential evaporation

NewCond_Tpot (float): potential transpiration

NewCond_Zroot (float): rooting depth

NewCond_th (numpy.array): soil water content

NewCond_DAP (int): days after planting

NewCond_TimeStepCounter (int): current simulation timestep

Crop (CropStructNT): Crop paramaters

Soil (SoilProfileNT): Soil object containing soil paramaters

growing_season (bool): is growing season (True or Flase)

Rain (float): daily precipitation mm

Runoff (float): surface runoff on current day

Returns:

NewCond_Depletion (float): soil water depletion

NewCond_TAW (float): total available water

NewCond_IrrCum (float): total irrigation aplpied so far

Irr (float): Irrigaiton applied on current day mm
Source code in aquacrop/solution/irrigation.py 
 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
def irrigation(
    IrrMngt_IrrMethod: int,
    IrrMngt_SMT: float,
    IrrMngt_AppEff: float,
    IrrMngt_MaxIrr: float,
    IrrMngt_IrrInterval: int,
    IrrMngt_Schedule: "ndarray",
    IrrMngt_depth: float,
    IrrMngt_MaxIrrSeason: float,
    NewCond_GrowthStage: float,
    NewCond_IrrCum: float,
    NewCond_Epot: float,
    NewCond_Tpot: float,
    NewCond_Zroot: float,
    NewCond_th: "ndarray",
    NewCond_DAP: int,
    NewCond_TimeStepCounter: int,
    Crop: "CropStructNT",
    prof: "SoilProfileNT",
    Soil_zTop: float,
    growing_season: bool,
    Rain: float,
    Runoff: float,
    ) -> Tuple[float,float,float, float]:
    """
    Function to get irrigation depth for current day


    <a href="https://www.fao.org/3/BR246E/br246e.pdf#page=40" target="_blank">Reference Manual: irrigation description</a> (pg. 31-32)


    Arguments:


        IrrMngt_IrrMethod (int): irrigation method

        IrrMngt_SMT (numpy.array): soil-moisture thresholds

        IrrMngt_AppEff (float): application efficiency

        IrrMngt_MaxIrr (float): max irrigation depth per event

        IrrMngt_IrrInterval (int): irrigation event interval (days)

        IrrMngt_Schedule (numpy.array): irrigation depth schedule

        IrrMngt_depth (float): depth to apply next day

        IrrMngt_MaxIrrSeason (float): max irrigation for the season

        NewCond_GrowthStage (float): crop growth stage

        NewCond_IrrCum (float): irrigation applied so far this season

        NewCond_Epot (float): potential evaporation

        NewCond_Tpot (float): potential transpiration

        NewCond_Zroot (float): rooting depth

        NewCond_th (numpy.array): soil water content

        NewCond_DAP (int): days after planting

        NewCond_TimeStepCounter (int): current simulation timestep

        Crop (CropStructNT): Crop paramaters

        Soil (SoilProfileNT): Soil object containing soil paramaters

        growing_season (bool): is growing season (True or Flase)

        Rain (float): daily precipitation mm

        Runoff (float): surface runoff on current day


    Returns:

        NewCond_Depletion (float): soil water depletion

        NewCond_TAW (float): total available water

        NewCond_IrrCum (float): total irrigation aplpied so far

        Irr (float): Irrigaiton applied on current day mm


"""
    ## Store intial conditions for updating ##
    # NewCond = InitCond

    ## Determine irrigation depth (mm/day) to be applied ##
    if growing_season == True:
        # Calculate root zone water content and depletion
        # TAW_ = TAW()
        # Dr_ = Dr()
        # thRZ = RootZoneWater()
        (
            WrAct,
            Dr_Zt,
            Dr_Rz,
            TAW_Zt,
            TAW_Rz,
            thRZ_Act,
            thRZ_S,
            thRZ_FC,
            thRZ_WP,
            thRZ_Dry,
            thRZ_Aer,
        ) = root_zone_water(
            prof,
            float(NewCond_Zroot),
            NewCond_th,
            Soil_zTop,
            float(Crop.Zmin),
            Crop.Aer,
        )
        # WrAct,Dr_,TAW_,thRZ = root_zone_water(prof,float(NewCond.z_root),NewCond.th,Soil_zTop,float(Crop.Zmin),Crop.Aer)
        # Use root zone depletions and taw only for triggering irrigation
        Dr = Dr_Rz
        taw = TAW_Rz

        # Determine adjustment for inflows and outflows on current day #
        if thRZ_Act > thRZ_FC:
            rootdepth = max(NewCond_Zroot, Crop.Zmin)
            AbvFc = (thRZ_Act - thRZ_FC) * 1000 * rootdepth
        else:
            AbvFc = 0

        WCadj = NewCond_Tpot + NewCond_Epot - Rain + Runoff - AbvFc

        NewCond_Depletion = Dr + WCadj
        NewCond_TAW = taw

        # Update growth stage if it is first day of a growing season
        if NewCond_DAP == 1:
            NewCond_GrowthStage = 1

        if IrrMngt_IrrMethod == 0:
            Irr = 0

        elif IrrMngt_IrrMethod == 1:

            Dr = NewCond_Depletion / NewCond_TAW
            index = int(NewCond_GrowthStage) - 1

            if Dr > 1 - IrrMngt_SMT[index] / 100:
                # Irrigation occurs
                IrrReq = max(0, NewCond_Depletion)
                # Adjust irrigation requirements for application efficiency
                EffAdj = ((100 - IrrMngt_AppEff) + 100) / 100
                IrrReq = IrrReq * EffAdj
                # Limit irrigation to maximum depth
                Irr = min(IrrMngt_MaxIrr, IrrReq)
                # Irr = 15 # hard-code in 15mm depth for tests
            else:
                Irr = 0

        elif IrrMngt_IrrMethod == 2:  # Irrigation - fixed interval

            Dr = NewCond_Depletion

            # Get number of days in growing season so far (subtract 1 so that
            # always irrigate first on day 1 of each growing season)
            nDays = NewCond_DAP - 1

            if nDays % IrrMngt_IrrInterval == 0:
                # Irrigation occurs
                IrrReq = max(0, Dr)
                # Adjust irrigation requirements for application efficiency
                EffAdj = ((100 - IrrMngt_AppEff) + 100) / 100
                IrrReq = IrrReq * EffAdj
                # Limit irrigation to maximum depth
                Irr = min(IrrMngt_MaxIrr, IrrReq)
            else:
                # No irrigation
                Irr = 0

        elif IrrMngt_IrrMethod == 3:  # Irrigation - pre-defined schedule
            # Get current date
            idx = NewCond_TimeStepCounter
            # Find irrigation value corresponding to current date
            Irr = IrrMngt_Schedule[idx]

            assert Irr >= 0

            Irr = min(IrrMngt_MaxIrr, Irr)

        elif IrrMngt_IrrMethod == 4:  # Irrigation - net irrigation
            # Net irrigation calculation performed after transpiration, so
            # irrigation is zero here

            Irr = 0

        elif IrrMngt_IrrMethod == 5:  # depth applied each day (usually specified outside of model)

            Irr = min(IrrMngt_MaxIrr, IrrMngt_depth)

        #         else:
        #             assert 1 ==2, f'somethings gone wrong in irrigation method:{IrrMngt.irrigation_method}'

        Irr = max(0, Irr)

    elif growing_season == False:
        # No irrigation outside growing season
        Irr = 0.
        NewCond_IrrCum = 0.
        NewCond_Depletion = 0.
        NewCond_TAW = 0.


    if NewCond_IrrCum + Irr > IrrMngt_MaxIrrSeason:
        Irr = max(0, IrrMngt_MaxIrrSeason - NewCond_IrrCum)

    # Update cumulative irrigation counter for growing season
    NewCond_IrrCum = NewCond_IrrCum + Irr

    return NewCond_Depletion,NewCond_TAW,NewCond_IrrCum, Irr

aquacrop.solution.pre_irrigation

pre_irrigation(prof, Crop, InitCond, growing_season, IrrMngt)

Function to calculate pre-irrigation when in net irrigation mode

Reference Manual: Net irrigation description (pg. 31)

Arguments:

prof (SoilProfile): Soil object containing soil paramaters

Crop (CropStruct): Crop object containing Crop paramaters

InitCond (InitialCondition): InitCond object containing model paramaters

growing_season (bool): is growing season (True or Flase)

IrrMngt (IrrMngtStruct): object containing irrigation management paramaters

Returns:

NewCond (InitialCondition): InitCond object containing updated model paramaters

PreIrr (float): Pre-Irrigaiton applied on current day mm
Source code in aquacrop/solution/pre_irrigation.py 
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
def pre_irrigation(
    prof: "SoilProfileNT",
    Crop: "CropStructNT",
    InitCond: "InitialCondition",
    growing_season: bool,
    IrrMngt: "IrrMngtStruct",
    ) -> Tuple["InitialCondition", float]:
    """
    Function to calculate pre-irrigation when in net irrigation mode

    <a href="https://www.fao.org/3/BR246E/br246e.pdf#page=40" target="_blank">Reference Manual: Net irrigation description</a> (pg. 31)


    Arguments:

        prof (SoilProfile): Soil object containing soil paramaters

        Crop (CropStruct): Crop object containing Crop paramaters

        InitCond (InitialCondition): InitCond object containing model paramaters

        growing_season (bool): is growing season (True or Flase)

        IrrMngt (IrrMngtStruct): object containing irrigation management paramaters



    Returns:

        NewCond (InitialCondition): InitCond object containing updated model paramaters

        PreIrr (float): Pre-Irrigaiton applied on current day mm




    """
    # Store initial conditions for updating ##
    NewCond = InitCond

    ## Calculate pre-irrigation needs ##
    if growing_season == True:
        if (IrrMngt.irrigation_method != 4) or (NewCond.dap != 1):
            # No pre-irrigation as not in net irrigation mode or not on first day
            # of the growing season
            PreIrr = 0
        else:
            # Determine compartments covered by the root zone
            rootdepth = round(max(NewCond.z_root, Crop.Zmin), 2)

            compRz = np.argwhere(prof.dzsum >= rootdepth).flatten()[0]

            PreIrr = 0
            for ii in range(int(compRz)):

                # Determine critical water content threshold
                thCrit = prof.th_wp[ii] + (
                    (IrrMngt.NetIrrSMT / 100) * (prof.th_fc[ii] - prof.th_wp[ii])
                )

                # Check if pre-irrigation is required
                if NewCond.th[ii] < thCrit:
                    PreIrr = PreIrr + ((thCrit - NewCond.th[ii]) * 1000 * prof.dz[ii])
                    NewCond.th[ii] = thCrit

    else:
        PreIrr = 0

    return NewCond, PreIrr

aquacrop.solution.rainfall_partition

rainfall_partition(precipitation, InitCond_th, NewCond_DaySubmerged, FieldMngt_SRinhb, FieldMngt_Bunds, FieldMngt_zBund, FieldMngt_CNadjPct, Soil_CN, Soil_AdjCN, Soil_zCN, Soil_nComp, prof)

Function to partition rainfall into surface runoff and infiltration using the curve number approach

Reference Manual: rainfall partition calculations (pg. 48-51)

Arguments:

precipitation (float): Percipitation on current day

InitCond_th (numpy.array): InitCond object containing model paramaters

NewCond_DaySubmerged (int): number of days submerged

FieldMngt_SRinhb (float): field management params

FieldMngt_Bunds (bool): field management params

FieldMngt_zBund (float): bund height

FieldMngt_CNadjPct (float): curve number adjustment percent

Soil_CN (float): curve number

Soil_AdjCN (float): adjusted curve number

Soil_zCN (float` :

Soil_nComp (float): number of compartments

prof (SoilProfile): Soil object

Returns:

Runoff (float): Total Suface Runoff

Infl (float): Total Infiltration

NewCond_DaySubmerged (float): number of days submerged
Source code in aquacrop/solution/rainfall_partition.py 
 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
@cc.export("rainfall_partition", (f8,f8[:],i8,f8,b1,f8,f8,f8,f8,f8,f8,SoilProfileNT_typ_sig))
def rainfall_partition(
    precipitation: float,
    InitCond_th: "ndarray",
    NewCond_DaySubmerged: int,
    FieldMngt_SRinhb: float,
    FieldMngt_Bunds: bool,
    FieldMngt_zBund: float,
    FieldMngt_CNadjPct: float,
    Soil_CN: float,
    Soil_AdjCN: float,
    Soil_zCN: float,
    Soil_nComp: int,
    prof: "SoilProfileNT",
) -> Tuple[float, float, float]:
    """
    Function to partition rainfall into surface runoff and infiltration using the curve number approach


    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=57" target="_blank">Reference Manual: rainfall partition calculations</a> (pg. 48-51)


    Arguments:

        precipitation (float): Percipitation on current day

        InitCond_th (numpy.array): InitCond object containing model paramaters

        NewCond_DaySubmerged (int): number of days submerged

        FieldMngt_SRinhb (float): field management params

        FieldMngt_Bunds (bool): field management params

        FieldMngt_zBund (float): bund height

        FieldMngt_CNadjPct (float): curve number adjustment percent

        Soil_CN (float): curve number

        Soil_AdjCN (float): adjusted curve number

        Soil_zCN (float` :

        Soil_nComp (float): number of compartments

        prof (SoilProfile): Soil object

    Returns:

        Runoff (float): Total Suface Runoff

        Infl (float): Total Infiltration

        NewCond_DaySubmerged (float): number of days submerged


    """

    # can probs make this faster by doing a if precipitation=0 loop

    ## Store initial conditions for updating ##
    # NewCond = InitCond

    ## Calculate runoff ##
    if (FieldMngt_SRinhb == False) and ((FieldMngt_Bunds == False) or (FieldMngt_zBund < 0.001)):
        # Surface runoff is not inhibited and no soil bunds are on field
        # Reset submerged days
        NewCond_DaySubmerged = 0
        # Adjust curve number for field management practices
        cn = Soil_CN * (1 + (FieldMngt_CNadjPct / 100))
        if Soil_AdjCN == 1:  # Adjust cn for antecedent moisture
            # Calculate upper and lowe curve number bounds
            CNbot = round(
                1.4 * (np.exp(-14 * np.log(10)))
                + (0.507 * cn)
                - (0.00374 * cn ** 2)
                + (0.0000867 * cn ** 3)
            )
            CNtop = round(
                5.6 * (np.exp(-14 * np.log(10)))
                + (2.33 * cn)
                - (0.0209 * cn ** 2)
                + (0.000076 * cn ** 3)
            )
            # Check which compartment cover depth of top soil used to adjust
            # curve number
            comp_sto_array = prof.dzsum[prof.dzsum >= Soil_zCN]
            if comp_sto_array.shape[0] == 0:
                comp_sto = int(Soil_nComp)
            else:
                comp_sto = int(Soil_nComp - comp_sto_array.shape[0])

            comp_sto+=1

            # Calculate weighting factors by compartment
            xx = 0
            wrel = np.zeros(comp_sto)
            for ii in range(comp_sto):
                if prof.dzsum[ii] > Soil_zCN:
                    prof.dzsum[ii] = Soil_zCN

                wx = 1.016 * (1 - np.exp(-4.16 * (prof.dzsum[ii] / Soil_zCN)))
                wrel[ii] = wx - xx
                if wrel[ii] < 0:
                    wrel[ii] = 0
                elif wrel[ii] > 1:
                    wrel[ii] = 1

                xx = wx

            # Calculate relative wetness of top soil
            wet_top = 0
            # prof = prof

            for ii in range(comp_sto):
                th = max(prof.th_wp[ii], InitCond_th[ii])
                wet_top = wet_top + (
                    wrel[ii] * ((th - prof.th_wp[ii]) / (prof.th_fc[ii] - prof.th_wp[ii]))
                )

            # Calculate adjusted curve number
            if wet_top > 1:
                wet_top = 1
            elif wet_top < 0:
                wet_top = 0

            cn = round(CNbot + (CNtop - CNbot) * wet_top)

        # Partition rainfall into runoff and infiltration (mm)
        S = (25400 / cn) - 254
        term = precipitation - ((5 / 100) * S)
        if term <= 0:
            Runoff = 0
            Infl = precipitation
        else:
            Runoff = (term ** 2) / (precipitation + (1 - (5 / 100)) * S)
            Infl = precipitation - Runoff

    else:
        # bunds on field, therefore no surface runoff
        Runoff = 0
        Infl = precipitation

    return Runoff, Infl, NewCond_DaySubmerged

aquacrop.solution.root_development

root_development(Crop, prof, NewCond_DAP, NewCond_Zroot, NewCond_DelayedCDs, NewCond_GDDcum, NewCond_DelayedGDDs, NewCond_TrRatio, NewCond_th, NewCond_CC, NewCond_CC_NS, NewCond_Germination, NewCond_rCor, NewCond_Tpot, NewCond_zGW, gdd, growing_season, water_table_presence)

Function to calculate root zone expansion

Reference Manual: root developement equations (pg. 37-41)

Arguments:

Crop (CropStructNT): crop params

prof (SoilProfileNT): soilv profile paramaters

NewCond_DAP (float): days after planting

NewCond_Zroot (float): root depth

NewCond_DelayedCDs (float): delayed calendar days

NewCond_GDDcum (float): cumulative growing degree days

NewCond_TrRatio (float): transpiration ratio

NewCond_CC (float): canopy cover

NewCond_CC_NS (float): canopy cover no-stress

NewCond_Germination (float): germination flag

NewCond_rCor (float):

NewCond_DAP (float): days after planting

NewCond_Tpot (float): potential transpiration

NewCond_zGW (float): groundwater depth

gdd (float): Growing degree days on current day

growing_season (bool): is growing season (True or Flase)

water_table_presence (int): water table present (True=1 or Flase=0)

Returns:

NewCond_Zroot (float): updated rooting depth
Source code in aquacrop/solution/root_development.py 
 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
@cc.export("root_development", (CropStructNT_type_sig,SoilProfileNT_typ_sig,f8,f8,f8,f8,f8,f8,f8[:],f8,f8,b1,f8,f8,f8,f8,b1,i8))
def root_development(
    Crop: "CropStructNT",
    prof: "SoilProfileNT",
    NewCond_DAP: float,
    NewCond_Zroot: float,
    NewCond_DelayedCDs: float,
    NewCond_GDDcum: float,
    NewCond_DelayedGDDs: float,
    NewCond_TrRatio: float,
    NewCond_th: "ndarray",
    NewCond_CC: float,
    NewCond_CC_NS: float,
    NewCond_Germination: bool,
    NewCond_rCor: float,
    NewCond_Tpot: float,
    NewCond_zGW: float,
    gdd: float,
    growing_season: bool,
    water_table_presence: int,
    ) -> float:
    """
    Function to calculate root zone expansion

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=46" target="_blank">Reference Manual: root developement equations</a> (pg. 37-41)


    Arguments:

        Crop (CropStructNT): crop params

        prof (SoilProfileNT): soilv profile paramaters

        NewCond_DAP (float): days after planting

        NewCond_Zroot (float): root depth

        NewCond_DelayedCDs (float): delayed calendar days

        NewCond_GDDcum (float): cumulative growing degree days

        NewCond_TrRatio (float): transpiration ratio

        NewCond_CC (float): canopy cover

        NewCond_CC_NS (float): canopy cover no-stress

        NewCond_Germination (float): germination flag

        NewCond_rCor (float): 

        NewCond_DAP (float): days after planting

        NewCond_Tpot (float): potential transpiration

        NewCond_zGW (float): groundwater depth

        gdd (float): Growing degree days on current day

        growing_season (bool): is growing season (True or Flase)

        water_table_presence (int): water table present (True=1 or Flase=0)


    Returns:

        NewCond_Zroot (float): updated rooting depth


    """
    # Store initial conditions for updating
    # NewCond = InitCond

    # save initial zroot
    Zroot_init = float(NewCond_Zroot) * 1.0
    Soil_nLayer = np.unique(prof.Layer).shape[0]

    # Calculate root expansion (if in growing season)
    if growing_season == True:
        # If today is first day of season, root depth is equal to minimum depth
        if NewCond_DAP == 1:
            NewCond_Zroot = float(Crop.Zmin) * 1.0
            Zroot_init = float(Crop.Zmin) * 1.0

        # Adjust time for any delayed development
        if Crop.CalendarType == 1:
            tAdj = NewCond_DAP - NewCond_DelayedCDs
        elif Crop.CalendarType == 2:
            tAdj = NewCond_GDDcum - NewCond_DelayedGDDs

        # Calculate root expansion #
        Zini = Crop.Zmin * (Crop.PctZmin / 100)
        t0 = round((Crop.Emergence / 2))
        tmax = Crop.MaxRooting
        if Crop.CalendarType == 1:
            tOld = tAdj - 1
        elif Crop.CalendarType == 2:
            tOld = tAdj - gdd

        # Potential root depth on previous day
        if tOld >= tmax:
            ZrOld = Crop.Zmax
        elif tOld <= t0:
            ZrOld = Zini
        else:
            X = (tOld - t0) / (tmax - t0)
            ZrOld = Zini + (Crop.Zmax - Zini) * np.power(X, 1 / Crop.fshape_r)

        if ZrOld < Crop.Zmin:
            ZrOld = Crop.Zmin

        # Potential root depth on current day
        if tAdj >= tmax:
            Zr = Crop.Zmax
        elif tAdj <= t0:
            Zr = Zini
        else:
            X = (tAdj - t0) / (tmax - t0)
            Zr = Zini + (Crop.Zmax - Zini) * np.power(X, 1 / Crop.fshape_r)

        if Zr < Crop.Zmin:
            Zr = Crop.Zmin

        # Store Zr as potential value
        ZrPot = Zr

        # Determine rate of change
        dZr = Zr - ZrOld

        # Adjust expansion rate for presence of restrictive soil horizons
        if Zr > Crop.Zmin:
            layeri = 1
            l_idx = np.argwhere(prof.Layer == layeri).flatten()
            Zsoil = prof.dz[l_idx].sum()
            while (round(Zsoil, 2) <= Crop.Zmin) and (layeri < Soil_nLayer):
                layeri = layeri + 1
                l_idx = np.argwhere(prof.Layer == layeri).flatten()
                Zsoil = Zsoil + prof.dz[l_idx].sum()

            soil_layer_dz = prof.dz[l_idx].sum()
            layer_comp = l_idx[0]
            # soil_layer = prof.Layer[layeri]
            ZrAdj = Crop.Zmin
            ZrRemain = Zr - Crop.Zmin
            deltaZ = Zsoil - Crop.Zmin
            EndProf = False
            while EndProf == False:
                ZrTest = ZrAdj + (ZrRemain * (prof.Penetrability[layer_comp] / 100))
                if (
                    (layeri == Soil_nLayer)
                    or (prof.Penetrability[layer_comp] == 0)
                    or (ZrTest <= Zsoil)
                ):
                    ZrOUT = ZrTest
                    EndProf = True
                else:
                    ZrAdj = Zsoil
                    ZrRemain = ZrRemain - (deltaZ / (prof.Penetrability[layer_comp] / 100))
                    layeri = layeri + 1
                    l_idx = np.argwhere(prof.Layer == layeri).flatten()
                    layer_comp = l_idx[0]
                    soil_layer_dz = prof.dz[l_idx].sum()
                    Zsoil = Zsoil + soil_layer_dz
                    deltaZ = soil_layer_dz

            # Correct Zr and dZr for effects of restrictive horizons
            Zr = ZrOUT
            dZr = Zr - ZrOld

        # Adjust rate of expansion for any stomatal water stress
        if NewCond_TrRatio < 0.9999:
            if Crop.fshape_ex >= 0:
                dZr = dZr * NewCond_TrRatio
            else:
                fAdj = (np.exp(NewCond_TrRatio * Crop.fshape_ex) - 1) / (np.exp(Crop.fshape_ex) - 1)
                dZr = dZr * fAdj

        # print(NewCond.dap,NewCond.th)

        # Adjust rate of root expansion for dry soil at expansion front
        if dZr > 0.001:
            # Define water stress threshold for inhibition of root expansion
            pZexp = Crop.p_up[1] + ((1 - Crop.p_up[1]) / 2)
            # Define potential new root depth
            ZiTmp = float(Zroot_init + dZr)
            # Find compartment that root zone will expand in to
            # compi_index = prof.dzsum[prof.dzsum>=ZiTmp].index[0] # have changed to index
            idx = np.argwhere(prof.dzsum >= ZiTmp).flatten()[0]
            prof = prof
            # Get taw in compartment
            layeri = prof.Layer[idx]
            TAWprof = prof.th_fc[idx] - prof.th_wp[idx]
            # Define stress threshold
            thThr = prof.th_fc[idx] - (pZexp * TAWprof)
            # Check for stress conditions
            if NewCond_th[idx] < thThr:
                # Root expansion limited by water content at expansion front
                if NewCond_th[idx] <= prof.th_wp[idx]:

                    # Expansion fully inhibited
                    dZr = 0
                else:
                    # Expansion partially inhibited
                    Wrel = (prof.th_fc[idx] - NewCond_th[idx]) / TAWprof
                    Drel = 1 - ((1 - Wrel) / (1 - pZexp))
                    Ks = 1 - (
                        (np.exp(Drel * Crop.fshape_w[1]) - 1) / (np.exp(Crop.fshape_w[1]) - 1)
                    )
                    dZr = dZr * Ks

        # Adjust for early senescence
        if (NewCond_CC <= 0) and (NewCond_CC_NS > 0.5):
            dZr = 0

        # Adjust root expansion for failure to germinate (roots cannot expand
        # if crop has not germinated)
        if NewCond_Germination == False:
            dZr = 0

        # Get new rooting depth
        NewCond_Zroot = float(Zroot_init + dZr)

        # Adjust root density if deepening is restricted due to dry subsoil
        # and/or restrictive layers
        if NewCond_Zroot < ZrPot:
            NewCond_rCor = (
                2 * (ZrPot / NewCond_Zroot) * ((Crop.SxTop + Crop.SxBot) / 2) - Crop.SxTop
            ) / Crop.SxBot

            if NewCond_Tpot > 0:
                NewCond_rCor = NewCond_rCor * NewCond_TrRatio
                if NewCond_rCor < 1:
                    NewCond_rCor = 1

        else:
            NewCond_rCor = 1

        # Limit rooting depth if groundwater table is present (roots cannot
        # develop below the water table)
        if (water_table_presence == 1) and (NewCond_zGW > 0):
            if NewCond_Zroot > NewCond_zGW:
                NewCond_Zroot = float(NewCond_zGW)
                if NewCond_Zroot < Crop.Zmin:
                    NewCond_Zroot = float(Crop.Zmin)

    else:
        # No root system outside of the growing season
        NewCond_Zroot = 0

    return NewCond_Zroot, NewCond_rCor

aquacrop.solution.root_zone_water

root_zone_water(prof, InitCond_Zroot, InitCond_th, Soil_zTop, Crop_Zmin, Crop_Aer)

Function to calculate actual and total available water in the rootzone at current time step

Reference Manual: root-zone water calculations (pg. 5-8)

Arguments:

prof (SoilProfile): jit class Object containing soil paramaters

InitCond_Zroot (float): Initial rooting depth

InitCond_th (np.array): Initial water content

Soil_zTop (float): Top soil depth

Crop_Zmin (float): crop minimum rooting depth

Crop_Aer (int): number of aeration stress days

Returns:

WrAct (float):  Actual rootzone water content

Dr_Zt (float):  topsoil depletion

Dr_Rz (float):  rootzone depletion

TAW_Zt (float):  topsoil total available water

TAW_Rz (float):  rootzone total available water

thRZ_Act (float):  Actual rootzone water content

thRZ_S (float):  rootzone water content at saturation

thRZ_FC (float):  rootzone water content at field capacity

thRZ_WP (float):  rootzone water content at wilting point

thRZ_Dry (float):  rootzone water content at air dry

thRZ_Aer (float):  rootzone water content at aeration stress threshold
Source code in aquacrop/solution/root_zone_water.py 
 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
@njit
@cc.export("root_zone_water", (SoilProfileNT_typ_sig,f8,f8[:],f8,f8,f8))
def root_zone_water(
    prof: "SoilProfileNT_typ_sig",
    InitCond_Zroot: float,
    InitCond_th: "ndarray",
    Soil_zTop: float,
    Crop_Zmin: float,
    Crop_Aer: float,
) -> Tuple[float, float, float, float, float, float, float, float, float, float, float]:
    """
    Function to calculate actual and total available water in the rootzone at current time step


    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=14" target="_blank">Reference Manual: root-zone water calculations</a> (pg. 5-8)


    Arguments:

        prof (SoilProfile): jit class Object containing soil paramaters

        InitCond_Zroot (float): Initial rooting depth

        InitCond_th (np.array): Initial water content

        Soil_zTop (float): Top soil depth

        Crop_Zmin (float): crop minimum rooting depth

        Crop_Aer (int): number of aeration stress days

    Returns:

        WrAct (float):  Actual rootzone water content

        Dr_Zt (float):  topsoil depletion

        Dr_Rz (float):  rootzone depletion

        TAW_Zt (float):  topsoil total available water

        TAW_Rz (float):  rootzone total available water

        thRZ_Act (float):  Actual rootzone water content

        thRZ_S (float):  rootzone water content at saturation

        thRZ_FC (float):  rootzone water content at field capacity

        thRZ_WP (float):  rootzone water content at wilting point

        thRZ_Dry (float):  rootzone water content at air dry

        thRZ_Aer (float):  rootzone water content at aeration stress threshold 


    """

    ## Calculate root zone water content and available water ##
    # Compartments covered by the root zone
    rootdepth = round(np.maximum(InitCond_Zroot, Crop_Zmin), 2)
    comp_sto = np.argwhere(prof.dzsum >= rootdepth).flatten()[0]

    # Initialise counters
    WrAct = 0
    WrS = 0
    WrFC = 0
    WrWP = 0
    WrDry = 0
    WrAer = 0
    for ii in range(comp_sto + 1):
        # Fraction of compartment covered by root zone
        if prof.dzsum[ii] > rootdepth:
            factor = 1 - ((prof.dzsum[ii] - rootdepth) / prof.dz[ii])
        else:
            factor = 1

        # Actual water storage in root zone (mm)
        WrAct = WrAct + round(factor * 1000 * InitCond_th[ii] * prof.dz[ii], 2)
        # Water storage in root zone at saturation (mm)
        WrS = WrS + round(factor * 1000 * prof.th_s[ii] * prof.dz[ii], 2)
        # Water storage in root zone at field capacity (mm)
        WrFC = WrFC + round(factor * 1000 * prof.th_fc[ii] * prof.dz[ii], 2)
        # Water storage in root zone at permanent wilting point (mm)
        WrWP = WrWP + round(factor * 1000 * prof.th_wp[ii] * prof.dz[ii], 2)
        # Water storage in root zone at air dry (mm)
        WrDry = WrDry + round(factor * 1000 * prof.th_dry[ii] * prof.dz[ii], 2)
        # Water storage in root zone at aeration stress threshold (mm)
        WrAer = WrAer + round(factor * 1000 * (prof.th_s[ii] - (Crop_Aer / 100)) * prof.dz[ii], 2)

    if WrAct < 0:
        WrAct = 0

    # define total available water, depletion, root zone water content
    # taw = TAW()
    # Dr = Dr()
    # thRZ = RootZoneWater()

    # Calculate total available water (m3/m3)
    TAW_Rz = max(WrFC - WrWP, 0.0)
    # Calculate soil water depletion (mm)
    Dr_Rz = min(WrFC - WrAct, TAW_Rz)

    # Actual root zone water content (m3/m3)
    thRZ_Act = WrAct / (rootdepth * 1000)
    # Root zone water content at saturation (m3/m3)
    thRZ_S = WrS / (rootdepth * 1000)
    # Root zone water content at field capacity (m3/m3)
    thRZ_FC = WrFC / (rootdepth * 1000)
    # Root zone water content at permanent wilting point (m3/m3)
    thRZ_WP = WrWP / (rootdepth * 1000)
    # Root zone water content at air dry (m3/m3)
    thRZ_Dry = WrDry / (rootdepth * 1000)
    # Root zone water content at aeration stress threshold (m3/m3)
    thRZ_Aer = WrAer / (rootdepth * 1000)

    # print('inside')

    # thRZ = thRZNT(
    # Act=thRZ_Act,
    # S=thRZ_S,
    # FC=thRZ_FC,
    # WP=thRZ_WP,
    # Dry=thRZ_Dry,
    # Aer=thRZ_Aer,
    # )
    # print(thRZ)


    ## Calculate top soil water content and available water ##
    if rootdepth > Soil_zTop:
        # Determine compartments covered by the top soil
        ztopdepth = round(Soil_zTop, 2)
        comp_sto = np.sum(prof.dzsum <= ztopdepth)
        # Initialise counters
        WrAct_Zt = 0
        WrFC_Zt = 0
        WrWP_Zt = 0
        # Calculate water storage in top soil
        assert comp_sto > 0

        for ii in range(comp_sto):

            # Fraction of compartment covered by root zone
            if prof.dzsum[ii] > ztopdepth:
                factor = 1 - ((prof.dzsum[ii] - ztopdepth) / prof.dz[ii])
            else:
                factor = 1

            # Actual water storage in top soil (mm)
            WrAct_Zt = WrAct_Zt + (factor * 1000 * InitCond_th[ii] * prof.dz[ii])
            # Water storage in top soil at field capacity (mm)
            WrFC_Zt = WrFC_Zt + (factor * 1000 * prof.th_fc[ii] * prof.dz[ii])
            # Water storage in top soil at permanent wilting point (mm)
            WrWP_Zt = WrWP_Zt + (factor * 1000 * prof.th_wp[ii] * prof.dz[ii])

        # Ensure available water in top soil is not less than zero
        if WrAct_Zt < 0:
            WrAct_Zt = 0

        # Calculate total available water in top soil (m3/m3)
        TAW_Zt = max(WrFC_Zt - WrWP_Zt, 0)
        # Calculate depletion in top soil (mm)
        Dr_Zt = min(WrFC_Zt - WrAct_Zt, TAW_Zt)
    else:
        # Set top soil depletions and taw to root zone values
        Dr_Zt = Dr_Rz
        TAW_Zt = TAW_Rz


    return (
        WrAct,
        Dr_Zt,
        Dr_Rz,
        TAW_Zt,
        TAW_Rz,
        thRZ_Act,
        thRZ_S,
        thRZ_FC,
        thRZ_WP,
        thRZ_Dry,
        thRZ_Aer,
    )

aquacrop.solution.soil_evaporation

soil_evaporation(ClockStruct_EvapTimeSteps, ClockStruct_SimOffSeason, ClockStruct_TimeStepCounter, prof, Soil_EvapZmin, Soil_EvapZmax, Soil_REW, Soil_Kex, Soil_fwcc, Soil_fWrelExp, Soil_fevap, Crop_CalendarType, Crop_Senescence, IrrMngt_IrrMethod, IrrMngt_WetSurf, FieldMngt_Mulches, FieldMngt_fMulch, FieldMngt_MulchPct, NewCond_DAP, NewCond_Wsurf, NewCond_EvapZ, NewCond_Stage2, NewCond_th, NewCond_DelayedCDs, NewCond_GDDcum, NewCond_DelayedGDDs, NewCond_CCxW, NewCond_CCadj, NewCond_CCxAct, NewCond_CC, NewCond_PrematSenes, NewCond_SurfaceStorage, NewCond_Wstage2, NewCond_Epot, et0, Infl, Rain, Irr, growing_season)

Function to calculate daily soil evaporation

Reference Manual: evaporation equations (pg. 73-81)

Arguments:

ClockStruct_EvapTimeSteps (int): number of evaportation time steps

ClockStruct_SimOffSeason (bool): simulate off season? (False=no, True=yes)

ClockStruct_TimeStepCounter (int): time step counter

prof (SoilProfileNT): soil profile object

Soil_EvapZmin (float): minimum evaporation depth (m)

Soil_EvapZmax (float): maximum evaporation depth (m)

Soil_REW (float): Readily Evaporable Water

Soil_Kex (float): Soil evaporation coefficient

Soil_fwcc (float):

Soil_fWrelExp (float):

Soil_fevap (float):

Crop_CalendarType (int): calendar type

Crop_Senescence (float):

IrrMngt_IrrMethod (int): irrigation method

IrrMngt_WetSurf (float): wet surface area

FieldMngt_Mulches (bool): mulch present? (0=no, 1=yes)

FieldMngt_fMulch (float): mulch factor

FieldMngt_MulchPct (float): mulch percentage

NewCond_DAP (int): days after planting

NewCond_Wsurf (float): wet surface area

NewCond_EvapZ (float): evaporation depth (m)

NewCond_Stage2 (float): stage 2 evaporation

NewCond_th (ndarray): soil water content

NewCond_DelayedCDs: delayed calendar days

NewCond_GDDcum (float): cumulative growing degree days

NewCond_DelayedGDDs (float): delayed growing degree days

NewCond_CCxW (float):

NewCond_CCadj (float): canopy cover adjusted

NewCond_CCxAct: max canopy cover actual

NewCond_CC (float): canopy cover

NewCond_PrematSenes (bool): prematurity senescence? (0=no, 1=yes)

NewCond_SurfaceStorage (float): surface storage

NewCond_Wstage2 (float): stage 2 water content

NewCond_Epot (float): potential evaporation

et0 (float): daily reference evapotranspiration

Infl (float): Infiltration on current day

Rain (float): daily precipitation mm

Irr (float): Irrigation applied on current day

growing_season (bool): is growing season (True or Flase)

Returns:

NewCond_Epot (float): Potential surface evaporation current day

NewCond_th (ndarray): updated soil water content

NewCond_Stage2 (bool): stage 2 soil evaporation

NewCond_Wstage2 (float): stage 2 soil evaporation

NewCond_Wsurf (float): updated surface water content

NewCond_SurfaceStorage (float): updated surface storage

NewCond_EvapZ (float): updated evaporation layer depth

EsAct (float): Actual surface evaporation current day

EsPot (float): Potential surface evaporation current day
Source code in aquacrop/solution/soil_evaporation.py 
 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
@cc.export(
    "soil_evaporation", (i8,i8,i8,SoilProfileNT_typ_sig,
    f8,f8,f8,f8,f8,f8,f8,i8,f8,i8,f8,b1,f8,f8,i8,f8,f8,f8,f8[:],f8,f8,f8,f8,f8,f8,
        f8,b1,f8,f8,f8,f8,f8,f8,f8,b1),
)
def soil_evaporation(
    ClockStruct_EvapTimeSteps: int,
    ClockStruct_SimOffSeason: bool,
    ClockStruct_TimeStepCounter: int,
    prof: "SoilProfileNT",
    Soil_EvapZmin: float,
    Soil_EvapZmax: float,
    Soil_REW: float,
    Soil_Kex: float,
    Soil_fwcc: float,
    Soil_fWrelExp: float,
    Soil_fevap: float,
    Crop_CalendarType: int,
    Crop_Senescence: float,
    IrrMngt_IrrMethod: int,
    IrrMngt_WetSurf: float,
    FieldMngt_Mulches: bool,
    FieldMngt_fMulch: float,
    FieldMngt_MulchPct: float,
    NewCond_DAP: int,
    NewCond_Wsurf: float,
    NewCond_EvapZ: float,
    NewCond_Stage2: float,
    NewCond_th: "ndarray",
    NewCond_DelayedCDs: float,
    NewCond_GDDcum: float,
    NewCond_DelayedGDDs: float,
    NewCond_CCxW: float,
    NewCond_CCadj: float,
    NewCond_CCxAct: float,
    NewCond_CC: float,
    NewCond_PrematSenes: bool,
    NewCond_SurfaceStorage: float,
    NewCond_Wstage2: float,
    NewCond_Epot: float,
    et0: float,
    Infl: float,
    Rain: float,
    Irr: float,
    growing_season: bool,
) -> Tuple[float, "ndarray", bool, float, float, float, float, float, float]:

    """
    Function to calculate daily soil evaporation

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=82" target="_blank">Reference Manual: evaporation equations</a> (pg. 73-81)


    Arguments:

        ClockStruct_EvapTimeSteps (int): number of evaportation time steps

        ClockStruct_SimOffSeason (bool): simulate off season? (False=no, True=yes)

        ClockStruct_TimeStepCounter (int): time step counter

        prof (SoilProfileNT): soil profile object

        Soil_EvapZmin (float): minimum evaporation depth (m)

        Soil_EvapZmax (float): maximum evaporation depth (m)

        Soil_REW (float): Readily Evaporable Water 

        Soil_Kex (float): Soil evaporation coefficient

        Soil_fwcc (float): 

        Soil_fWrelExp (float): 

        Soil_fevap (float): 

        Crop_CalendarType (int): calendar type 

        Crop_Senescence (float):

        IrrMngt_IrrMethod (int): irrigation method

        IrrMngt_WetSurf (float): wet surface area

        FieldMngt_Mulches (bool): mulch present? (0=no, 1=yes)

        FieldMngt_fMulch (float): mulch factor

        FieldMngt_MulchPct (float): mulch percentage

        NewCond_DAP (int): days after planting

        NewCond_Wsurf (float): wet surface area

        NewCond_EvapZ (float): evaporation depth (m)

        NewCond_Stage2 (float): stage 2 evaporation

        NewCond_th (ndarray): soil water content

        NewCond_DelayedCDs: delayed calendar days

        NewCond_GDDcum (float): cumulative growing degree days

        NewCond_DelayedGDDs (float): delayed growing degree days

        NewCond_CCxW (float): 

        NewCond_CCadj (float): canopy cover adjusted

        NewCond_CCxAct: max canopy cover actual

        NewCond_CC (float): canopy cover

        NewCond_PrematSenes (bool): prematurity senescence? (0=no, 1=yes)

        NewCond_SurfaceStorage (float): surface storage

        NewCond_Wstage2 (float): stage 2 water content

        NewCond_Epot (float): potential evaporation

        et0 (float): daily reference evapotranspiration

        Infl (float): Infiltration on current day

        Rain (float): daily precipitation mm

        Irr (float): Irrigation applied on current day

        growing_season (bool): is growing season (True or Flase)


    Returns:

        NewCond_Epot (float): Potential surface evaporation current day

        NewCond_th (ndarray): updated soil water content

        NewCond_Stage2 (bool): stage 2 soil evaporation

        NewCond_Wstage2 (float): stage 2 soil evaporation 

        NewCond_Wsurf (float): updated surface water content

        NewCond_SurfaceStorage (float): updated surface storage

        NewCond_EvapZ (float): updated evaporation layer depth

        EsAct (float): Actual surface evaporation current day

        EsPot (float): Potential surface evaporation current day

    """

    # Wevap = WaterEvaporation()

    ## Store initial conditions in new structure that will be updated ##
    # NewCond = InitCond

    ## Prepare stage 2 evaporation (rew gone) ##
    # Only do this if it is first day of simulation, or if it is first day of
    # growing season and not simulating off-season
    if (ClockStruct_TimeStepCounter == 0) or (
        (NewCond_DAP == 1) and (ClockStruct_SimOffSeason == False)
    ):
        # Reset storage in surface soil layer to zero
        NewCond_Wsurf = 0
        # Set evaporation depth to minimum
        NewCond_EvapZ = Soil_EvapZmin
        # Trigger stage 2 evaporation
        NewCond_Stage2 = True
        # Get relative water content for start of stage 2 evaporation
        Wevap_Sat, Wevap_Fc, Wevap_Wp, Wevap_Dry, Wevap_Act = evap_layer_water_content(
            NewCond_th,
            NewCond_EvapZ,
            prof,
        )
        NewCond_Wstage2 = round(
            (Wevap_Act - (Wevap_Fc - Soil_REW)) / (Wevap_Sat - (Wevap_Fc - Soil_REW)), 2
        )
        if NewCond_Wstage2 < 0:
            NewCond_Wstage2 = 0

    ## Prepare soil evaporation stage 1 ##
    # Adjust water in surface evaporation layer for any infiltration
    if (Rain > 0) or ((Irr > 0) and (IrrMngt_IrrMethod != 4)):
        # Only prepare stage one when rainfall occurs, or when irrigation is
        # trigerred (not in net irrigation mode)
        if Infl > 0:
            # Update storage in surface evaporation layer for incoming
            # infiltration
            NewCond_Wsurf = Infl
            # Water stored in surface evaporation layer cannot exceed rew
            if NewCond_Wsurf > Soil_REW:
                NewCond_Wsurf = Soil_REW

            # Reset variables
            NewCond_Wstage2 = 0
            NewCond_EvapZ = Soil_EvapZmin
            NewCond_Stage2 = False

    ## Calculate potential soil evaporation rate (mm/day) ##
    if growing_season == True:
        # Adjust time for any delayed development
        if Crop_CalendarType == 1:
            tAdj = NewCond_DAP - NewCond_DelayedCDs
        elif Crop_CalendarType == 2:
            tAdj = NewCond_GDDcum - NewCond_DelayedGDDs

        # Calculate maximum potential soil evaporation
        EsPotMax = Soil_Kex * et0 * (1 - NewCond_CCxW * (Soil_fwcc / 100))
        # Calculate potential soil evaporation (given current canopy cover
        # size)
        EsPot = Soil_Kex * (1 - NewCond_CCadj) * et0

        # Adjust potential soil evaporation for effects of withered canopy
        if (tAdj > Crop_Senescence) and (NewCond_CCxAct > 0):
            if NewCond_CC > (NewCond_CCxAct / 2):
                if NewCond_CC > NewCond_CCxAct:
                    mult = 0
                else:
                    mult = (NewCond_CCxAct - NewCond_CC) / (NewCond_CCxAct / 2)

            else:
                mult = 1

            EsPot = EsPot * (1 - NewCond_CCxAct * (Soil_fwcc / 100) * mult)
            CCxActAdj = (
                (1.72 * NewCond_CCxAct) - (NewCond_CCxAct ** 2) + 0.3 * (NewCond_CCxAct ** 3)
            )
            EsPotMin = Soil_Kex * (1 - CCxActAdj) * et0
            if EsPotMin < 0:
                EsPotMin = 0

            if EsPot < EsPotMin:
                EsPot = EsPotMin
            elif EsPot > EsPotMax:
                EsPot = EsPotMax

        if NewCond_PrematSenes == True:
            if EsPot > EsPotMax:
                EsPot = EsPotMax

    else:
        # No canopy cover outside of growing season so potential soil
        # evaporation only depends on reference evapotranspiration
        EsPot = Soil_Kex * et0

    ## Adjust potential soil evaporation for mulches and/or partial wetting ##
    # mulches
    if NewCond_SurfaceStorage < 0.000001:
        if not FieldMngt_Mulches:
            # No mulches present
            EsPotMul = EsPot
        elif FieldMngt_Mulches:
            # mulches present
            EsPotMul = EsPot * (1 - FieldMngt_fMulch * (FieldMngt_MulchPct / 100))

    else:
        # Surface is flooded - no adjustment of potential soil evaporation for
        # mulches
        EsPotMul = EsPot

    # Partial surface wetting by irrigation
    if (Irr > 0) and (IrrMngt_IrrMethod != 4):
        # Only apply adjustment if irrigation occurs and not in net irrigation
        # mode
        if (Rain > 1) or (NewCond_SurfaceStorage > 0):
            # No adjustment for partial wetting - assume surface is fully wet
            EsPotIrr = EsPot
        else:
            # Adjust for proprtion of surface area wetted by irrigation
            EsPotIrr = EsPot * (IrrMngt_WetSurf / 100)

    else:
        # No adjustment for partial surface wetting
        EsPotIrr = EsPot

    # Assign minimum value (mulches and partial wetting don't combine)
    EsPot = min(EsPotIrr, EsPotMul)

    ## Surface evaporation ##
    # Initialise actual evaporation counter
    EsAct = 0
    # Evaporate surface storage
    if NewCond_SurfaceStorage > 0:
        if NewCond_SurfaceStorage > EsPot:
            # All potential soil evaporation can be supplied by surface storage
            EsAct = EsPot
            # Update surface storage
            NewCond_SurfaceStorage = NewCond_SurfaceStorage - EsAct
        else:
            # Surface storage is not sufficient to meet all potential soil
            # evaporation
            EsAct = NewCond_SurfaceStorage
            # Update surface storage, evaporation layer depth, stage
            NewCond_SurfaceStorage = 0
            NewCond_Wsurf = Soil_REW
            NewCond_Wstage2 = 0
            NewCond_EvapZ = Soil_EvapZmin
            NewCond_Stage2 = False

    ## stage 1 evaporation ##
    # Determine total water to be extracted
    ToExtract = EsPot - EsAct
    # Determine total water to be extracted in stage one (limited by surface
    # layer water storage)
    ExtractPotStg1 = min(ToExtract, NewCond_Wsurf)
    # Extract water
    if ExtractPotStg1 > 0:
        # Find soil compartments covered by evaporation layer
        comp_sto = np.sum(prof.dzsum < Soil_EvapZmin) + 1
        comp = -1
        # prof = Soil_Profile
        while (ExtractPotStg1 > 0) and (comp < comp_sto):
            # Increment compartment counter
            comp = comp + 1
            # Specify layer number
            # Determine proportion of compartment in evaporation layer
            if prof.dzsum[comp] > Soil_EvapZmin:
                factor = 1 - ((prof.dzsum[comp] - Soil_EvapZmin) / prof.dz[comp])
            else:
                factor = 1

            # Water storage (mm) at air dry
            Wdry = 1000 * prof.th_dry[comp] * prof.dz[comp]
            # Available water (mm)
            W = 1000 * NewCond_th[comp] * prof.dz[comp]
            # Water available in compartment for extraction (mm)
            AvW = (W - Wdry) * factor
            if AvW < 0:
                AvW = 0

            if AvW >= ExtractPotStg1:
                # Update actual evaporation
                EsAct = EsAct + ExtractPotStg1
                # Update depth of water in current compartment
                W = W - ExtractPotStg1
                # Update total water to be extracted
                ToExtract = ToExtract - ExtractPotStg1
                # Update water to be extracted from surface layer (stage 1)
                ExtractPotStg1 = 0
            else:
                # Update actual evaporation
                EsAct = EsAct + AvW
                # Update water to be extracted from surface layer (stage 1)
                ExtractPotStg1 = ExtractPotStg1 - AvW
                # Update total water to be extracted
                ToExtract = ToExtract - AvW
                # Update depth of water in current compartment
                W = W - AvW

            # Update water content
            NewCond_th[comp] = W / (1000 * prof.dz[comp])

        # Update surface evaporation layer water balance
        NewCond_Wsurf = NewCond_Wsurf - EsAct
        if (NewCond_Wsurf < 0) or (ExtractPotStg1 > 0.0001):
            NewCond_Wsurf = 0

        # If surface storage completely depleted, prepare stage 2
        if NewCond_Wsurf < 0.0001:
            # Get water contents (mm)
            Wevap_Sat, Wevap_Fc, Wevap_Wp, Wevap_Dry, Wevap_Act = evap_layer_water_content(
                NewCond_th,
                NewCond_EvapZ,
                prof,
            )
            # Proportional water storage for start of stage two evaporation
            NewCond_Wstage2 = round(
                (Wevap_Act - (Wevap_Fc - Soil_REW)) / (Wevap_Sat - (Wevap_Fc - Soil_REW)), 2
            )
            if NewCond_Wstage2 < 0:
                NewCond_Wstage2 = 0

    ## stage 2 evaporation ##
    # Extract water
    if ToExtract > 0:
        # Start stage 2
        NewCond_Stage2 = True
        # Get sub-daily evaporative demand
        Edt = ToExtract / ClockStruct_EvapTimeSteps
        # Loop sub-daily steps
        for jj in range(int(ClockStruct_EvapTimeSteps)):
            # Get current water storage (mm)
            Wevap_Sat, Wevap_Fc, Wevap_Wp, Wevap_Dry, Wevap_Act = evap_layer_water_content(
                NewCond_th,
                NewCond_EvapZ,
                prof,
            )
            # Get water storage (mm) at start of stage 2 evaporation
            Wupper = NewCond_Wstage2 * (Wevap_Sat - (Wevap_Fc - Soil_REW)) + (Wevap_Fc - Soil_REW)
            # Get water storage (mm) when there is no evaporation
            Wlower = Wevap_Dry
            # Get relative depletion of evaporation storage in stage 2
            Wrel = (Wevap_Act - Wlower) / (Wupper - Wlower)
            # Check if need to expand evaporation layer
            if Soil_EvapZmax > Soil_EvapZmin:
                Wcheck = Soil_fWrelExp * (
                    (Soil_EvapZmax - NewCond_EvapZ) / (Soil_EvapZmax - Soil_EvapZmin)
                )
                while (Wrel < Wcheck) and (NewCond_EvapZ < Soil_EvapZmax):
                    # Expand evaporation layer by 1 mm
                    NewCond_EvapZ = NewCond_EvapZ + 0.001
                    # Update water storage (mm) in evaporation layer
                    Wevap_Sat, Wevap_Fc, Wevap_Wp, Wevap_Dry, Wevap_Act = evap_layer_water_content(
                        NewCond_th,
                        NewCond_EvapZ,
                        prof,
                    )
                    Wupper = NewCond_Wstage2 * (Wevap_Sat - (Wevap_Fc - Soil_REW)) + (
                        Wevap_Fc - Soil_REW
                    )
                    Wlower = Wevap_Dry
                    # Update relative depletion of evaporation storage
                    Wrel = (Wevap_Act - Wlower) / (Wupper - Wlower)
                    Wcheck = Soil_fWrelExp * (
                        (Soil_EvapZmax - NewCond_EvapZ) / (Soil_EvapZmax - Soil_EvapZmin)
                    )

            # Get stage 2 evaporation reduction coefficient
            Kr = (np.exp(Soil_fevap * Wrel) - 1) / (np.exp(Soil_fevap) - 1)
            if Kr > 1:
                Kr = 1

            # Get water to extract (mm)
            ToExtractStg2 = Kr * Edt

            # Extract water from compartments
            comp_sto = np.sum(prof.dzsum < NewCond_EvapZ) + 1
            comp = -1
            # prof = Soil_Profile
            while (ToExtractStg2 > 0) and (comp < comp_sto):
                # Increment compartment counter
                comp = comp + 1
                # Specify layer number
                # Determine proportion of compartment in evaporation layer
                if prof.dzsum[comp] > NewCond_EvapZ:
                    factor = 1 - ((prof.dzsum[comp] - NewCond_EvapZ) / prof.dz[comp])
                else:
                    factor = 1

                # Water storage (mm) at air dry
                Wdry = 1000 * prof.th_dry[comp] * prof.dz[comp]
                # Available water (mm)
                W = 1000 * NewCond_th[comp] * prof.dz[comp]
                # Water available in compartment for extraction (mm)
                AvW = (W - Wdry) * factor
                if AvW >= ToExtractStg2:
                    # Update actual evaporation
                    EsAct = EsAct + ToExtractStg2
                    # Update depth of water in current compartment
                    W = W - ToExtractStg2
                    # Update total water to be extracted
                    ToExtract = ToExtract - ToExtractStg2
                    # Update water to be extracted from surface layer (stage 1)
                    ToExtractStg2 = 0
                else:
                    # Update actual evaporation
                    EsAct = EsAct + AvW
                    # Update depth of water in current compartment
                    W = W - AvW
                    # Update water to be extracted from surface layer (stage 1)
                    ToExtractStg2 = ToExtractStg2 - AvW
                    # Update total water to be extracted
                    ToExtract = ToExtract - AvW

                # Update water content
                NewCond_th[comp] = W / (1000 * prof.dz[comp])

    ## Store potential evaporation for irrigation calculations on next day ##
    NewCond_Epot = EsPot

    return (
        NewCond_Epot,
        NewCond_th,
        NewCond_Stage2,
        NewCond_Wstage2,
        NewCond_Wsurf,
        NewCond_SurfaceStorage,
        NewCond_EvapZ,
        EsAct,
        EsPot,
    )

aquacrop.solution.temperature_stress

temperature_stress(Crop, temp_max, temp_min)

Function to get irrigation depth for current day

Reference Manual: temperature stress (pg. 14)

Arguments:

Crop (Crop): Crop object containing Crop paramaters

temp_max (float): max tempatature on current day (celcius)

temp_min (float): min tempature on current day (celcius)

Returns:

Kst_PolH (float): heat stress coefficient for current day

Kst_PolC (float): cold stress coefficient for current day
Source code in aquacrop/solution/temperature_stress.py 
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
@cc.export("temperature_stress", (CropStructNT_type_sig,f8,f8))
def temperature_stress(
    Crop: "CropStructNT",
    temp_max: float,
    temp_min: float,
    ) -> Tuple[float,float]:
    # Function to calculate temperature stress coefficients
    """
    Function to get irrigation depth for current day

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=23" target="_blank">Reference Manual: temperature stress</a> (pg. 14)



    Arguments:

        Crop (Crop): Crop object containing Crop paramaters

        temp_max (float): max tempatature on current day (celcius)

        temp_min (float): min tempature on current day (celcius)


    Returns:

        Kst_PolH (float): heat stress coefficient for current day

        Kst_PolC (float): cold stress coefficient for current day







    """

    ## Calculate temperature stress coefficients affecting crop pollination ##
    # Get parameters for logistic curve
    KsPol_up = 1
    KsPol_lo = 0.001

    # Kst = Kst()

    # Calculate effects of heat stress on pollination
    if Crop.PolHeatStress == 0:
        # No heat stress effects on pollination
        Kst_PolH = 1
    elif Crop.PolHeatStress == 1:
        # Pollination affected by heat stress
        if temp_max <= Crop.Tmax_lo:
            Kst_PolH = 1
        elif temp_max >= Crop.Tmax_up:
            Kst_PolH = 0
        else:
            Trel = (temp_max - Crop.Tmax_lo) / (Crop.Tmax_up - Crop.Tmax_lo)
            Kst_PolH = (KsPol_up * KsPol_lo) / (
                KsPol_lo + (KsPol_up - KsPol_lo) * np.exp(-Crop.fshape_b * (1 - Trel))
            )

    # Calculate effects of cold stress on pollination
    if Crop.PolColdStress == 0:
        # No cold stress effects on pollination
        Kst_PolC = 1
    elif Crop.PolColdStress == 1:
        # Pollination affected by cold stress
        if temp_min >= Crop.Tmin_up:
            Kst_PolC = 1
        elif temp_min <= Crop.Tmin_lo:
            Kst_PolC = 0
        else:
            Trel = (Crop.Tmin_up - temp_min) / (Crop.Tmin_up - Crop.Tmin_lo)
            Kst_PolC = (KsPol_up * KsPol_lo) / (
                KsPol_lo + (KsPol_up - KsPol_lo) * np.exp(-Crop.fshape_b * (1 - Trel))
            )

    return (Kst_PolH,Kst_PolC)

aquacrop.solution.transpiration

transpiration(Soil_Profile, Soil_nComp, Soil_zTop, Crop, IrrMngt_IrrMethod, IrrMngt_NetIrrSMT, InitCond, et0, CO2, growing_season, gdd)

Function to calculate crop transpiration on current day

Reference Manual: transpiration equations (pg. 82-91)

Arguments:

Soil_Profile (SoilProfileNT): Soil profile params

Soil_nComp (int): number of soil components

Soil_zTop (float): depth of topsoil

Crop (Crop): Crop params

IrrMngt_IrrMethod (int): irrigation method

IrrMngt_NetIrrSMT (float): net irrigation soil-moisture target

InitCond (InitialCondition): InitCond object

et0 (float): reference evapotranspiration

CO2 (CO2): CO2

gdd (float): Growing Degree Days

growing_season (bool): is it currently within the growing season (True, Flase)

Returns:

TrAct (float): Actual Transpiration on current day

TrPot_NS (float): Potential Transpiration on current day with no water stress

TrPot0 (float): Potential Transpiration on current day

NewCond (InitialCondition): updated InitCond object

IrrNet (float): Net Irrigation (if required)
Source code in aquacrop/solution/transpiration.py 
 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
def transpiration(
    Soil_Profile: "SoilProfileNT",
    Soil_nComp: int,
    Soil_zTop: float,
    Crop: "CropStructNT",
    IrrMngt_IrrMethod: int,
    IrrMngt_NetIrrSMT: float,
    InitCond: "InitialCondition",
    et0: float,
    CO2: "CO2",
    growing_season: bool,
    gdd: float,
) -> Tuple[float,float,float,"InitialCondition",float]:

    """
    Function to calculate crop transpiration on current day

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=91" target="_blank">Reference Manual: transpiration equations</a> (pg. 82-91)



    Arguments:


        Soil_Profile (SoilProfileNT): Soil profile params

        Soil_nComp (int): number of soil components

        Soil_zTop (float): depth of topsoil

        Crop (Crop): Crop params

        IrrMngt_IrrMethod (int): irrigation method 

        IrrMngt_NetIrrSMT (float): net irrigation soil-moisture target

        InitCond (InitialCondition): InitCond object

        et0 (float): reference evapotranspiration

        CO2 (CO2): CO2

        gdd (float): Growing Degree Days

        growing_season (bool): is it currently within the growing season (True, Flase)

    Returns:


        TrAct (float): Actual Transpiration on current day

        TrPot_NS (float): Potential Transpiration on current day with no water stress

        TrPot0 (float): Potential Transpiration on current day

        NewCond (InitialCondition): updated InitCond object

        IrrNet (float): Net Irrigation (if required)







    """

    ## Store initial conditions ##
    NewCond = InitCond

    InitCond_th = InitCond.th

    prof = Soil_Profile

    ## Calculate transpiration (if in growing season) ##
    if growing_season == True:
        ## Calculate potential transpiration ##
        # 1. No prior water stress
        # Update ageing days counter
        DAPadj = NewCond.dap - NewCond.delayed_cds
        if DAPadj > Crop.MaxCanopyCD:
            NewCond.age_days_ns = DAPadj - Crop.MaxCanopyCD

        # Update crop coefficient for ageing of canopy
        if NewCond.age_days_ns > 5:
            Kcb_NS = Crop.Kcb - ((NewCond.age_days_ns - 5) * (Crop.fage / 100)) * NewCond.ccx_w_ns
        else:
            Kcb_NS = Crop.Kcb

        # Update crop coefficient for CO2 concentration
        CO2CurrentConc = CO2.current_concentration
        CO2RefConc = CO2.ref_concentration
        if CO2CurrentConc > CO2RefConc:
            Kcb_NS = Kcb_NS * (1 - 0.05 * ((CO2CurrentConc - CO2RefConc) / (550 - CO2RefConc)))

        # Determine potential transpiration rate (no water stress)
        TrPot_NS = Kcb_NS * (NewCond.canopy_cover_adj_ns) * et0

        # Correct potential transpiration for dying green canopy effects
        if NewCond.canopy_cover_ns < NewCond.ccx_w_ns:
            if (NewCond.ccx_w_ns > 0.001) and (NewCond.canopy_cover_ns > 0.001):
                TrPot_NS = TrPot_NS * ((NewCond.canopy_cover_ns / NewCond.ccx_w_ns) ** Crop.a_Tr)

        # 2. Potential prior water stress and/or delayed development
        # Update ageing days counter
        DAPadj = NewCond.dap - NewCond.delayed_cds
        if DAPadj > Crop.MaxCanopyCD:
            NewCond.age_days = DAPadj - Crop.MaxCanopyCD

        # Update crop coefficient for ageing of canopy
        if NewCond.age_days > 5:
            Kcb = Crop.Kcb - ((NewCond.age_days - 5) * (Crop.fage / 100)) * NewCond.ccx_w
        else:
            Kcb = Crop.Kcb

        # Update crop coefficient for CO2 concentration
        if CO2CurrentConc > CO2RefConc:
            Kcb = Kcb * (1 - 0.05 * ((CO2CurrentConc - CO2RefConc) / (550 - CO2RefConc)))

        # Determine potential transpiration rate
        TrPot0 = Kcb * (NewCond.canopy_cover_adj) * et0
        # Correct potential transpiration for dying green canopy effects
        if NewCond.canopy_cover < NewCond.ccx_w:
            if (NewCond.ccx_w > 0.001) and (NewCond.canopy_cover > 0.001):
                TrPot0 = TrPot0 * ((NewCond.canopy_cover / NewCond.ccx_w) ** Crop.a_Tr)

        # 3. Adjust potential transpiration for cold stress effects
        # Check if cold stress occurs on current day
        if Crop.TrColdStress == 0:
            # Cold temperature stress does not affect transpiration
            KsCold = 1
        elif Crop.TrColdStress == 1:
            # Transpiration can be affected by cold temperature stress
            if gdd >= Crop.GDD_up:
                # No cold temperature stress
                KsCold = 1
            elif gdd <= Crop.GDD_lo:
                # Transpiration fully inhibited by cold temperature stress
                KsCold = 0
            else:
                # Transpiration partially inhibited by cold temperature stress
                # Get parameters for logistic curve
                KsTr_up = 1
                KsTr_lo = 0.02
                fshapeb = (-1) * (
                    np.log(((KsTr_lo * KsTr_up) - 0.98 * KsTr_lo) / (0.98 * (KsTr_up - KsTr_lo)))
                )
                # Calculate cold stress level
                GDDrel = (gdd - Crop.GDD_lo) / (Crop.GDD_up - Crop.GDD_lo)
                KsCold = (KsTr_up * KsTr_lo) / (
                    KsTr_lo + (KsTr_up - KsTr_lo) * np.exp(-fshapeb * GDDrel)
                )
                KsCold = KsCold - KsTr_lo * (1 - GDDrel)

        # Correct potential transpiration rate (mm/day)
        TrPot0 = TrPot0 * KsCold
        TrPot_NS = TrPot_NS * KsCold

        # print(TrPot0,NewCond.dap)

        ## Calculate surface layer transpiration ##
        if (NewCond.surface_storage > 0) and (NewCond.day_submerged < Crop.LagAer):

            # Update submergence days counter
            NewCond.day_submerged = NewCond.day_submerged + 1
            # Update anerobic conditions counter for each compartment
            for ii in range(int(Soil_nComp)):
                # Increment aeration days counter for compartment ii
                NewCond.aer_days_comp[ii] = NewCond.aer_days_comp[ii] + 1
                if NewCond.aer_days_comp[ii] > Crop.LagAer:
                    NewCond.aer_days_comp[ii] = Crop.LagAer

            # Reduce actual transpiration that is possible to account for
            # aeration stress due to extended submergence
            fSub = 1 - (NewCond.day_submerged / Crop.LagAer)
            if NewCond.surface_storage > (fSub * TrPot0):
                # Transpiration occurs from surface storage
                NewCond.surface_storage = NewCond.surface_storage - (fSub * TrPot0)
                TrAct0 = fSub * TrPot0
            else:
                # No transpiration from surface storage
                TrAct0 = 0

            if TrAct0 < (fSub * TrPot0):
                # More water can be extracted from soil profile for transpiration
                TrPot = (fSub * TrPot0) - TrAct0
                # print('now')

            else:
                # No more transpiration possible on current day
                TrPot = 0
                # print('here')

        else:

            # No surface transpiration occurs
            TrPot = TrPot0
            TrAct0 = 0

        # print(TrPot,NewCond.dap)

        ## Update potential root zone transpiration for water stress ##
        # Determine root zone and top soil depletion, and root zone water
        # content

        taw = TAW()
        water_root_depletion = Dr()
        thRZ = RootZoneWater()
        (
            _,
            water_root_depletion.Zt,
            water_root_depletion.Rz,
            taw.Zt,
            taw.Rz,
            thRZ.Act,
            thRZ.S,
            thRZ.FC,
            thRZ.WP,
            thRZ.Dry,
            thRZ.Aer,
        ) = root_zone_water(
            prof,
            float(NewCond.z_root),
            NewCond.th,
            Soil_zTop,
            float(Crop.Zmin),
            Crop.Aer,
        )

        class_args = {key:value for key, value in thRZ.__dict__.items() if not key.startswith('__') and not callable(key)}
        thRZ = thRZNT(**class_args)

        # _,water_root_depletion,taw,thRZ = root_zone_water(Soil_Profile,float(NewCond.z_root),NewCond.th,Soil_zTop,float(Crop.Zmin),Crop.Aer)
        # Check whether to use root zone or top soil depletions for calculating
        # water stress
        if (water_root_depletion.Rz / taw.Rz) <= (water_root_depletion.Zt / taw.Zt):
            # Root zone is wetter than top soil, so use root zone value
            water_root_depletion = water_root_depletion.Rz
            taw = taw.Rz
        else:
            # Top soil is wetter than root zone, so use top soil values
            water_root_depletion = water_root_depletion.Zt
            taw = taw.Zt

        # Calculate water stress coefficients
        beta = True
        water_stress_coef = Ksw()
        water_stress_coef.exp, water_stress_coef.sto, water_stress_coef.sen, water_stress_coef.pol, water_stress_coef.sto_lin = water_stress(
            Crop.p_up,
            Crop.p_lo,
            Crop.ETadj,
            Crop.beta,
            Crop.fshape_w,
            NewCond.t_early_sen,
            water_root_depletion,
            taw,
            et0,
            beta,
        )
        # water_stress_coef = water_stress(Crop, NewCond, water_root_depletion, taw, et0, beta)

        # Calculate aeration stress coefficients
        Ksa_Aer, NewCond.aer_days = aeration_stress(NewCond.aer_days, Crop.LagAer, thRZ)
        # Maximum stress effect
        Ks = min(water_stress_coef.sto_lin, Ksa_Aer)
        # Update potential transpiration in root zone
        if IrrMngt_IrrMethod != 4:
            # No adjustment to TrPot for water stress when in net irrigation mode
            TrPot = TrPot * Ks

        ## Determine compartments covered by root zone ##
        # Compartments covered by the root zone
        rootdepth = round(max(float(NewCond.z_root), float(Crop.Zmin)), 2)
        comp_sto = min(np.sum(Soil_Profile.dzsum < rootdepth) + 1, int(Soil_nComp))
        RootFact = np.zeros(int(Soil_nComp))
        # Determine fraction of each compartment covered by root zone
        for ii in range(comp_sto):
            if Soil_Profile.dzsum[ii] > rootdepth:
                RootFact[ii] = 1 - ((Soil_Profile.dzsum[ii] - rootdepth) / Soil_Profile.dz[ii])
            else:
                RootFact[ii] = 1

        ## Determine maximum sink term for each compartment ##
        SxComp = np.zeros(int(Soil_nComp))
        if IrrMngt_IrrMethod == 4:
            # Net irrigation mode
            for ii in range(comp_sto):
                SxComp[ii] = (Crop.SxTop + Crop.SxBot) / 2

        else:
            # Maximum sink term declines linearly with depth
            SxCompBot = Crop.SxTop
            for ii in range(comp_sto):
                SxCompTop = SxCompBot
                if Soil_Profile.dzsum[ii] <= rootdepth:
                    SxCompBot = Crop.SxBot * NewCond.r_cor + (
                        (Crop.SxTop - Crop.SxBot * NewCond.r_cor)
                        * ((rootdepth - Soil_Profile.dzsum[ii]) / rootdepth)
                    )
                else:
                    SxCompBot = Crop.SxBot * NewCond.r_cor

                SxComp[ii] = (SxCompTop + SxCompBot) / 2

        # print(TrPot,NewCond.dap)
        ## Extract water ##
        ToExtract = TrPot
        comp = -1
        TrAct = 0
        while (ToExtract > 0) and (comp < comp_sto - 1):
            # Increment compartment
            comp = comp + 1
            # Specify layer number

            # Determine taw (m3/m3) for compartment
            thTAW = prof.th_fc[comp] - prof.th_wp[comp]
            if Crop.ETadj == 1:
                # Adjust stomatal stress threshold for et0 on current day
                p_up_sto = Crop.p_up[1] + (0.04 * (5 - et0)) * (np.log10(10 - 9 * Crop.p_up[1]))

            # Determine critical water content at which stomatal closure will
            # occur in compartment
            thCrit = prof.th_fc[comp] - (thTAW * p_up_sto)

            # Check for soil water stress
            if NewCond.th[comp] >= thCrit:
                # No water stress effects on transpiration
                KsComp = 1
            elif NewCond.th[comp] > prof.th_wp[comp]:
                # Transpiration from compartment is affected by water stress
                Wrel = (prof.th_fc[comp] - NewCond.th[comp]) / (prof.th_fc[comp] - prof.th_wp[comp])
                pRel = (Wrel - Crop.p_up[1]) / (Crop.p_lo[1] - Crop.p_up[1])
                if pRel <= 0:
                    KsComp = 1
                elif pRel >= 1:
                    KsComp = 0
                else:
                    KsComp = 1 - (
                        (np.exp(pRel * Crop.fshape_w[1]) - 1) / (np.exp(Crop.fshape_w[1]) - 1)
                    )

                if KsComp > 1:
                    KsComp = 1
                elif KsComp < 0:
                    KsComp = 0

            else:
                # No transpiration is possible from compartment as water
                # content does not exceed wilting point
                KsComp = 0

            # Adjust compartment stress factor for aeration stress
            if NewCond.day_submerged >= Crop.LagAer:
                # Full aeration stress - no transpiration possible from
                # compartment
                AerComp = 0
            elif NewCond.th[comp] > (prof.th_s[comp] - (Crop.Aer / 100)):
                # Increment aeration stress days counter
                NewCond.aer_days_comp[comp] = NewCond.aer_days_comp[comp] + 1
                if NewCond.aer_days_comp[comp] >= Crop.LagAer:
                    NewCond.aer_days_comp[comp] = Crop.LagAer
                    fAer = 0
                else:
                    fAer = 1

                # Calculate aeration stress factor
                AerComp = (prof.th_s[comp] - NewCond.th[comp]) / (
                    prof.th_s[comp] - (prof.th_s[comp] - (Crop.Aer / 100))
                )
                if AerComp < 0:
                    AerComp = 0

                AerComp = (fAer + (NewCond.aer_days_comp[comp] - 1) * AerComp) / (
                    fAer + NewCond.aer_days_comp[comp] - 1
                )
            else:
                # No aeration stress as number of submerged days does not
                # exceed threshold for initiation of aeration stress
                AerComp = 1
                NewCond.aer_days_comp[comp] = 0

            # Extract water
            ThToExtract = (ToExtract / 1000) / Soil_Profile.dz[comp]
            if IrrMngt_IrrMethod == 4:
                # Don't reduce compartment sink for stomatal water stress if in
                # net irrigation mode. Stress only occurs due to deficient
                # aeration conditions
                Sink = AerComp * SxComp[comp] * RootFact[comp]
            else:
                # Reduce compartment sink for greatest of stomatal and aeration
                # stress
                if KsComp == AerComp:
                    Sink = KsComp * SxComp[comp] * RootFact[comp]
                else:
                    Sink = min(KsComp, AerComp) * SxComp[comp] * RootFact[comp]

            # Limit extraction to demand
            if ThToExtract < Sink:
                Sink = ThToExtract

            # Limit extraction to avoid compartment water content dropping
            # below air dry
            if (InitCond_th[comp] - Sink) < prof.th_dry[comp]:
                Sink = InitCond_th[comp] - prof.th_dry[comp]
                if Sink < 0:
                    Sink = 0

            # Update water content in compartment
            NewCond.th[comp] = InitCond_th[comp] - Sink
            # Update amount of water to extract
            ToExtract = ToExtract - (Sink * 1000 * prof.dz[comp])
            # Update actual transpiration
            TrAct = TrAct + (Sink * 1000 * prof.dz[comp])

        ## Add net irrigation water requirement (if this mode is specified) ##
        if (IrrMngt_IrrMethod == 4) and (TrPot > 0):
            # Initialise net irrigation counter
            IrrNet = 0
            # Get root zone water content

            taw = TAW()
            water_root_depletion = Dr()
            thRZ = RootZoneWater()
            (
                _,
                water_root_depletion.Zt,
                water_root_depletion.Rz,
                taw.Zt,
                taw.Rz,
                thRZ.Act,
                thRZ.S,
                thRZ.FC,
                thRZ.WP,
                thRZ.Dry,
                thRZ.Aer,
            ) = root_zone_water(
                prof,
                float(NewCond.z_root),
                NewCond.th,
                Soil_zTop,
                float(Crop.Zmin),
                Crop.Aer,
            )

            # _,_Dr,_TAW,thRZ = root_zone_water(Soil_Profile,float(NewCond.z_root),NewCond.th,Soil_zTop,float(Crop.Zmin),Crop.Aer)
            NewCond.depletion = water_root_depletion.Rz
            NewCond.taw = taw.Rz
            # Determine critical water content for net irrigation
            thCrit = thRZ.WP + ((IrrMngt_NetIrrSMT / 100) * (thRZ.FC - thRZ.WP))
            # Check if root zone water content is below net irrigation trigger
            if thRZ.Act < thCrit:
                # Initialise layer counter
                prelayer = 0
                for ii in range(comp_sto):
                    # Get soil layer
                    layeri = Soil_Profile.Layer[ii]
                    if layeri > prelayer:
                        # If in new layer, update critical water content for
                        # net irrigation
                        thCrit = prof.th_wp[ii] + (
                            (IrrMngt_NetIrrSMT / 100) * (prof.th_fc[ii] - prof.th_wp[ii])
                        )
                        # Update layer counter
                        prelayer = layeri

                    # Determine necessary change in water content in
                    # compartments to reach critical water content
                    dWC = RootFact[ii] * (thCrit - NewCond.th[ii]) * 1000 * prof.dz[ii]
                    # Update water content
                    NewCond.th[ii] = NewCond.th[ii] + (dWC / (1000 * prof.dz[ii]))
                    # Update net irrigation counter
                    IrrNet = IrrNet + dWC

            # Update net irrigation counter for the growing season
            NewCond.irr_net_cum = NewCond.irr_net_cum + IrrNet
        elif (IrrMngt_IrrMethod == 4) and (TrPot <= 0):
            # No net irrigation as potential transpiration is zero
            IrrNet = 0
        else:
            # No net irrigation as not in net irrigation mode
            IrrNet = 0
            NewCond.irr_net_cum = 0

        ## Add any surface transpiration to root zone total ##
        TrAct = TrAct + TrAct0

        ## Feedback with canopy cover development ##
        # If actual transpiration is zero then no canopy cover growth can occur
        if ((NewCond.canopy_cover - NewCond.cc_prev) > 0.005) and (TrAct == 0):
            NewCond.canopy_cover = NewCond.cc_prev

        ## Update transpiration ratio ##
        if TrPot0 > 0:
            if TrAct < TrPot0:
                NewCond.tr_ratio = TrAct / TrPot0
            else:
                NewCond.tr_ratio = 1

        else:
            NewCond.tr_ratio = 1

        if NewCond.tr_ratio < 0:
            NewCond.tr_ratio = 0
        elif NewCond.tr_ratio > 1:
            NewCond.tr_ratio = 1

    else:
        # No transpiration if not in growing season
        TrAct = 0
        TrPot0 = 0
        TrPot_NS = 0
        # No irrigation if not in growing season
        IrrNet = 0
        NewCond.irr_net_cum = 0

    ## Store potential transpiration for irrigation calculations on next day ##
    NewCond.t_pot = TrPot0

    return TrAct, TrPot_NS, TrPot0, NewCond, IrrNet

aquacrop.solution.update_CCx_CDC

update_CCx_CDC(cc_prev, CDC, CCx, dt)

Function to update CCx and CDC parameter valyes for rewatering in late season of an early declining canopy

Reference Manual: canopy_cover stress response (pg. 27-33)

Arguments:

cc_prev (float): Canopy Cover at previous timestep.

CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

CCx (float): Maximum canopy cover (fraction of soil cover)

dt (float): Time delta of canopy growth (1 calander day or ... gdd)

Returns:

CCxAdj (float): updated CCxAdj

CDCadj (float): updated CDCadj
Source code in aquacrop/solution/update_CCx_CDC.py 
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
@cc.export("update_CCx_CDC", "(f8,f8,f8,f8)")
def update_CCx_CDC(
    cc_prev: float,
    CDC: float,
    CCx: float,
    dt: float,
    ) -> Tuple[float,float]:
    """
    Function to update CCx and CDC parameter valyes for rewatering in late season of an early declining canopy

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=36" target="_blank">Reference Manual: canopy_cover stress response</a> (pg. 27-33)


    Arguments:


        cc_prev (float): Canopy Cover at previous timestep.

        CDC (float): Canopy decline coefficient (fraction per gdd/calendar day)

        CCx (float): Maximum canopy cover (fraction of soil cover)

        dt (float): Time delta of canopy growth (1 calander day or ... gdd)


    Returns:

        CCxAdj (float): updated CCxAdj

        CDCadj (float): updated CDCadj





    """

    ## Get adjusted CCx ##
    CCXadj = cc_prev / (1 - 0.05 * (np.exp(dt * ((CDC * 3.33) / (CCx + 2.29))) - 1))

    ## Get adjusted CDC ##
    CDCadj = CDC * ((CCXadj + 2.29) / (CCx + 2.29))

    return CCXadj, CDCadj

aquacrop.solution.water_stress

water_stress(Crop_p_up, Crop_p_lo, Crop_ETadj, Crop_beta, Crop_fshape_w, InitCond_tEarlySen, Dr, taw, et0, beta)

Function to calculate water stress coefficients

Reference Manual: water stress equations (pg. 9-13)

Arguments:

Crop_p_up (ndarray): water stress thresholds for start of water stress

Crop_p_lo (ndarray): water stress thresholds for maximum water stress

Crop_ETadj (float):

Crop_beta (float):

Crop_fshape_w (ndarray): shape factors for water stress

InitCond_tEarlySen (float): days in early senesence

Dr (Dr): rootzone depletion

taw (TAW): root zone total available water

et0 (float): Reference Evapotranspiration

beta (float): Adjust senescence threshold if early sensescence is triggered

Returns:

Ksw (Ksw): Ksw object containint water stress coefficients
Source code in aquacrop/solution/water_stress.py 
 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
@njit
@cc.export("water_stress", "(f8[:],f8[:],f8,f8,f8[:],f8,f8,f8,f8,f8)")
def water_stress(
    Crop_p_up: "ndarray",
    Crop_p_lo: "ndarray",
    Crop_ETadj: float,
    Crop_beta: float,
    Crop_fshape_w: "ndarray",
    InitCond_tEarlySen: float,
    Dr: float,
    taw: float,
    et0: float,
    beta: float,
) -> Tuple[float, float, float, float, float]:
    """
    Function to calculate water stress coefficients

    <a href="https://www.fao.org/3/BR248E/br248e.pdf#page=18" target="_blank">Reference Manual: water stress equations</a> (pg. 9-13)


    Arguments:


        Crop_p_up (ndarray): water stress thresholds for start of water stress

        Crop_p_lo (ndarray): water stress thresholds for maximum water stress

        Crop_ETadj (float): 

        Crop_beta (float): 

        Crop_fshape_w (ndarray): shape factors for water stress

        InitCond_tEarlySen (float): days in early senesence

        Dr (Dr): rootzone depletion

        taw (TAW): root zone total available water

        et0 (float): Reference Evapotranspiration

        beta (float): Adjust senescence threshold if early sensescence is triggered


    Returns:

        Ksw (Ksw): Ksw object containint water stress coefficients


    """

    ## Calculate relative root zone water depletion for each stress type ##
    # Number of stress variables
    nstress = len(Crop_p_up)

    # Store stress thresholds
    p_up = np.ones(nstress) * Crop_p_up
    p_lo = np.ones(nstress) * Crop_p_lo
    if Crop_ETadj == 1:
        # Adjust stress thresholds for et0 on currentbeta day (don't do this for
        # pollination water stress coefficient)

        for ii in range(3):
            p_up[ii] = p_up[ii] + (0.04 * (5 - et0)) * (np.log10(10 - 9 * p_up[ii]))
            p_lo[ii] = p_lo[ii] + (0.04 * (5 - et0)) * (np.log10(10 - 9 * p_lo[ii]))

    # Adjust senescence threshold if early sensescence is triggered
    if (beta == True) and (InitCond_tEarlySen > 0):
        p_up[2] = p_up[2] * (1 - Crop_beta / 100)

    # Limit values
    p_up = np.maximum(p_up, np.zeros(4))
    p_lo = np.maximum(p_lo, np.zeros(4))
    p_up = np.minimum(p_up, np.ones(4))
    p_lo = np.minimum(p_lo, np.ones(4))

    # Calculate relative depletion
    Drel = np.zeros(nstress)
    for ii in range(nstress):
        if Dr <= (p_up[ii] * taw):
            # No water stress
            Drel[ii] = 0
        elif (Dr > (p_up[ii] * taw)) and (Dr < (p_lo[ii] * taw)):
            # Partial water stress
            Drel[ii] = 1 - ((p_lo[ii] - (Dr / taw)) / (p_lo[ii] - p_up[ii]))
        elif Dr >= (p_lo[ii] * taw):
            # Full water stress
            Drel[ii] = 1

    ## Calculate root zone water stress coefficients ##
    Ks = np.ones(3)
    for ii in range(3):
        Ks[ii] = 1 - ((np.exp(Drel[ii] * Crop_fshape_w[ii]) - 1) / (np.exp(Crop_fshape_w[ii]) - 1))

    # Ksw = Ksw()

    # Water stress coefficient for leaf expansion
    Ksw_Exp = Ks[0]
    # Water stress coefficient for stomatal closure
    Ksw_Sto = Ks[1]
    # Water stress coefficient for senescence
    Ksw_Sen = Ks[2]
    # Water stress coefficient for pollination failure
    Ksw_Pol = 1 - Drel[3]
    # Mean water stress coefficient for stomatal closure
    Ksw_StoLin = 1 - Drel[1]

    return Ksw_Exp, Ksw_Sto, Ksw_Sen, Ksw_Pol, Ksw_StoLin