!-----------------------------------------------------------------------------!
!             Modeling and Development Division - DMD/CPTEC/INPE              !
!-----------------------------------------------------------------------------!
!BOI
! !TITLE: humNphy
!
! !AUTHORS: Luiz F. Sapucci based on fct2anl_trans.f90 from João Gerd Zell de Mattos
!
! !AFFILIATION: Modeling and Development Division, CPTEC/INPE
!
! !DATE:  03 de abril de 2017
!
! !DESCRIPTION:
!      Compiler of the evaluation protocol of the humidity not physical values
!
! ! INOUT FILE:
!  input argument list:
!     BAM.fct.hh.bin      - binary file of the BAM forecast of hh with field required  
!     BAM.fct.06.bin.ctl  - ctl file for read the binary file created.
!
!   output files
!
! !REMARKS:
!     (1)
!
! !USES:
!
!
!EOI
!-----------------------------------------------------------------------------!
!


program hunNphy

!
! !USES:
!

   use sigio_BAMMod, only : BAM_Open,            & ! Open a BAM file
                            BAM_GetField,        & ! Read BAM field from a file
                            BAM_GetDims,         & ! Return BAM field dimensions
                            BAM_GetTimeInfo,     & ! Return time info
                            BAM_GetNlevels,      & ! Return the number of leves of a field
                            BAM_WriteAnlHeader,  & ! Write the Anl file header
                            BAM_GetWCoord,       & ! Return the vector of cordenate of a field
                            BAM_WriteField,      & ! Write a BAM field (mpi or serial)
                            BAMFile,             & ! BAM data type structure
                            BAM_GetVerticalCoord,& ! Informations vertical coordinates from BAM model.
                            r8,                  & ! Kind for 64-bits Real Numbers
                            r4                     ! Kind for 64-bits Real Numbers

   USE scam_MetForm                                ! Module to convert meteorological variables

   implicit none

   integer :: istat

   integer :: IMax, i, j,k,d
   integer :: JMax
   integer :: KMax
   integer :: Mend
   integer :: Nvar

   real    :: mwave2
   integer :: mnwv2, mnwv3
   
   integer :: mydate(5) 
   
   real    :: URnegCONTdomi(5),URsatCONTdomi(5),URnegVALdomi(5),URsatVALdomi(5)

   integer :: ivar
   integer :: ilev
   integer :: Nlevs
   integer :: NVars,n

   type(BAMfile) :: BAM

   real(kind=8) :: Rm,Esat,PA,ur,TVirt,Uesp,tempA,TVconf,UespConf,eeconf,esconf
   real :: LnPaMax,LnPaMin
 
   real(r8), allocatable :: field(:),umimax(:),umimin(:),BK(:)
   real(r4), allocatable :: grid(:,:),LnPa(:,:),TempVir(:,:,:),UmiEsp(:,:,:),UmiRela(:,:,:),press(:,:)

   real(r4), allocatable :: URnegCONT(:,:),URsatCONT(:,:),URnegVAL(:,:),URsatVAL(:,:)
   real(r4), allocatable :: URnegV(:),URsatV(:),rlat(:)
  
   integer, allocatable :: ar1(:)
   integer, dimension(6) :: maskDomi

   integer, parameter :: stderr = 0
   integer, parameter :: stdinp = 5
   integer, parameter :: stdout = 6
   parameter(maskDomi=[-90,-60,-20,20,60,90])

   character(len=20) :: myname_='humNphy'
   character(len=20) :: spcFILEname,dirFILEname,ctlTemplate
   character(len=25) :: binName,ctlName,binNameHum,ctlNameHum
   character(len=82)  :: mensa = '#'
   character(len=40), dimension(30) :: VName,Vlable

   ! Open files config (the same used spc2grb.x)
      
   OPEN (UNIT=9,FILE='spc2grd.config',STATUS='old')

   !
   ! lendo o arquivo de configuracao e carregando os parametros
   !

   do while (mensa(1:1) == '#' .OR. mensa(1:1) == '$' .OR. mensa(1:1) == ' ')
     read(9,'(A40)')mensa
