!----------------------------------------------------------------------!
!                        BAM/CPTEC/INPE - 2017                         !
!----------------------------------------------------------------------!
!BOP
!
! !DESCRIPTION:
!      Programa para abrir e imprimir os valores do arquivo 
!      global_berror.l64y386.f77.gcv ou de qualquer outra matriz no 
!      mesmo formato.
!      Vide https://projetos.cptec.inpe.br/projects/berror/repository/show/bins
!      para uma lista de matrizes disponiveis em diferentes resolucoes
!
! !INTERFACE:
!      gfortran -fconvert=big-endian adjBerr.f90 -o adjBerr.exe
!
! !REVISION HISTORY:
!      01-06-2015 - Carlos Frederico Bastarz - Codigo Inicial
!      02-06-2015 - Carlos Frederico Bastarz - Adicao de subrotinas para 
!                   verificacao da leitura e tamanhos dos records
!      12-06-2015 - Correcao de algumas informacoes sobre os records
!      26-06-2017 - Inclusao deste programa no pacote DIAG do EVAL
!      27-06-2017 - Criacao dos arquivos CTLs para as amplitudes 
!                   de acordo com as dimensoes da matriz lida.
!
! !REMARKS:
!      Esta sendo utilizada a subrotina gauss2lats para calcular
!      as latitudes para o arquivo CTL.
!      A versao atual deste programa esta escrevendo os CTLs
!      apenas para as amplitudes (desvios-padrao) das quantidades
!      representadas.
!
! !TODO:
!      Remover a subrotina gauss2ctl para um modulo;
!      Generalizar o calculo dos niveis sigma;
!      Incluir escrita arquivo CTL para matrizes de projecao e
!      comprimentos de escala.
!
! !BUGS:
!      Nenhum ate o momento.
!
!----------------------------------------------------------------------!
!EOP

program read_gsi_berror

  implicit none

  integer :: i, isig, inerr,   &  ! i: contador 
             nlat, nlon, nsig, &  ! inerr: codigo de status de leitura (0=ok) 
             rlen                 ! nlat/nlon/nsig: numero de latitudes, longitudes e niveis sigma da matriz B
                                  ! rlen: dimensao do record (record length) 

  integer, parameter :: berror=10 ! unidade de leitura da matriz B

  ! Classificacao dos Arrays:
  ! Tipo                                                 Arrays
  ! ===================================================  ===============================
  ! 1) Raiz quadrada da variancia (desvio padrao)        corpin, corq2in, corsst, corzin
  ! 2) Comprimento de Escala Horizontal (da correlacao)  hscalespin, hsstin, hscalesin
  ! 3) Comprimento de Escala Vertical (da correlacao)    vscalesin
  ! 4) Coeficiente de Regressao (da variavel)            wgvin, bgvin, agvin

  ! Classificacao das Variaveis:
  ! Indice  Variavel                       Dimensao             Arrays
  ! ======  =============================  ===================  =====================================
  ! i=1     Stream Function (sf)           3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin
  ! i=2     Velocity Potential (vp)        3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin
  ! i=3     Absolute Temperature (t)       3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin
  ! i=4     Relative Humidity (q)          3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin, corq2in
  ! i=5     Ozone (oz)                     3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin
  ! i=6     Cloud Water (cw)               3D (nlat,nsig,nsig)  corzin, hscalesin, vscalesin      
  ! i=7     Surface Pressure (ps)          2D (nlat,nsig)       corpin, hscalesin
  ! i=8     Sea Surface Temperature (sst)  2D (nlat,nsig)       hsstin, corsstin

  ! Arrays da matriz B
  ! - Arrays dependentes de nlat:
  real(4), allocatable, dimension(:)     :: corpin, hscalespin   ! corpin: raiz quadrada da variancia da ps desbalanceada
                                                                 ! hscalespin: comprimento de escala horizontal da ps desbalanceada
  ! - Arrays dependentes de nlat, nsig
  real(4), allocatable, dimension(:,:)   :: corq2in,      &      ! corq2in: raiz quadrada da variancia da q
                                            wgvin, bgvin, &      ! wgvin: coeficientes de regressao da funcao da sp e ps (sf->ps)
                                            corsstin, hsstin     ! corsstin: raiz quadrada da variancia da sst
                                                                 ! hsstin: comprimento de escala horizontal da sst
                                                                 ! bgvin: coeficientes de regressao da sf e vp (sf->vp)
  ! - Arrays dependentes de nlat, nsig e nsig
  real(4), allocatable, dimension(:,:,:) :: agvin, corzin, &     ! corzin: raiz quadrada da variancia da sf, vp desbalanceada, t desbalanceada e q
                                            hscalesin, vscalesin ! hscalesin: comprimento de escala horizontal da sf, vp desbalanceada, t desbalanceada e q
                                                                 ! vscalesin: comprimento de escala vertical da sf, vp desbalanceada, t desbalanceada e q
                                                                 ! agvin: coeficientes de regressao da sf e t (sf->t)

  character(5)   :: var(40)       ! var: nome da variavel lida

  character(100) :: fmt1, fmt2    ! fmt1/fmt2: formatos para escrita na tela (apenas)

  character(1)   :: hflagin, &    ! hflagin: flag lida da linha de comando para escrever (T) ou nao (F) os records da matriz
                    hflagon, T, F ! hflagon: flag utilizada para atribuir os valores logicos T ou F para a variavel hflag \
                                  ! T: valor logico true \
                                  ! F: valor logico false

  character(100) :: filein        ! filein: nome do arquivo a ser lido

  logical :: hflag                ! hflag: flag utilizada dentro do programa para escrever ou nao os records da matriz

  integer :: ninput, iargc        ! ninput: numero de argumentos lidos da linha de comando

