! Replica as equacoes de momentum portadas (KFPARA.f90 KFMX) com perfis
! sinteticos consistentes com a continuidade do KF, e testa a conservacao de
! momentum da coluna: SUM[(UPA-U0)*EMS] deve ser ~0 (conveccao so redistribui).
      PROGRAM CONS
      IMPLICIT NONE
      INTEGER, PARAMETER :: KX=50
      INTEGER :: NK,NK1,NU,NU1,ND,ND1,NTC,NSTEP,LTOP,LET,KSTART,LFS,LDT,LDB
      REAL :: UMF(KX),UER(KX),UDR(KX),DMF(KX),DER(KX),DDR(KX),U0(KX),V0(KX)
      REAL :: UMFV(KX),UERV(KX),UDRV(KX),DMFV(KX),DERV(KX),DDRV(KX)
      REAL :: UUP(KX),VUP(KX),UDW(KX),VDW(KX),UPA(KX),VPA(KX),UG(KX),VG(KX)
      REAL :: UFXIN(KX),UFXOUT(KX),VFXIN(KX),VFXOUT(KX),DOMGVDP(KX),FXMV(KX)
      REAL :: OMGV(KX+1),DP(KX),EMS(KX),EMSD(KX),P0(KX)
      REAL :: G,DXSQ,DTIME,TIMEC,SM0,SMF,SCALE,DUMF,DPDD,UCONST,DPTHMX,BAL
      INTEGER :: MODE
      CALL GET_MODE(MODE)
      UCONST=10.
      G=9.81; DXSQ=20000.**2; TIMEC=1800.; NSTEP=20; DTIME=TIMEC/REAL(NSTEP)
! --- grade e ambiente (NK cresce p/ CIMA, como no KFPARA) ---
      DO NK=1,KX
        P0(NK)=100000.-(NK-1)*1800.          ! 1000 -> 118 hPa
        DP(NK)=1800.
        EMS(NK)=DP(NK)*DXSQ/G
        EMSD(NK)=1./EMS(NK)
        IF(MODE == 1)THEN
          U0(NK)=UCONST; V0(NK)=UCONST      ! campo CONSTANTE: operador conservativo
        ELSE                                 ! deve devolve-lo intacto
          U0(NK)=5.+0.5*(NK-1); V0(NK)=-2.+0.1*(NK-1)
        ENDIF
        UMF(NK)=0.; UER(NK)=0.; UDR(NK)=0.
        DMF(NK)=0.; DER(NK)=0.; DDR(NK)=0.
      ENDDO
      KSTART=5; LET=30; LTOP=40; LFS=32; LDT=KSTART-1; LDB=1
! --- updraft: UMF cresce por entranhamento ate LET, decai por detrainment ---
! camada-fonte: a massa da base do updraft e' ENTRANHADA (UER) em 1..KSTART,
! como no KFPARA:1230-1233 (UER=VMFLCL*DP/DPTHMX; UMF(NK)=UMF(NK-1)+UER(NK)).
! Sem isto a coluna nao fecha o balanco de massa.
      DPTHMX=0.
      DO NK=1,KSTART
        DPTHMX=DPTHMX+DP(NK)
      ENDDO
      DO NK=1,KSTART
        UER(NK)=1.0E8*DP(NK)/DPTHMX
        IF(NK == 1)THEN
          UMF(NK)=UER(NK)
        ELSE
          UMF(NK)=UMF(NK-1)+UER(NK)
        ENDIF
      ENDDO
      DO NK=KSTART+1,LTOP
        IF(NK <= LET)THEN
          UER(NK)=0.06*UMF(NK-1); UDR(NK)=0.01*UMF(NK-1)
        ELSE
          UER(NK)=0.0;            UDR(NK)=0.22*UMF(NK-1)
        ENDIF
        UMF(NK)=UMF(NK-1)+UER(NK)-UDR(NK)
        IF(NK == LTOP)THEN
          UDR(NK)=UMF(NK-1); UER(NK)=0.; UMF(NK)=0.
        ENDIF
      ENDDO