!     write(*,*)'dentro loop',mensa
   enddo 
 
   backspace(9)
   
   read(9,'(A13,A15)')mensa,spcFILEname ! Nome do arquivo de entrada espectral previsao ou analise
   read(9,'(A13,A15)')mensa,dirFILEname ! Nome do arquivo de directivas do arquivo espectral de previsao ou analise
   read(9,'(A13,A31)')mensa,ctlTemplate ! Nome do arquivo ctl template nao usado aqui
   
   read(9,'(A6,I4)')mensa,NVars

   do Nvar = 1, NVars
     read(9,'(A4,2x,A26)')Vlable(Nvar),VName(Nvar)
   enddo

   ! Make input and output file name and open them  

   N = 1
   do while (spcFILEname(N:N) /= ' ')
      N=N+1
   enddo 

   binName=spcFILEname(1:N-1)//'.bin'
   ctlName=spcFILEname(1:N-1)//'.bin.ctl' 

   OPEN (UNIT=13,FILE=binName,status='old',form='unformatted')
   OPEN (UNIT=15,FILE=ctlName,status='old')
  
   binNameHum=spcFILEname(1:N-1)//'.RH.bin'
   ctlNameHum=spcFILEname(1:N-1)//'.RH.bin.ctl'

   OPEN (UNIT=16,FILE=binNameHum,status='unknown',form='unformatted')
   OPEN (UNIT=17,FILE=ctlNameHum,status='unknown')

   !
   ! Files to read and get first-guess or analysis spectral format and read header
   !

   BAM%ffct = spcFILEname ! Arquivo de spectral a ser lido anl ou fct
   BAM%fdir = dirFILEname ! Arquivo descritor do fct

   write (*,*)
   write (*,*) "-------------------------------------"
   write (*,*) " Evaluation GSI: rotine humHphy:  "
   write (*,*) "-------------------------------------"
   write (*,*)


   call BAM_Open(BAM, 'r', istat=istat)
   if (istat.ne.0)then
      write(*,*)'Problem to read BAM files!'
   endif

   ! Get BAM dimensions
   call BAM_GetDims(BAM, imax, jmax, kmax, Mend, istat)
   
   mnwv2 = (Mend+1)*(Mend+2)
   mnwv3 = (mnWV2+2)*(Mend+1)
   write(stdout,*)'IMAx',IMax
   write(stdout,*)'JMax',JMax
   write(stdout,*)'KMax',KMax
   write(stdout,*)'Mend',Mend
   write(stdout,*)'mnWV2',mnwv2
   write(stdout,*)

   allocate(BK(Kmax))
   call BAM_GetVerticalCoord(BAM, 'BK', BK, istat)

   !
   ! Set time from fct file to put in IC header
   !
   mydate(1:8) = 0.0
   ! Analysis date
   mydate(1) = BAM_GetTimeInfo(BAM,'iyr') ! Year of Century
   mydate(2) = BAM_GetTimeInfo(BAM,'imo') ! Month of year
   mydate(3) = BAM_GetTimeInfo(BAM,'idy') ! Day of month
   mydate(4) = BAM_GetTimeInfo(BAM,'ihr') ! Hour of day
   ! Forecast date
   mydate(5) = BAM_GetTimeInfo(BAM,'fyr') ! Year of Century
   mydate(6) = BAM_GetTimeInfo(BAM,'fmo') ! Month of year
   mydate(7) = BAM_GetTimeInfo(BAM,'fdy') ! Day of month
   mydate(8) = BAM_GetTimeInfo(BAM,'fhr') ! Hour of day

   write (*,*) " Input files:"
   write (*,*) "CONFIG file   : ",'spc2grd.config'
   write (*,*) "BIN File      : ",binName
   write (*,*) "CTL File      : ",ctlName
   write (*,*) "SPC File info : ",spcFILEname
   write (*,*)
   write (*,*) " Output files:"
   write (*,*) "BIN File output: ",binNameHum
   write (*,*) "CTL File output: ",ctlNameHum
   write (*,*)
   write (*,'(A30,I4,3(I2))') " Analysis date(YYYYMMDDHH): ",mydate(1:4)
   write (*,'(A30,I4,3(I2))') " Forecast date(YYYYMMDDHH): ",mydate(5:8)
   write (*,*)

