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