! --- downdraft: DMF<0, comeca no LFS, some no LDB ---
      DMF(LFS)=-0.4*UMF(KSTART)
      DER(LFS)=DMF(LFS)   ! KFPARA:1486 - massa inicial do downdraft e' entranhada no LFS
! in-cloud: LFS-1 -> KSTART, so entranhamento (DER<0), DMF(ND)=DMF(ND1)+DER(ND)
      DO ND=LFS-1,KSTART,-1
        DER(ND)=0.030*DMF(LFS)*EMS(ND)/EMS(LFS); DDR(ND)=0.
        DMF(ND)=DMF(ND+1)+DER(ND)
      ENDDO
! abaixo da base: LDT -> LDB, so detrainment (DDR>0), DMF(ND)=DMF(ND1)+DDR(ND)
      DPDD=0.
      DO ND=LDB,LDT
        DPDD=DPDD+DP(ND)
      ENDDO
      DO ND=LDT,LDB,-1
        DDR(ND)=-DMF(KSTART)*DP(ND)/DPDD; DER(ND)=0.
        DMF(ND)=DMF(ND+1)+DDR(ND)
      ENDDO
!======================= BLOCO PORTADO (KFPARA.f90) =======================
      DO NK=1,KX
        UMFV(NK)=0.; UERV(NK)=0.; UDRV(NK)=0.
        DMFV(NK)=0.; DERV(NK)=0.; DDRV(NK)=0.
        UDW(NK)=0.;  VDW(NK)=0.;  UUP(NK)=0.; VUP(NK)=0.
        OMGV(NK)=0.; FXMV(NK)=0.
      ENDDO
      OMGV(KX+1)=0.
! downdraft UDW/VDW
      IF(ABS(DMF(KSTART)) > 1.E-3)THEN
        UDW(LFS)=U0(LFS); VDW(LFS)=V0(LFS)
        DO ND=LFS-1,KSTART,-1
          ND1=ND+1
          IF(ND == (LFS-1))THEN
            DMFV(ND1)=DMF(LFS); DERV(ND1)=DER(LFS); DDRV(ND1)=DDR(LFS)
          ENDIF
          DMFV(ND)=DMF(ND); DERV(ND)=DER(ND); DDRV(ND)=DDR(ND)
          IF(ABS(DMFV(ND)) > 1.E-3)THEN
            UDW(ND)=(UDW(ND1)*DMFV(ND1)+0.5*(U0(ND)+U0(ND1))*DERV(ND))/DMFV(ND)
            VDW(ND)=(VDW(ND1)*DMFV(ND1)+0.5*(V0(ND)+V0(ND1))*DERV(ND))/DMFV(ND)
          ELSE
            UDW(ND)=UDW(ND1); VDW(ND)=VDW(ND1)
          ENDIF
        ENDDO
        DO ND=LDT,LDB,-1
          ND1=ND+1
          IF(ND == LDT)THEN
            DMFV(ND1)=DMF(KSTART); DDRV(ND1)=DDR(KSTART); DERV(ND1)=DER(KSTART)
          ENDIF
          DMFV(ND)=DMF(ND); DDRV(ND)=DDR(ND); DERV(ND)=DER(ND)
          IF(ND == LDB)THEN
            DDRV(ND)=DMFV(ND)-DMFV(ND1); DMFV(ND)=0.
          ENDIF
        ENDDO
        DO ND=LDT,LDB,-1
          ND1=ND+1
          IF(ABS(2.*DMFV(ND)-DDRV(ND)) >= 1.E-3)THEN
            UDW(ND)=UDW(ND1)*((2.*DMFV(ND1)+DDRV(ND))/(2.*DMFV(ND)-DDRV(ND)))
            VDW(ND)=VDW(ND1)*((2.*DMFV(ND1)+DDRV(ND))/(2.*DMFV(ND)-DDRV(ND)))
          ELSE
            UDW(ND)=UDW(ND1); VDW(ND)=VDW(ND1)
          ENDIF
        ENDDO
      ENDIF