! -- Inicio da leitura da matriz B

  ! Formatos
  fmt1="(I2,A10,I2)"
  fmt2="(100A)"

  ! Leitura dos argumentos da linha de comando
  ninput=iargc()

  if(ninput.ne.2)then

    print *, "Input arguments: filein flag"
    print *, "filein: Berror file name"
    print *, "flag..: Flag to write output files (T) or just to dump general information (F)"
    print *, ""
    stop "Error reading input arguments"

  else

    call getarg(1,filein)
    call getarg(2,hflagin)

    ! Flag para dump do header
    read(hflagin,*) hflagon
    
    if(hflagon.eq.'T')then
      hflag=.TRUE.
    else if(hflagon.eq.'F')then
      hflag=.FALSE.
    else
      stop "hflag must be either T or F"
    end if

    ! Abre a matriz B para leitura
    open(berror,file=filein,form='unformatted',status='old')
  
    ! Posiciona a leitura para o comeco do arquivo
    rewind(berror)

    ! Le a primeira linha do arquivo
    read(berror,iostat=inerr) nsig, nlat, nlon
    call check_inerr(inerr,"nsig,nlat,nlon")

    ! Imprime as informacoes lidas
    write(*,fmt2) "- Grid Dimensions (nsig, nlat, nlon):"
    print *, nsig, nlat, nlon
  
    ! Aloca as matrizes a serem lidas
    ! corzin, hscalesin e vscales in possuem zdim=6 porque as variavels ps e sst nao \
    ! fazem parte deste record
    allocate(agvin(nlat,nsig,nsig))
    allocate(wgvin(nlat,nsig),bgvin(nlat,nsig))
    allocate(corzin(nlat,nsig,6),corpin(nlat),corq2in(nlat,nsig))
    allocate(hscalesin(nlat,nsig,6),hscalespin(nlat))
    allocate(vscalesin(nlat,nsig,6))
    allocate(corsstin(nlat,nlon),hsstin(nlat,nlon))
  
    print *
  
    ! Reposiciona a leitura para o inicio do arquivo
    rewind(berror)

    ! Le novamente as primeira linha
    read(berror,iostat=inerr) nsig, nlat, nlon
    call check_inerr(inerr,"nsig,nlat,nlon")

    ! Le as amplitudes agvin, bgvin e wgvin
    write(*,fmt2) "- Error Amplitudes:"
    read(berror,iostat=inerr) agvin, bgvin, wgvin ! E razoavel que alguns dos valores destes arrays contenham valores negativos
    call save_3d(agvin,nlat,nsig,nsig,"","agvin",hflag)
    call save_2d(bgvin,nlat,nsig,"","bgvin",hflag)
    call save_2d(wgvin,nlat,nsig,"","wgvin",hflag)

    ! Leitura dos records das variaveis, na seguinte ordem:
    ! i=1: Stream Function
    ! i=2: Velocity Potential
    ! i=3: Absolute Temperature
    ! i=4: Specific Humidity
    ! i=5: Ozone
    ! i=6: Cloud Water
    write(*,fmt2) "- Correlations and Lenght Scales:"
    do i=1,6 
      call print_info(i,var(i),nsig)

      read(berror,iostat=inerr) var(i),nsig ! Na matriz B do NCEP, estas informacoes estao em branco, na matriz do CPTEC, temos o nome da variavel lida
      call check_inerr(inerr,"var,nsig")

      if (i==4) then ! Quando a variavel for i=4 (Specific Humidity), le tambem o array corq2in (que tambem depende apenas de nlat e nlon)
        read(berror,iostat=inerr) corzin(:,:,i),corq2in
        call check_inerr(inerr,"corzin,corq2in")
        call save_3d(corzin(:,:,i),nlat,nsig,i,var(i),"corzin",hflag)

        call create_ctl(i,nlat,nlon,nsig,var(i),"corzin",hflag)

        call save_2d(corq2in,nlat,nsig,var(i),"corqin",hflag)

        read(berror,iostat=inerr) hscalesin(:,:,i)
        call check_inerr(inerr,"hscalesin")
        call save_3d(hscalesin(:,:,i),nlat,nsig,i,var(i),"hscalesin",hflag)

        read(berror,iostat=inerr) vscalesin(:,:,i)
        call check_inerr(inerr,"vscalesin")
        call save_3d(vscalesin(:,:,i),nlat,nsig,i,var(i),"vscalesin",hflag)
      else
        read(berror,iostat=inerr) corzin(:,:,i)
        call check_inerr(inerr,"corzin")
        call save_3d(corzin(:,:,i),nlat,nsig,i,var(i),"corzin",hflag)

        call create_ctl(i,nlat,nlon,nsig,var(i),"corzin",hflag)

        read(berror,iostat=inerr) hscalesin(:,:,i)
        call check_inerr(inerr,"hscalesin")
        call save_3d(hscalesin(:,:,i),nlat,nsig,i,var(i),"hscalesin",hflag)

        read(berror,iostat=inerr) vscalesin(:,:,i)
        call check_inerr(inerr,"vscalesin")
        call save_3d(vscalesin(:,:,i),nlat,nsig,i,var(i),"vscalesin",hflag)

      end if
    end do
  
    isig=1
  
    ! i=7: Surface Pressure
    i=7
    read(berror,iostat=inerr) var(7),isig
    call check_inerr(inerr,"var,isig")
    call print_info(i,var(i),isig)

    read(berror,iostat=inerr) corpin
    call check_inerr(inerr,"corpin")
    call save_1d(corpin,nlat,var(i),"corpin",hflag)

    call create_ctl(i,nlat,nlon,1,var(i),"corpin",hflag)

    read(berror,iostat=inerr) hscalespin
    call check_inerr(inerr,"hscalespin")
    call save_1d(hscalespin,nlat,var(i),"hscalespin",hflag)
  
    ! i=8: Sea Surface Temperature
    i=8
    read(berror,iostat=inerr) var(8),isig
    call check_inerr(inerr,"var,isig")
    call print_info(i,var(i),isig)

    read(berror,iostat=inerr) corsstin
    call check_inerr(inerr,"corsstin")
    call save_2d(corsstin,nlat,nsig,var(8),"corsstin",hflag)

    call create_ctl(i,nlat,nlon,1,var(i),"corsstin",hflag)

    read(berror,iostat=inerr) hsstin
    call check_inerr(inerr,"hsstin")
    call save_2d(hsstin,nlat,nsig,var(8),"hsstin",hflag)
  
    ! Fecha o arquivo
    close(berror)
  
  ! -- Fim da leitura da matriz B
  
    ! Desaloca todas as matrizes alocadas
    deallocate(agvin,wgvin,bgvin)
    deallocate(corzin,corpin,corq2in)
    deallocate(hscalesin,vscalesin,hscalespin)
    deallocate(corsstin,hsstin)

  end if