! Generating the new ctl file for Relative humidity bin file

   write(17,'(A6,A20)')'DSET ^',binNameHum
   read(15,'(A82)')mensa   
   do while (mensa(1:4) /= 'TDEF')
     read(15,'(A82)')mensa
     write(17,'(A82)')mensa
   ENDDO
 
   write(17,'(A4,1X,I2)')'VARS',8

   !
   ! Lendo dos dados LnPres, TV, Q no arqvivo grib
   !

   ! Pulando os dados da topografia no arquivo grib
   read(13) 

   ! Lendo os dados da pressão a superficie 
   ivar   = 2
   allocate(LnPa(imax,jmax))
   allocate(Press(imax,jmax))
   read(13) LnPa
   LnPaMax=MAXVAL(LnPa)
   LnPaMin=MINVAL(LnPa)  
   write(*,*)
   write(*,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,'(A20,A20,A25)')'%%%%%%%%%%%%%%     ',trim(VName(ivar)),'          %%%%%%%%%%%%%%%'
   write(stdout,*)"         level       Maximum values             Minimum  values "
   write(stdout,*) ivar,LnPaMax,LnPaMin
   LnPaMax=EXP(MAXVAL(LnPa))*10
   LnPaMin=EXP(MINVAL(LnPa))*10
   write(stdout,'(A62)')'%%%%%%%%%%%%%%    Pressao Atmosférica          %%%%%%%%%%%%%%%'
   write(stdout,*) 2,LnPaMax,LnPaMin

   PA=EXP(LnPa(100,50))*10

   write(16)EXP(LnPa)*10
   write(17,'(A4,I3,I3,2X,A26)')Vlable(ivar),1,99,'Presure at surface [hPa]' ! Variable information in the ctl file

   ! Lendo os dados da temperatura virtual 
   ivar   = 3

   allocate(umimax(KMax))
   allocate(umimin(KMax))

   allocate(TempVir(imax,jmax,kmax))
   allocate(grid(imax,jmax))

   do ilev = 1, KMax
  
       read(13)grid
       write(16) grid
   
       TempVir(:,:,ilev)=grid

       umimax(ilev)=MAXVAL(TempVir(:,:,ilev))
       umimin(ilev)=MINVAL(TempVir(:,:,ilev))      

       if (ivar .eq. 3 .AND. ilev .eq. 1)then; TVirt=grid(100,50); endif

   enddo

   write(*,*)
   write(*,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,'(A20,A20,A25)')'%%%%%%%%%%%%%%     ',trim(VName(ivar)),'          %%%%%%%%%%%%%%%'
   write(stdout,*)"         level       Maximum values             Minimum  values "

   do ilev = 1, KMax 
       write(stdout,*) ilev,umimax(ilev),umimin(ilev)      
   enddo

   write(stdout,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,*)

   deallocate(umimax)
   deallocate(umimin)

   deallocate(grid)
   deallocate(press)

   write(17,'(A4,I3,I3,2X,A26)')Vlable(ivar),KMax,99,'Virtual temperature [K]' ! Variable information in the ctl file

   ! Pulando os kMax niveis da divergencia e vorticidade
   do ilev = 1, KMax*2
       read(13)
   enddo

   ! Lendo os dados da Umidade Especifica 
   ivar   = 6

   allocate(umimax(KMax))
   allocate(umimin(KMax))

   allocate(UmiEsp(imax,jmax,kmax))
   allocate(grid(imax,jmax))

   do ilev = 1, KMax
  
       read(13)grid
       write(16)grid
   
       UmiEsp(:,:,ilev)=grid

       umimax(ilev)=MAXVAL(UmiEsp(:,:,ilev))
       umimin(ilev)=MINVAL(UmiEsp(:,:,ilev))      

       if (ivar .eq. 6 .AND. ilev .eq. 1)then; Uesp=grid(100,50); endif

   enddo

   write(*,*)
   write(*,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,'(A20,A20,A25)')'%%%%%%%%%%%%%%     ',trim(VName(ivar)),'          %%%%%%%%%%%%%%%'
   write(stdout,*)"         level       Maximum values             Minimum  values "

   do ilev = 1, KMax 
       write(stdout,*) ilev,umimax(ilev),umimin(ilev)      
   enddo

   write(stdout,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,*)

   deallocate(umimax)
   deallocate(umimin)

   deallocate(grid)

   write(17,'(A4,I3,I3,2X,A26)')Vlable(ivar),KMax,99,'Specific Humidity [kg/kg]' ! Variable information in the ctl file

   !
   ! Calculando os valores da umidade relativa usando TV, Q e LN(PRES) lidos do grid
   !

   ! Calcula a umidade relativa usando a funçao de conversão do SCANTEC
   ! rh    Umidade relativa [%]
   ! w     Razão de Mistura [-]
   ! es    Presão de Vapor Saturada [Pa]
   ! pres  Presão Atmosférica [Pa]
   ! t     Temperatura [Graus Celsius][C]
   ! q     Umidade especifica [kg/kg]
   ! tv    Temperatura Virtual [C]
   ! Funções usadas do modulo scam_MetForm
   ! rh(w,es,pres)
   ! w(q)
   ! ta(tv,q)
   ! es(t)
   !tv2(t,rh,pres)
   !q3(pres,t,rh)

   Rm=w(Uesp)
   tempA=ta(TVirt,Uesp)
   Esat=es(tempA-273.16)
   
   write(stdout,*)"lnPA TV Uesp Rm tempA Esat (100,50,super): ",PA, TVirt, Uesp,Rm,tempA-273.16,Esat/100

   ur=rh(Rm,Esat,PA*100)

   write(stdout,*)"Rm,Esat,PA,ur (100,50,super): ",Rm,Esat,PA*100,ur

   TVconf=tv(tempA-273.16,ur,PA*100,UespConf,esconf,eeconf)

   write(stdout,*)"tempA,ur,PA/100,TVconf TVlido(100,50,super) TVconf,UespConf,esconf ,eeconf: ",tempA-273.16,ur,PA*100,TVconf+273.16,TVirt,UespConf,esconf,eeconf

   UespConf=q(PA*100,tempA-273.16,ur)

   write(stdout,*)"PA/100,tempA,ur,UespConf UespLIda(100,50,super): ",PA*100,tempA-273.16,ur,UespConf,Uesp

   allocate(umimax(KMax))
   allocate(umimin(KMax))
   allocate(UmiRela(imax,jmax,kmax))

   do k = 1, KMax
     do j = 1, JMax
       do i = 1, IMax

         Rm=w(real(UmiEsp(i,j,k),r8))
         tempA=ta(real(TempVir(i,j,k),r8),real(UmiEsp(i,j,k),r8))
         Esat=es(tempA-273.16)
         PA=EXP(real(LnPa(i,j),r8))*10*BK(k)
         UmiRela(i,j,k)=rh(Rm,Esat,PA*100)

       enddo
     enddo
     write(16) UmiRela(:,:,k)
     write(*,*)k,BK(k)
 
     umimax(k)=MAXVAL(UmiRela(:,:,k))
     umimin(k)=MINVAL(UmiRela(:,:,k)) 

   enddo

   write(*,*)
   write(*,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,'(A60)')'%%%%%%%%%%%%%%     Umidade Relativa          %%%%%%%%%%%%%%%'
   write(stdout,*)"         level       Maximum values             Minimum  values ",UmiRela(100,50,1)

   do ilev = 1, KMax 
       write(stdout,*) ilev,umimax(ilev),umimin(ilev)      
   enddo

   write(stdout,*)"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   write(stdout,*)

   deallocate(umimax)
   deallocate(umimin)

  !
  ! Calculo das estatisticas da ocorrencia da umidade nao fisica
  !

  ! campos espacial de somatórios de pontos com UR não fisica no perfil 
   allocate(URnegCONT(imax,jmax))
   allocate(URsatCONT(imax,jmax))

   ! Campos de somatórios dos valores de UR não fisica  no perfil
   allocate(URnegVAL(imax,jmax))
   allocate(URsatVAL(imax,jmax))

  ! Perfis verticias do somatórios totais de pontos com UR não ficica na horizontal
   allocate(URnegV(kmax))
   allocate(URsatV(kmax))

   do j = 1, JMax
     do i = 1, IMax 
       URnegCONT(i,j)=count(UmiRela(i,j,:).lt.0)
       URsatCONT(i,j)=count(UmiRela(i,j,:).gt.100)
       URnegVAL(i,j)=sum(UmiRela(i,j,:),mask=UmiRela(i,j,:).lt.0)
       URsatVAL(i,j)=sum(UmiRela(i,j,:)-100,mask=UmiRela(i,j,:).gt.100)
     enddo
   enddo

   do k = 1, Kmax
       URnegV(k)=count(UmiRela(:,:,k).lt.0)
       URsatV(k)=count(UmiRela(:,:,k).gt.100)
   enddo

   write (stdout,*)sum(URnegCONT),sum(URsatCONT),sum(URnegVAL),sum(URsatVAL)

   !
   ! Calculando os valores totais por regiões do globo em latitude
   !

   ! Chamando o módulo do BAM_io para tomar o vertor de latitudes do modelo
   allocate(rlat(jmax))
   allocate(ar1(6))
   call BAM_GetWCoord(BAM,'rlat',rlat,istat)
