cvsp_integrate_volume.f

Go to the documentation of this file.

!           ***********************************************
            DOUBLE PRECISION FUNCTION cvsp_integrate_volume
!           ***********************************************
!
     &(j,i,z_high,z_low,a)
!
!***********************************************************************
! SISYPHE   V7P2                                   16/05/2017
!***********************************************************************
!
!brief   INTEGRATES THE VOLUME OF A FRACTION WITHIN THE
!+       VERTICAL SORTING PROFIL BETWEEN 2 DEPTH Z-COORDINATES Z_HIGH & Z_LOW
!+
!
!history UWE MERKEL
!+        2011
!+        V6P2
!+
!
!history  P. A. TASSI (EDF R&D, LNHE)
!+        12/03/2013
!+        V6P3
!+   Cleaning, cosmetic
!
!history UWE MERKEL, R. KOPMANN
!+        2016, 2017
!+        V6P3, V7P2
!+        improved stability
!
!history  R. KOPMANN
!+        2018
!+        V7P2
!+        max. number of fractions
!
!history  R. KOPMANN (BAW)
!+        25/02/2019
!+        V7P2
!+   Removing 1/NSICLA
!+   no correction of truncation error moved to new subroutine cvsp_check_mass_bilan
!+   no check_f
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| J              |<--| INDEX OF A POINT IN MESH
!| I              |<--| INDEX OF A FRACTION IN VERTICAL SORTING PROFILE
!| Z_HIGH         |<--| HIGHER DEPTH COORDINATE
!| Z_LOW          |<--| LOWER  DEPTH COORDINATE
!| A1...A_NSICLA  |<--| INTEGRATED VOLUMES PER FRACTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe,
     & ex_cvsp_integrate_volume => cvsp_integrate_volume
      USE declarations_sisyphe
      USE cvsp_outputfiles, ONLY: cp
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,          INTENT(IN)  :: J
      INTEGER,          INTENT(IN)  :: I
      DOUBLE PRECISION, INTENT(IN)  :: Z_HIGH
      DOUBLE PRECISION, INTENT(IN)  :: Z_LOW
      DOUBLE PRECISION, INTENT(OUT) :: A(nsicla)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION TEMP,TEMPOLD, DHIG, DLOW, AT, SUMUP, FLOW, FHIG
      DOUBLE PRECISION TEMP2, TEMP2MAX, SUMUP2,TEMP3, TEMP3MAX, SUMUP3
      DOUBLE PRECISION CHSUM
      INTEGER L_CNT, F_CNT, REVCNT, HELPER, LASTCASE, JG
      INTEGER(KIND=K8) :: MYCASE
!
!-----------------------------------------------------------------------
!
      jg = j
      IF (ncsize.GT.1) jg = mesh%KNOLG%I(j)
!
      at = dt*lt/percou
      sumup = 0.d0
      sumup2  = 0.d0
      sumup3  = 0.d0
      mycase  = 0
!
      DO f_cnt = 1, nsicla
        a(f_cnt) = 0.d0
      END DO
      cvsp_integrate_volume = a(i)
      IF ((z_high-z_low).LT.zero) THEN
        RETURN
      ENDIF
!-----------------------------------------------------------------------
! DOING ALL FRACTIONS
!-----------------------------------------------------------------------
      DO f_cnt = 1, nsicla
        temp = 0.d0
        temp2 = 0.d0
        temp2max = 0.d0
        temp3 = 0.d0
        temp3max = 0.d0
        chsum = 0.d0
!-----------------------------------------------------------------------
! GOING THROUGH ALL SECTIONS
!-----------------------------------------------------------------------
        DO l_cnt = 0,(pro_max(j)-2)
          tempold = temp
          revcnt = pro_max(j)-l_cnt
!-----------------------------------------------------------------------
! DEPTH COORDINATES OF THE SECTION TO CHECK
!-----------------------------------------------------------------------
          dhig = pro_d(j,revcnt,f_cnt)
          dlow = pro_d(j,revcnt-1,f_cnt)
!-----------------------------------------------------------------------
! SORTING PROFILE SECTION TOTALLY INSIDE (CASE ZDDZ)
!-----------------------------------------------------------------------
          IF ( (dhig <= z_high ) .AND.
     &         (dlow >= z_low  ) ) THEN
            flow = pro_f(j,revcnt-1,f_cnt)
            fhig = pro_f(j,revcnt,f_cnt)
            mycase  = mycase  + 1
            lastcase = 1
            temp2 = 0.5d0*(fhig+flow)*(dhig-dlow)  + temp2
            temp2max = z_high - dlow
            DO helper = 1, nsicla
              chsum = pro_f(j,revcnt, helper) + chsum
            ENDDO
            chsum = 1.d0 - chsum
!-----------------------------------------------------------------------
! SORTING PROFILE SECTION PARTIALLY LOWER (CASE ZDZD)
!-----------------------------------------------------------------------
          ELSEIF ((dhig <= z_high) .AND.
     &            (dhig >  z_low ) .AND.
     &            (dlow <  z_low ) )  THEN

            fhig = pro_f(j,revcnt,f_cnt)
            flow = pro_f(j,revcnt-1,f_cnt)
            flow = - ((fhig-flow)/(dhig-dlow))*(dhig-z_low) + fhig
!-----------------------------------------------------------------------
! CUT THE SECTION
!-----------------------------------------------------------------------
            dlow = z_low
            mycase  = mycase  + 1000
            lastcase = 2
            temp3 = 0.5d0*(fhig+flow)*(dhig-dlow)  + temp3
            temp3max = dhig - dlow