end program read_gsi_berror

subroutine check_inerr(code,varin)

  integer, intent(in) :: code

  character(*), intent(in) :: varin

  if(code.ne.0)then

    print *, "Record: ", varin
    stop "Erro na leitura!"

  end if

end subroutine check_inerr

subroutine print_info(i,varin,zdim)

  integer, intent(in) :: zdim

  integer, intent(inout) :: i

  character(25) :: varin(25), varout(25)

  character(100) :: fmt

  fmt="(I1,A,1x,2A,I2,1x,A)"

  if(i.eq.1)then; varout="Stream Function"; endif
  if(i.eq.2)then; varout="Velocity Potential"; endif
  if(i.eq.3)then; varout="Temperature"; endif
  if(i.eq.4)then; varout="Specific Humidity"; endif
  if(i.eq.5)then; varout="Ozone"; endif
  if(i.eq.6)then; varout="Cloud Water"; endif
  if(i.eq.7)then; varout="Surface Pressure"; endif
  if(i.eq.8)then; varout="Sea Surface Temperature"; endif
  
  write(*,fmt) i,".",trim(varout(i))," - ",zdim,"level(s)"

end subroutine print_info

subroutine save_1d(work,xdim,varin,fileout,flag)

  integer :: rdim 

  integer, intent(in) :: xdim

  real(4), dimension(xdim), intent(in) :: work

  character(*), intent(in) :: fileout, varin

  character(100) :: fmt1, fmt2, fmt3, record, filename

  logical :: flag

  fmt1="(100A)"
  fmt2="(A,F20.10)"
  fmt3="(A7,I18)"

  inquire(iolength=rdim) work

  if(len(varin).eq.0)then
    record=trim(fileout)
    filename=trim(fileout)//".dat"
  else
    record=trim(fileout)//"-"//trim(varin)
    filename=trim(fileout)//"-"//trim(varin)//".dat"
  end if

  write(*,fmt1) "Record: ", trim(record)
  write(*,fmt3) "Length: ", rdim
  write(*,fmt2) "max: ", maxval(work)
  write(*,fmt2) "min: ", minval(work)

  if(flag)then
    write(*,fmt1) "Output: ", trim(filename)
    open(100,file=filename,form='unformatted',status='unknown')
    write(100) work
    close(100)
  end if

  print *

