Opened 5 years ago
Last modified 2 years ago
#450 accepted defect
Negative soil moisture values
| Reported by: | maignan | Owned by: | aducharne |
|---|---|---|---|
| Priority: | major | Milestone: | Not scheduled yet |
| Component: | Physical processes | Version: | trunc |
| Keywords: | Cc: |
Description
- This summer during the preparation of the TRENDY simulations with the TRUNK, we had negative soil moisture values. We didn't go further as we soon discovered there was a problem with the time step of the forcing files.
- Albert is now debugging (in MICT) and has again negative soil moisture values.
Could we devise a test ensuring that soil moisture stay positive?
Change History (6)
comment:1 Changed 5 years ago by ajornet
comment:2 follow-up: ↓ 6 Changed 5 years ago by aducharne
Hello,
I have a couple of questions to better apprehend the pb:
- does negative mcl affect the first layer only?
- do you use dynroot ? soil freezing ? can't the pb come from these options?
- does the pb happen with the trunk?
Then, in case it is needed to prevent mc to be negative, it could be introduced in hydrol_soil_smooth_under_mcr. Thre used to be a "smoothing" to move the deficit under mcr to soil layers with more moisture. This feature has been removed when we cleaned hydrol 2 years ago (ticket #191, rev 3402).
Agnès
comment:3 Changed 5 years ago by ajornet
Yes, it does only affect the first layer (nstm=1 and nslm=1).
DYNROOT=y, when changed to false, mcl can be negative
FREEZE=y, when changed to false, mcl can be negative
Yes, it does also happen in the trunk.
Some extra info:
mcl is negative after the call to hydrol_soil_tridiag
!! 2.7 Solves diffusion equations, but only in grid-cells where resolv is true, i.e. everywhere (cf 2.2) !! The result is an updated mcl profile CALL hydrol_soil_tridiag(kjpindex,jst)
Some more examples, printed values inside hydrol_soil_tridiag:
hydrol_soil_tridiag:: 1st i,j,v= 52 1 1 hydrol_soil_tridiag:: 1st, rhs= -7.264605390340498E-002 hydrol_soil_tridiag:: 1st, bet= 75.9814400817341 hydrol_soil_tridiag:: 1st, mcl= -9.561026196036661E-004 hydrol_soil_tridiag:: 1st, mcl= 0.296287108532727 hydrol_soil_tridiag:: 1st i,j,v= 53 1 1 hydrol_soil_tridiag:: 1st, rhs= -0.120159074728964 hydrol_soil_tridiag:: 1st, bet= 2.61584278570163 hydrol_soil_tridiag:: 1st, mcl= -4.593512858867570E-002 hydrol_soil_tridiag:: 1st, mcl= 0.320481824185928
Some values when rhs is negative:
hydrol_soil:: 1st ji= 52 hydrol_soil:: 1st fp= 0.733137830250000 hydrol_soil:: 1st mcl= 0.296803248260190 hydrol_soil:: 1st gp= 0.244379276750000 hydrol_soil:: 1st fp*mcl= 0.217219287880730 hydrol_soil:: 1st gp*mcl= 7.240642929357657E-002 hydrol_soil:: 1st flux_top= 0.520784622231542 hydrol_soil:: 1st b(1)= -7.60256243042550 hydrol_soil:: 1st b(2)= -7.60256243042550 hydrol_soil:: 1st b(1)+b(2)= -0.158386717300531 hydrol_soil:: 1st rootsink= 0.000000000000000E+000 hydrol_soil:: 1st rhs= -7.264605390340498E-002 hydrol_soil:: 1st ji= 53 hydrol_soil:: 1st fp= 0.733137830250000 hydrol_soil:: 1st mcl= 0.317498002996160 hydrol_soil:: 1st gp= 0.244379276750000 hydrol_soil:: 1st fp*mcl= 0.234957349218234 hydrol_soil:: 1st gp*mcl= 7.831911640607786E-002 hydrol_soil:: 1st flux_top= 0.432751703097998 hydrol_soil:: 1st b(1)= -2.176646833432279E-003 hydrol_soil:: 1st b(2)= -2.176646833432279E-003 hydrol_soil:: 1st b(1)+b(2)= -4.534680902983915E-005 hydrol_soil:: 1st rootsink= 0.000000000000000E+000 hydrol_soil:: 1st rhs= -0.120159074728964
I need to take a deeper look at the revision [3402]. It's not obvious to me where is the modification/
comment:4 Changed 5 years ago by aducharne
- Milestone changed from ORCHIDEE 3.0 to ORCHIDEE 4.0
- Owner changed from somebody to aducharne
- Status changed from new to accepted
comment:5 Changed 3 years ago by luyssaert
- Milestone changed from ORCHIDEE 4.0 to Not scheduled yet
comment:6 in reply to: ↑ 2 Changed 2 years ago by aducharne
Replying to aducharne:
Hello,
I have a couple of questions to better apprehend the pb:
- does negative mcl affect the first layer only?
- do you use dynroot ? soil freezing ? can't the pb come from these options?
- does the pb happen with the trunk?
Then, in case it is needed to prevent mc to be negative, it could be introduced in hydrol_soil_smooth_under_mcr. Thre used to be a "smoothing" to move the deficit under mcr to soil layers with more moisture. This feature has been removed when we cleaned hydrol 2 years ago (ticket #191, rev 3402).
Agnès
Correction on July 20, 2021: hydrol_smooth_under_mcr is still active, with an active smoothing negative (mc-mcr) towards layers where mc> mcr; it is hydrol_smooth_over_mcs which does not include smoothing anymore.
For more information:
In the subroutine hydrol_soil_tridiag mcl turns negative for the first time:
DO ji = 1, kjpindex IF (resolv(ji)) THEN bet(ji) = tmat(ji,1,2) mcl(ji,1,ins) = rhs(ji,1)/bet(ji) ENDIF ENDDO IF (ANY(mcl<0)) THEN ! introduced for debugging purposes DO jsl = 1, nslm DO ji=1, kjpindex IF (mcl(ji,jsl,ins) < 0) THEN WRITE(numout,*) "hydrol_soil_tridiag:: (1st) ji,jsl,jst=",ji,jsl,ins WRITE(numout,*) "hydrol_soil_tridiag:: (1st) bet(ji)=",bet(ji) WRITE(numout,*) "hydrol_soil_tridiag:: (1st) rhs(ji,1)=",rhs(ji,1) WRITE(numout,*) "hydrol_soil_tridiag:: (1st) mcl=",mcl(ji,jsl,ins) ENDIF ENDDO ENDDO CALL ipslerr_p(3,'hydrol_soil_tridiag','mcl has negative values','1st','') ENDIFWith the follwing output:
This mcl negatives values are then moved to mc.
Here is how rhs is calculated.
!- First layer DO ji = 1, kjpindex tmat(ji,1,1) = zero tmat(ji,1,2) = f(ji,1) tmat(ji,1,3) = g1(ji,1) rhs(ji,1) = fp(ji,1) * mcl(ji,1,jst) + gp(ji,1)*mcl(ji,2,jst) & & - flux_top(ji) - (b(ji,1)+b(ji,2))/deux *(dt_sechiba/one_day) - rootsink(ji,1,jst) ENDDO