! updraft UUP/VUP + init UPA/VPA
      DO NK=1,LTOP
        UPA(NK)=U0(NK); VPA(NK)=V0(NK)
        IF(NK < KSTART)THEN
          UUP(NK)=U0(NK); VUP(NK)=V0(NK)
          UMFV(NK)=UMF(NK); UERV(NK)=UER(NK); UDRV(NK)=UDR(NK)
        ELSE
          UUP(NK)=U0(KSTART); VUP(NK)=V0(KSTART)
        ENDIF
      ENDDO
      DO NK=LTOP+1,KX
        UFXIN(NK)=0.; UFXOUT(NK)=0.; VFXIN(NK)=0.; VFXOUT(NK)=0.
        UPA(NK)=U0(NK); VPA(NK)=V0(NK); UG(NK)=U0(NK); VG(NK)=V0(NK)
      ENDDO
      IF(ABS(UMF(KSTART)) > 1.E-3)THEN
        DO NK=KSTART,LET-1
          NK1=NK+1
          UMFV(NK)=UMF(NK); UERV(NK)=UER(NK); UDRV(NK)=UDR(NK)
          IF(NK == (LET-1))THEN
            UMFV(NK1)=UMF(LET); UERV(NK1)=UER(LET); UDRV(NK1)=UDR(LET)
          ENDIF
        ENDDO
        DO NK=KSTART,LET-1
          NK1=NK+1
          IF(ABS(2.*UMFV(NK1)+UDRV(NK1)) < 1.E-3)THEN
            UUP(NK1)=UUP(NK); VUP(NK1)=VUP(NK)
          ELSE
            UUP(NK1)=(UUP(NK)*(2.*UMFV(NK)-UDRV(NK1))+UERV(NK1)*(U0(NK)+U0(NK1)))/(2.*UMFV(NK1)+UDRV(NK1))
            VUP(NK1)=(VUP(NK)*(2.*UMFV(NK)-UDRV(NK1))+UERV(NK1)*(V0(NK)+V0(NK1)))/(2.*UMFV(NK1)+UDRV(NK1))
          ENDIF
        ENDDO
        DO NU=LET,LTOP-1
          NU1=NU+1
          UMFV(NU)=UMF(NU); UDRV(NU)=UDR(NU); UERV(NU)=UER(NU)
          IF(NU == LTOP-1)THEN
            UDRV(NU1)=UMFV(NU)-UMFV(NU1); UMFV(NU1)=0.
          ENDIF
          IF(ABS(2.*UMFV(NU1)+UDRV(NU1)) < 1.E-3)THEN
            UUP(NU1)=UUP(NU); VUP(NU1)=VUP(NU)
          ELSE
            UUP(NU1)=(UUP(NU)*(2.*UMFV(NU)-UDRV(NU1))+UERV(NU1)*(U0(NU)+U0(NU1)))/(2.*UMFV(NU1)+UDRV(NU1))
            VUP(NU1)=(VUP(NU)*(2.*UMFV(NU)-UDRV(NU1))+UERV(NU1)*(V0(NU)+V0(NU1)))/(2.*UMFV(NU1)+UDRV(NU1))
          ENDIF
        ENDDO
      ENDIF