end subroutine save_1d

subroutine save_2d(work,xdim,ydim,varin,fileout,flag)

  integer :: rdim

  integer, intent(in) :: xdim, ydim

  real(4), dimension(xdim,ydim), intent(in) :: work

  character(*), intent(in) :: fileout, varin

  character(100) :: fmt1, fmt2, fmt3, record, filename

  logical :: flag

  fmt1="(100A)"
  fmt2="(A,10F20.10)"
  fmt3="(A7,I18)"

  inquire(iolength=rdim) work

  if(len(varin).eq.0)then
    record=trim(fileout)
    filename=trim(fileout)//".dat"
  else
    record=trim(fileout)//"-"//trim(varin)
    filename=trim(fileout)//"-"//trim(varin)//".dat"
  end if

  write(*,fmt1) "Record: ", trim(record)
  write(*,fmt3) "Length: ", rdim
  write(*,fmt2) "max: ", maxval(work)
  write(*,fmt2) "min: ", minval(work)

  if(flag)then
    write(*,fmt1) "Output: ", trim(filename)
    open(200,file=filename,form='unformatted',status='unknown')
    write(200) work
    close(200)
  end if

  print *

end subroutine save_2d

subroutine save_3d(work,xdim,ydim,zdim,varin,fileout,flag)

  integer :: rdim

  integer, intent(in) :: xdim, ydim, zdim

  real(4), dimension(xdim,ydim,zdim), intent(in) :: work

  character(*), intent(in) :: fileout, varin

  character(100) :: fmt1, fmt2, fmt3, record, filename

  logical :: flag

  fmt1="(100A)"
  fmt2="(A,10F20.10)"
  fmt3="(A7,I18)"

  inquire(iolength=rdim) work

  if(len(varin).eq.0)then
    record=trim(fileout)
    filename=trim(fileout)//".dat"
  else
    record=trim(fileout)//"-"//trim(varin)
    filename=trim(fileout)//"-"//trim(varin)//".dat"
  end if

  write(*,fmt1) "Record: ", trim(record)
  write(*,fmt3) "Length: ", rdim
  write(*,fmt2) "max: ", maxval(work)
  write(*,fmt2) "min: ", minval(work)

  if(flag)then
    write(*,fmt1) "Output: ", trim(filename)
    open(300,file=filename,form='unformatted',status='unknown')
    write(300) work
    close(300)
  end if

  print *