!-----------------------------------------------------------------------
! SORTING PROFILE SECTION PARTIALLY HIGHER (CASE DZDZ)
!-----------------------------------------------------------------------
!
          ELSEIF ((dhig > z_high) .AND.
     &            (dlow >= z_low) .AND.
     &            (dlow <  z_high) ) THEN
            flow = pro_f(j,revcnt-1,f_cnt)
            fhig = pro_f(j,revcnt,f_cnt)
            fhig = ((fhig-flow)/(dhig-dlow))*(z_high-dlow) + flow
!-----------------------------------------------------------------------
!INSERT
!-----------------------------------------------------------------------
!CUT THE SECTION
            dhig = z_high
            mycase  = mycase + 1000000
            lastcase = 3
!-----------------------------------------------------------------------
! LAYER TOTALLY INSIDE ONE SECTION (CASE DZZD)
!-----------------------------------------------------------------------
          ELSEIF ((dhig >= z_high) .AND.
     &            (dlow <= z_low) ) THEN
            fhig =
     &         - (pro_f(j,revcnt,f_cnt)-pro_f(j,revcnt-1,f_cnt)) /
     &         (dhig-dlow) *
     &         (dhig-z_high)
     &         + pro_f(j,revcnt,f_cnt)
            flow =
     &         - (pro_f(j,revcnt,f_cnt)-pro_f(j,revcnt-1,f_cnt)) /
     &         (dhig-dlow) *
     &         (dhig-z_low)
     &         + pro_f(j,revcnt,f_cnt)
!-----------------------------------------------------------------------
! CUT THE SECTION
!-----------------------------------------------------------------------
            dhig = z_high
            dlow = z_low
            lastcase = 8
            mycase  = mycase + 100000000
!-----------------------------------------------------------------------
! SECTION WITH 0 STRENGTH
!-----------------------------------------------------------------------
          ELSEIF (dhig == dlow) THEN
            flow = 0.d0
            fhig = 0.d0
            mycase  = mycase + 1000000000
            lastcase = 4
!A TRUE BUG!
          ELSEIF (z_low.GE.z_high) THEN
            flow = 0.d0
            fhig = 0.d0
            IF(cp)WRITE(lu,*)'Z_LOW>=Z_HIGH',dhig,dlow,z_high,z_low
            CALL cvsp_p('./','ZLOHI',jg)
            mycase  = mycase + 1000000000
            lastcase = 5
!A TRUE BUG!
          ELSEIF (dhig < dlow) THEN
            flow = 0.d0
            fhig = 0.d0

            IF(cp)WRITE(lu,*)'DHIG<=DLOW',j,dhig,dlow,z_high,z_low
            CALL cvsp_p('./','DLOHI',jg)
            mycase  = mycase + 1000000000
            lastcase = 6
!A SECTION THAT IS NOT INVOLVED / NOT A BUG!!
          ELSE
            flow = 0.d0
            fhig = 0.d0
            mycase  = mycase + 1000000000
            lastcase = 7
          ENDIF
!-----------------------------------------------------------------------
! TRAPEZOID FORMULA
!-----------------------------------------------------------------------
          temp = 0.5d0*(fhig+flow)*(dhig-dlow) + temp
!DEBUG
          IF (0.5d0*(fhig+flow)*(dhig-dlow) < 0.d0) THEN
            WRITE(lu,fmt='(A,1X,2(I11,1X),11(G20.10,1X),1X,I11)')
     &         'INTEGRATE_VOL_ER_TMP:<0:'
     &         ,jg, i, at, fhig,flow,dhig,dlow, dhig-dlow, revcnt,
     &         pro_f(j,revcnt-1,f_cnt),pro_f(j,revcnt,f_cnt),
     &         pro_d(j,revcnt-1,f_cnt),pro_d(j,revcnt,f_cnt),
     &         lastcase
            CALL cvsp_p('./','IVKT',jg)
            CALL plante(1)
          ENDIF
!
        ENDDO                 !SECTION
!-----------------------------------------------------------------------
! ADDING UP FRACTIONS FOR DEBUGGING PURPOSES
!-----------------------------------------------------------------------
        sumup = temp + sumup
        sumup2 = temp2 + sumup2
        sumup3 = temp3 + sumup3
        a(f_cnt) = temp
      ENDDO                    !FRACTION
      cvsp_integrate_volume = a(i)
!
!     DEBUG
      IF (cvsp_integrate_volume < 0.d0) THEN
        CALL cvsp_p('./','IVK0',jg)
        WRITE(*,fmt='(A,2(I11),14(G20.10))')'INTEGRATE_VOL_ER:<0:'
     &       ,jg, i, at, cvsp_integrate_volume, mycase ,z_high,z_low,
     &       sumup,sumup2,sumup3,chsum,a(1),a(2),a(3),a(4),a(5)
        CALL plante(1)
      ENDIF
!-----------------------------------------------------------------------
! CHECKSUM OVER ALL FRACTIONS AND LAYERS
!-----------------------------------------------------------------------
      sumup  = sumup  - abs(z_high-z_low)
      sumup2 = sumup2 - abs(temp2max)
      sumup3 = sumup3 - abs(temp3max)
      chsum =  chsum / 5.d0
!-----------------------------------------------------------------------
! BEWARE: 1.D-5 IS UP TO 10G OF A 1000KG BUT HIGHER ACCURACY
! LEADS AGAIN TO FLOATING POINT TRUNCATION ERRORS ...
!-----------------------------------------------------------------------
      IF (abs(sumup).GT.1.d-5) THEN
        IF(cp) CALL cvsp_p('./','IV_E',jg)
        IF(cp)WRITE(lu,*) 'INTEGRATE VOLUME ACCURRACY!!!',sumup,jg
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END FUNCTION cvsp_integrate_volume