!-----------------------------------------------------------------------------! ! 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