end subroutine save_3d

subroutine ctl_info(i,varout)

  integer, intent(in) :: i

  character(25), intent(out) :: varout

  character(100) :: fmt

  fmt="(I1,A,1x,2A,I2,1x,A)"

  if(i.eq.1)then; varout="Stream Function"; endif
  if(i.eq.2)then; varout="Velocity Potential"; endif
  if(i.eq.3)then; varout="Temperature"; endif
  if(i.eq.4)then; varout="Specific Humidity"; endif
  if(i.eq.5)then; varout="Ozone"; endif
  if(i.eq.6)then; varout="Cloud Water"; endif
  if(i.eq.7)then; varout="Surface Pressure"; endif
  if(i.eq.8)then; varout="Sea Surface Temperature"; endif
  
end subroutine ctl_info

subroutine create_ctl(i,ydim,xdim,zdim,varin,varname,flag)

  integer, intent(in) :: i, xdim, ydim, zdim

  character(*), intent(in) :: varname, varin

  character(25) :: varout

  real(4), allocatable, dimension(:) :: lats
  real(4), allocatable, dimension(:) :: sigs

  character(100) :: fmt1, record, filenameCTL, filenameDAT

  logical :: flag

  fmt1="(100A)"

  if(len(varin).eq.0)then
    record=trim(varname)
    filenameCTL=trim(varname)//".ctl"
    filenameDAT=trim(varname)//".dat"
  else
    record=trim(varname)//"-"//trim(varin)
    filenameCTL=trim(varname)//"-"//trim(varin)//".ctl"
    filenameDAT=trim(varname)//"-"//trim(varin)//".dat"
  end if

  if(flag)then

    allocate(lats(ydim),sigs(zdim))

    lats = 0.
    sigs = 0.

    call lat_sig(varin,ydim,zdim,lats,sigs)
    call ctl_info(i,varout)

    write(*,fmt1) "Output: ", trim(filenameCTL)
    open(400,file=filenameCTL,form='formatted',status='replace')
    write(400,100) filenameDAT
    if(varin.ne."ps".and.varin.ne."sst")then
    write(400,105) varout
    else if(varin.eq."ps".or.varin.eq."sst")then
    write(400,110) varout
    else
    write(400,115) varout
    end if
    write(400,120)
    write(400,125)
    write(400,130) xdim/360.
    write(400,135) ydim
    write(400,140) lats
    if(varin.eq."ps".or.varin.eq."sst")then
    write(400,145) zdim,sigs
    else
    write(400,150) zdim
    write(400,155) sigs 
    end if
    write(400,160)
    write(400,165)
    write(400,170) varname,zdim,varout
    write(400,175)

    close(400)

100 format ("dset ^",A)
    if(varname.eq."corq2in".or.varname.eq."corzin")then
105 format ("title Sigma Level x Latitude Standard Deviation Distribution for ",A)
    else if(varname.eq."corpin".or.varname.eq."corsstin")then
110 format ("title Latitudinal Standard Deviation Distribution for ",A)
    else
115 format ("title Standard Deviation for ",A)
    end if
