Opened 3 years ago

Last modified 14 months 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 (5)

comment:1 Changed 3 years ago by ajornet

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','')
    ENDIF

With the follwing output:

 hydrol_soil_tridiag:: (1st) ji,jsl,jst=           2           1           1
 hydrol_soil_tridiag:: (1st) bet(ji)=   17.7660142491691
 hydrol_soil_tridiag:: (1st) rhs(ji,1)= -0.149703576324796
 hydrol_soil_tridiag:: (1st) mcl= -8.426401905638281E-003
 hydrol_soil_tridiag:: (1st) ji,jsl,jst=          35           1           1
 hydrol_soil_tridiag:: (1st) bet(ji)=   9.42310026602285
 hydrol_soil_tridiag:: (1st) rhs(ji,1)= -0.193420292415905
 hydrol_soil_tridiag:: (1st) mcl= -2.052618426584357E-002
 hydrol_soil_tridiag:: (1st) ji,jsl,jst=          36           1           1
 hydrol_soil_tridiag:: (1st) bet(ji)=   11.8512022942869
 hydrol_soil_tridiag:: (1st) rhs(ji,1)= -0.177970189346963
 hydrol_soil_tridiag:: (1st) mcl= -1.501705775731777E-002
 hydrol_soil_tridiag:: (1st) ji,jsl,jst=          57           1           1
 hydrol_soil_tridiag:: (1st) bet(ji)=   6.27096361655699
 hydrol_soil_tridiag:: (1st) rhs(ji,1)= -0.133282335692959
 hydrol_soil_tridiag:: (1st) mcl= -2.125388438565630E-002
 hydrol_soil_tridiag:: (1st) ji,jsl,jst=          58           1           1
 hydrol_soil_tridiag:: (1st) bet(ji)=   41.4191914613879
 hydrol_soil_tridiag:: (1st) rhs(ji,1)= -7.470394780738840E-002
 hydrol_soil_tridiag:: (1st) mcl= -1.803607100274508E-003

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

comment:2 Changed 3 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 3 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 2 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

Agnès for further analysis:

  • check if normal that b(1) and b(2) are negative
  • ask help form Pascal Maugis regarding the numerics
  • can it be related to teh pbs of ae_ns (Tickets #438 and #325)?
  • can it be alleviated with rsoil?

comment:5 Changed 14 months ago by luyssaert

  • Milestone changed from ORCHIDEE 4.0 to Not scheduled yet
Note: See TracTickets for help on using tickets.