! laco temporal
      DO NTC=1,NSTEP
        DO NK=1,LTOP
          UFXIN(NK)=0.; UFXOUT(NK)=0.; VFXIN(NK)=0.; VFXOUT(NK)=0.
        ENDDO
        DO NK=1,LTOP
          IF((NK >= KSTART) .AND. (NK <= LFS))THEN
            DDRV(NK)=DDR(NK); DERV(NK)=DER(NK)
          ENDIF
          IF(NK <= LET)THEN
            UDRV(NK)=UDR(NK); UERV(NK)=UER(NK)
          ENDIF
          DOMGVDP(NK)=-(UERV(NK)-DERV(NK)-UDRV(NK)-DDRV(NK))*EMSD(NK)
          IF(NK > 1)THEN
            OMGV(NK)=OMGV(NK-1)-DP(NK-1)*DOMGVDP(NK-1)
          ELSE
            OMGV(NK)=0.
          ENDIF
          FXMV(NK)=OMGV(NK)*DXSQ/G
        ENDDO
        DO NK=LTOP+1,KX
          FXMV(NK)=0.
        ENDDO
        DO NK=2,LTOP
          IF(OMGV(NK) <= 0.)THEN
            UFXIN(NK)=-FXMV(NK)*UPA(NK-1); VFXIN(NK)=-FXMV(NK)*VPA(NK-1)
            UFXOUT(NK-1)=UFXOUT(NK-1)+UFXIN(NK); VFXOUT(NK-1)=VFXOUT(NK-1)+VFXIN(NK)
          ELSE
            UFXOUT(NK)=FXMV(NK)*UPA(NK); VFXOUT(NK)=FXMV(NK)*VPA(NK)
            UFXIN(NK-1)=UFXIN(NK-1)+UFXOUT(NK); VFXIN(NK-1)=VFXIN(NK-1)+VFXOUT(NK)
          ENDIF
        ENDDO
        DO NK=1,LTOP
          UPA(NK)=UPA(NK)+(UFXIN(NK)+UDRV(NK)*UUP(NK)+DDRV(NK)*UDW(NK) &
                  -UFXOUT(NK)-(UERV(NK)-DERV(NK))*U0(NK))*DTIME*EMSD(NK)
          VPA(NK)=VPA(NK)+(VFXIN(NK)+UDRV(NK)*VUP(NK)+DDRV(NK)*VDW(NK) &
                  -VFXOUT(NK)-(UERV(NK)-DERV(NK))*V0(NK))*DTIME*EMSD(NK)
        ENDDO
      ENDDO
      DO NK=1,LTOP
        UG(NK)=UPA(NK); VG(NK)=VPA(NK)
      ENDDO
!========================= DIAGNOSTICO ==================================
      BAL=0.
      DO NK=1,LTOP
        BAL=BAL+(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))
      ENDDO
      WRITE(*,'(A,ES12.4,A,ES10.2,A)') ' balanco de massa da coluna SUM(UER-DER-UDR-DDR) = ', &
            BAL,'   (escala UMF base = ',UMF(KSTART),')'
      SM0=0.; SMF=0.; SCALE=0.
      DO NK=1,KX
        SM0=SM0+U0(NK)*EMS(NK)
        SMF=SMF+UG(NK)*EMS(NK)
        SCALE=SCALE+ABS(UG(NK)-U0(NK))*EMS(NK)
      ENDDO
      WRITE(*,'(A,ES14.6)') ' momentum U inicial (kg m/s) = ', SM0
      WRITE(*,'(A,ES14.6)') ' momentum U final   (kg m/s) = ', SMF
      WRITE(*,'(A,ES12.4,A,F9.5,A)') ' residuo = ', SMF-SM0, &
            '   (', 100.*ABS(SMF-SM0)/SCALE, ' % do transporte bruto)'
      WRITE(*,'(A)') ' perfil (NK, p_hPa, U0, UG, dU):'
      DO NK=LTOP+2,1,-4
        IF(NK <= KX) WRITE(*,'(I4,F9.1,3F9.3)') NK,P0(NK)*0.01,U0(NK),UG(NK),UG(NK)-U0(NK)
      ENDDO
      WRITE(*,'(A,2F12.6)') ' dU max/min = ', MAXVAL(UG(1:LTOP)-U0(1:LTOP)), MINVAL(UG(1:LTOP)-U0(1:LTOP))
      IF(MODE == 1)THEN
        WRITE(*,'(A,ES12.4,A)') ' TESTE DE PRESERVACAO DE CONSTANTE: desvio max = ', &
              MAXVAL(ABS(UG(1:LTOP)-UCONST)), ' m/s  (deve ser ~0)'
        DO NK=1,LTOP
          IF(ABS(UG(NK)-UCONST) > 1.E-3) WRITE(*,'(A,I4,A,F12.5)') '   NK=',NK,'  UG=',UG(NK)
        ENDDO
      ENDIF
      END PROGRAM CONS
      SUBROUTINE GET_MODE(M)
      INTEGER M
      CHARACTER(8) A
      M=1
      IF(COMMAND_ARGUMENT_COUNT() >= 1)THEN
        CALL GET_COMMAND_ARGUMENT(1,A); READ(A,*) M
      ENDIF
      END SUBROUTINE GET_MODE