120 format ("undef -9.99e33")
125 format ("options big_endian")
130 format ("xdef 1 linear 0 ",F6.4,"")
135 format ("ydef ",I6," levels")
140 format (8F11.6)
    if(varin.eq."ps".or.varin.eq."sst")then
145 format ("zdef ",I6," levels",8F11.6)
    else
150 format ("zdef ",I6," levels")
155 format (8F11.6)
    end if
160 format ("tdef 1 linear 00z01jan2013 6hr")
165 format ("vars 1")
    if(varname.eq."corpin".or.varname.eq."corq2in".or. &
       varname.eq."corsstin".or.varname.eq."corzin")then
170 format (A," ",I6," 99 Standard Deviation for ",A)
    end if
175 format ("endvars")

    deallocate(lats,sigs)

  end if

  print *

end subroutine create_ctl

subroutine lat_sig(varin,ydim,zdim,lats,sigs)

  implicit none

  integer :: j

  real(4) :: lat, latr, pi

  integer, intent(in) :: ydim, zdim

  real(4), dimension(zdim), intent(inout) :: sigs

  real(4), dimension(ydim), intent(inout) :: lats
  real(4), dimension(ydim) :: latweight

  character(*), intent(in) :: varin

  pi=4.0d0*datan(1.0d0)

  call gauss2lats(ydim,lats)

  do j=1, ydim

    lat=lats(j)

    latr=(pi*lat)/180.0
    latweight(j)=cos(latr)

    if (latweight(j) .lt. 0.0) latweight(j) = 0.0

  end do

  ! Os valores abaixo foram calculados a partir do Chopping_serial do
  ! pre-processamento do BAM, utilizando os arquivos DeltaSigma.XXX.

  if(varin.eq."ps".or.varin.eq."sst")then
    sigs = (/ 1 /)
  else
    if(zdim.eq.9)then
      sigs = (/ 0.994997, 0.981491, 0.946550, 0.850494, 0.654281, &
                0.393673, 0.178303, 0.062615, 0.012449 /)
    elseif(zdim.eq.18)then
      sigs = (/ 0.994997, 0.981491, 0.959472, 0.925558, 0.875610, &  
                0.805666, 0.714115, 0.603990, 0.483907, 0.366042, &
                0.261754, 0.177767, 0.115254, 0.071490, 0.042190, &
                0.023171, 0.011284, 0.002905 /)
    else if(zdim.eq.28)then
      sigs = (/ 0.994997, 0.982082, 0.964373, 0.942547, 0.915917, &  
                0.883844, 0.845785, 0.801416, 0.750756, 0.694265, & 
                0.632896, 0.568091, 0.501681, 0.435679, 0.372048, &  
                0.312479, 0.258231, 0.210057, 0.168230, 0.132611, & 
                0.102782, 0.078152, 0.058048, 0.041795, 0.028750, & 
                0.018336, 0.010056, 0.002726 /)
    else if(zdim.eq.42)then 
      sigs = (/ 0.995983, 0.987352, 0.977447, 0.966105, 0.953154, &   
                0.938427, 0.921760, 0.902982, 0.881943, 0.858509, & 
                0.832589, 0.804137, 0.773155, 0.739722, 0.703998, & 
                0.666213, 0.626683, 0.585798, 0.544008, 0.501804, &  
                0.459701, 0.418210, 0.377801, 0.338903, 0.301882, & 
                0.267023, 0.234526, 0.204515, 0.177045, 0.152100, & 
                0.129606, 0.109457, 0.091512, 0.075612, 0.061597, &
                0.049290, 0.038517, 0.029117, 0.020938, 0.013831, &
                0.007647, 0.002042 /)
    else if(zdim.eq.64)then
      sigs = (/ 0.997335, 0.991650, 0.985215, 0.977940, 0.969732, & 
                0.960491, 0.950113, 0.938491, 0.925513, 0.911071, & 
                0.895062, 0.877389, 0.857970, 0.836740, 0.813660, &
                0.788719, 0.761942, 0.733396, 0.703189, 0.671477, &
                0.638458, 0.604374, 0.569502, 0.534148, 0.498634, & 
                0.463287, 0.428431, 0.394369, 0.361378, 0.329698, &
                0.299528, 0.271022, 0.244288, 0.219390, 0.196353, &
                0.175166, 0.155788, 0.138155, 0.122183, 0.107777, &
                0.094832, 0.083239, 0.072890, 0.063674, 0.055489, &
                0.048234, 0.041817, 0.036150, 0.031152, 0.026750, &
                0.022879, 0.019477, 0.016490, 0.013870, 0.011573, &
                0.009560, 0.007798, 0.006255, 0.004904, 0.003724, &
                0.002691, 0.001787, 0.000994, 0.000266 /)
    else if(zdim.eq.96)then
      sigs = (/ 0.997335, 0.991650, 0.985411, 0.978752, 0.971651, &  
                0.964096, 0.956069, 0.947550, 0.938518, 0.928955, &  
                0.918845, 0.908171, 0.896918, 0.885075, 0.872630, & 
                0.859578, 0.845912, 0.831632, 0.816741, 0.801245, &  
                0.785154, 0.768483, 0.751253, 0.733487, 0.715215, & 
                0.696470, 0.677291, 0.657722, 0.637808, 0.617600, &
                0.597153, 0.576523, 0.555768, 0.534949, 0.514127, &
                0.493364, 0.472717, 0.452248, 0.432012, 0.412065, & 
                0.392456, 0.373235, 0.354443, 0.336121, 0.318302, &
                0.301016, 0.284289, 0.268141, 0.252586, 0.237637, &
                0.223299, 0.209577, 0.196469, 0.183971, 0.172076, & 
                0.160774, 0.150052, 0.139897, 0.130292, 0.121219, &
                0.112662, 0.104599, 0.097011, 0.089879, 0.083180, &
                0.076896, 0.071006, 0.065490, 0.060328, 0.055501, &
                0.050990, 0.046777, 0.042846, 0.039179, 0.035760, &
                0.032575, 0.029608, 0.026846, 0.024275, 0.021884, &
                0.019661, 0.017593, 0.015672, 0.013887, 0.012229, & 
                0.010689, 0.009260, 0.007933, 0.006701, 0.005559, & 
                0.004498, 0.003514, 0.002602, 0.001764, 0.000994, & 
                0.000266 /)
    else
      stop "Wrong number of sigma levels."
    end if
  end if