!DEBUG   do d=1,jmax
!DEBUG     write(stdout,*)d,rlat(d),rlat(d)+90
!DEBUG   enddo
   ! Usando minloc para tomar os indices associados com as latitudes MaskDomi=[-90,-60,-20,20,60,90]
   ar1(1)=0   
   do d=2,5
     ar1(d)=minloc(rlat, MASK = rlat .GT. MaskDomi(d),dim=1)
     write(stdout,*)d,ar1(d),rlat(ar1(d)),MaskDomi(d)
   enddo
   ar1(6)=jmax
   write(stdout,*)d,ar1(6),rlat(ar1(6)),MaskDomi(6)

   ! Calculando os valores para os recortes em latitude usando ar1 associado ao MaskDomi=[-90,-60,-20,20,60,90]
   do d=2,6
      ! Contador de casos no recorte
      URnegCONTdomi(d-1)=sum(URnegCONT(:,ar1(d-1)+1:ar1(d)))
   write(stdout,*)ar1(d-1),ar1(d)
      URsatCONTdomi(d-1)=sum(URsatCONT(:,ar1(d-1)+1:ar1(d)))
      ! Valor somatorio dos valores não físicos no recorte
      URnegVALdomi(d-1)=sum(URnegVAL(:,ar1(d-1)+1:ar1(d)))
      URsatVALdomi(d-1)=sum(URsatVAL(:,ar1(d-1)+1:ar1(d)))
   enddo
   
   
   
   ! Escrevendo os valores do perfil e total em recortes no arquivo de log para edicão via script
   ! Valores por niveis rever para outras resolucoes verticais a serem usadas
   if (kmax.eq.28)then
     write(stdout,'(A8,A7,I5,3(I3),28(F10.2))')"#URnegH_",binName,mydate(1:4),URnegV(:)
     write(stdout,'(A8,A7,I5,3(I3),28(F10.2))')"#URsatH_",binName,mydate(1:4),URsatV(:)
     write(*,*)'Plotando 28 niveis do modelo. Rever codigo se resolucao vertical errada!'
   else
     write(stdout,'(A8,A7,I5,3(I3),64(F10.2))')"#URnegH_",binName,mydate(1:4),URnegV(:)
     write(stdout,'(A8,A7,I5,3(I3),64(F10.2))')"#URsatH_",binName,mydate(1:4),URsatV(:)
     write(*,*)'Plotando 64 niveis do modelo. Rever codigo se resolucao vertical errada!'
   endif

   ! Numeros de casos por dominio
      write(stdout,'( A)')"Variable  File_Name   YYYY MM DD HH  Antartico  HemiSul  Tropical  HemiNorte   Artico"
   write(stdout,'(A14,A7,I5,3(I3),5(F10.2))')"#URnegCONdomi_",binName,mydate(1:4),URnegCONTdomi(:)
   write(stdout,'(A14,A7,I5,3(I3),5(F10.2))')"#URsatCONdomi_",binName,mydate(1:4),URsatCONTdomi(:)
   ! Valor total por dominio
   write(stdout,'(A14,A7,I5,3(I3),5(F10.2))')"#URnegVALdomi_",binName,mydate(1:4),URnegVALdomi(:)
   write(stdout,'(A14,A7,I5,3(I3),5(F10.2))')"#URsatVALdomi_",binName,mydate(1:4),URsatVALdomi(:)

   ! Escrevendo resultados no arquivo binario
   write(16)URnegCONT
   write(16)URsatCONT
   write(16)URnegVAL
   write(16)URsatVAL

   ! UR Variables information in the ctl file
   write(17,'(A4,I3,I3,2X,A26)')'Rehu',KMax,99,'Relative humidity [%]' 
   write(17,'(A4,I3,I3,2X,A46)')'Urnc',1,99,'Total Count of Relative humidity negative [n]' 
   write(17,'(A4,I3,I3,2X,A46)')'Ursc',1,99,'Total Count of Relative humidity saturate [n]' 
   write(17,'(A4,I3,I3,2X,A37)')'Urnv',1,99,'Sum of Relative humidity negative [%]' 
   write(17,'(A4,I3,I3,2X,A36)')'Ursv',1,99,'Sum of Relative humidity saturate[%]' 
   write(17,'(A7)')'ENDVARS'  ! Finalizing the ctl file 

   ! Finalizing the routine
   close(13)
   close(14)   
   close(15)
   close(16)   
   close(17)

   deallocate(UmiRela)
   deallocate(URnegCONT)
   deallocate(URsatCONT)
   deallocate(URnegVAL)
   deallocate(URsatVAL)
   deallocate(URnegV)
   deallocate(URsatV)

   write(stdout,'(A)')'Program hunNphy ends  normaly'
   
end program