end subroutine lat_sig

subroutine gauss2lats(nlat,xlats)

  implicit none

  integer, intent(in) :: nlat

  integer :: i, nzero

  real :: acon, fi, fi1, a, b, c, d, ylat, gord, delta, pi

  real, parameter :: xlim=1.0e-6 ! Convervence criterion for iteration of cos latitude

  real, dimension(:), allocatable :: cosc, gwt, sinc, colat, wos2

  real, dimension(:), allocatable :: dlat, g, gm, gp, gt

  real, dimension(nlat/2) :: xlat, xlatm

  real, dimension(nlat) :: xlats

  pi=4.0d0*datan(1.0d0)

  acon=180.0d0/pi

  allocate(cosc(nlat+1))
  allocate(sinc(nlat+1))
  allocate(colat(nlat+1))

  allocate(gwt(nlat))
  allocate(wos2(nlat))

  ! Initialize arrays
  cosc=0.0d0
  gwt=0.0d0
  sinc=0.0d0
  colat=0.0d0
  wos2=0.0d0

  ! The number of zeros between pole and Equator
  nzero=nlat/2

  allocate(dlat(nzero))

  ! Set the first guess for cos(colat)
  do i=1, nzero
    cosc(i)=sin((i-0.5d0)*pi/nlat+pi*0.5d0)
  end do

  ! Constants for determining the derivative of the polynomial
  fi=nlat
  fi1=fi+1.0d0
  a=fi*fi1/sqrt(4.0d0*fi1*fi1-1.0d0)
  b=fi1*fi/sqrt(4.0d0*fi*fi-1.0d0)

  allocate(g(nlat))
  allocate(gm(nlat-1))
  allocate(gt(nlat+1))
  allocate(gp(nlat+1))

  ! Loop over latitudes, iterating the search for each root
  do i=1, nzero
    ! Determine the value of the ordinary Legendre polynomial for the current guess root
    g=gord(nlat,cosc(i))
    ! Determine the derivative of the polynomial at this point
    gm=gord(nlat-1,cosc(i))
    gp=gord(nlat+1,cosc(i)) 
    gt=(cosc(i)*cosc(i)-1.0d0)/(a*gp-b*gm)
    ! Update the estimate of the root
    delta=g(i)*gt(i)
    cosc(i)=cosc(i)-delta
    ! If convergence criterion has not been met, keep trying
    do while(abs(delta).gt.xlim)
      g=gord(nlat,cosc(i))
      gm=gord(nlat-1,cosc(i))
      gp=gord(nlat+1,cosc(i))
      gt=(cosc(i)*cosc(i)-1)/(a*gp-b*gm)
      delta=g(i)*gt(i)
      cosc(i)=cosc(i)-delta
    end do
    ! Determine the Gaussian weights
    c=2.0d0*(1.0d0-cosc(i)*cosc(i))
    d=gord(nlat-1,cosc(i))
    d=d*d*fi*fi
    gwt(i)=c*(fi-0.5d0)/d
  end do

  ! Determine the colatitudes and sin(colat) and weights over sin
  do i=1,nzero
    colat(i)=acos(cosc(i))
    sinc(i)=sin(colat(i))
    wos2(i)=gwt(i)/(sinc(i)*sinc(i))
  end do
  
  ! If nlat is odd, set values at the equator
  if (mod(nlat,2).ne.0) then
    i=nzero+1
    cosc(i)=0.0d0
    c=2.0d0
    d=gord(nlat-1,cosc(i))
    d=d*d*fi*fi
    gwt(i)=c*(fi-0.5d0)/d
    colat(i)=pi*0.5d0
    sinc(i)=1.0d0
    wos2(i)=gwt(i)
  end if
  
  ! Determine the southern hemisphere values by symmetry
!  do i=nlat-nzero+1,nlat
!    cosc(i)=-cosc(nlat+1-i)
!    gwt(i)=gwt(nlat+1-i)
!    colat(i)=pi-colat(nlat+1-i)
!    sinc(i)=sinc(nlat+1-i)
!    wos2(i)=wos2(nlat+1-i)
!  end do

  ylat=90.0d0

  ! Calculate latitudes and latitude spacing
  do i=1,nzero
    xlat(i)=acos(sinc(i))*acon
    dlat(i)=xlat(i)-ylat
    ylat=xlat(i)
  end do

!  print *, 'xlat=', xlat

  ! Reverse latitudes for south hemisphere and assembly full array 
  xlatm=-1.0d0*xlat

!  print *, 'xlatm=', xlatm

  xlats=0.0d0

  do i=1,nzero
    xlats(i)=xlatm(i)
  end do

  xlat=xlat(nzero:1:-1)

  do i=nzero+1,nlat
    xlats(i)=xlat(i-nzero)
  end do

!  print *, 'xlats=', xlats

  deallocate(cosc)
  deallocate(gwt)
  deallocate(sinc)
  deallocate(colat)
  deallocate(wos2)

  deallocate(g)
  deallocate(gm)
  deallocate(gt)
  deallocate(gp)
  deallocate(dlat)

end subroutine gauss2lats

real function gord(n,x)

  ! Determine the colatitude
  colat=acos(x)
  
  c1=sqrt(2.0d0)
  
  do i=1,n
    c1=c1*sqrt(1.0d0-1.0d0/(4*i*i))
  end do
  
  fn=n
  ang=fn*colat
  s1=0.0d0
  c4=1.0d0
  a=-1.0d0
  b=0.0d0
  
  do k=0,n,2
    if (k.eq.n) then
      c4=0.5d0*c4
    end if
    s1=s1+c4*cos(ang)
    a=a+2.0d0
    b=b+1.0d0
    fk=k
    ang=colat*(fn-fk-2.0d0)
    c4=(a*(fn-b+1.0d0)/(b*(fn+fn-a)))*c4
  end do

  gord=s1*c1

  return

end function gord
