diff --git a/CHANGELOG.md b/CHANGELOG.md index c0a6ece..1604f9e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +28Jan2026: +- bring back changes from PR #209 + +28Jan2026: +- bug fix: add ABI tlapmean back (had lost them when backing out PR#209) +- another miss (qcmod) wrt ABI and IR back out + +12Jan2026: +- bug fixes in interfacing GSI w/ libSP - these are zero-diff fixes + +09Jan2026: +- Back out changes from PR 209 (IR adjustments) since this is + aimed at FPP and the update must be zero-diff. +- optional write out in setupw + 07Jan2026: - minor revisions of CG tolerance (consistency between pcgsoi and bicg) diff --git a/GEOSaana_GridComp/GSI_GridComp/general_specmod.f90 b/GEOSaana_GridComp/GSI_GridComp/general_specmod.f90 index 439e26e..266d3d8 100644 --- a/GEOSaana_GridComp/GSI_GridComp/general_specmod.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/general_specmod.f90 @@ -7,16 +7,16 @@ module general_specmod ! abstract: copy of specmod, introducing structure variable spec_vars, so ! spectral code can be used for arbitrary resolutions. ! -! program history log: +! program history log: ! 2003-11-24 treadon ! 2004-04-28 d. kokron, updated SGI's fft to use scsl ! 2004-05-18 kleist, documentation -! 2004-08-27 treadon - add/initialize variables/arrays needed by -! splib routines for grid <---> spectral +! 2004-08-27 treadon - add/initialize variables/arrays needed by +! splib routines for grid <---> spectral ! transforms ! 2007-04-26 yang - based on idrt value xxxx descriptionxxx ! 2010-02-18 parrish - copy specmod to general_specmod and add structure variable spec_vars. -! remove all *_b variables, since now not necessary to have two +! remove all *_b variables, since now not necessary to have two ! resolutions. any number of resolutions can be now contained in ! type(spec_vars) variables passed in through init_spec_vars. also ! remove init_spec, since not really necessary. @@ -127,7 +127,7 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) ! program history log: ! 2003-11-24 treadon ! 2004-05-18 kleist, new variables and documentation -! 2004-08-27 treadon - add call to sptranf0 and associated arrays, +! 2004-08-27 treadon - add call to sptranf0 and associated arrays, ! remove del21 and other unused arrays/variables ! 2006-04-06 middlecoff - remove jc=ncpus() since not used ! 2008-04-11 safford - rm unused vars @@ -137,7 +137,7 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) ! 2013-10-23 el akkraoui - initialize lats to zero (otherwise point is undefined) ! ! input argument list: -! sp - type(spec_vars) variable +! sp - type(spec_vars) variable ! jcap - target resolution ! jcap_test - test resolution, used to construct mask which will zero out coefs ! with total wavenumber n in range jcap_test < n <= jcap @@ -161,13 +161,15 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) integer(i_kind) ,intent(in ) :: jcap,jcap_test,nlat_a,nlon_a logical,optional,intent(in ) :: eqspace -! Declare local variables +! Declare local variables integer(i_kind) i,ii1,j,l,m,jhe,n integer(i_kind) :: ldafft real(r_kind) :: dlon_a,half_pi,two_pi real(r_kind),dimension(nlat_a-2) :: wlatx,slatx real(r_kind) :: epsi0(0:jcap) ! epsilon factor for m=0 real(r_kind) :: fnum, fden + real(r_kind), dimension(:, :), allocatable :: dummy_w + real(r_kind), dimension(:, :), allocatable :: dummy_g ! Set constants used in transforms for analysis grid sp%jcap=jcap @@ -191,7 +193,6 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) sp%je=(sp%jmax+1)/2 - ! Allocate and initialize fact arrays if(sp%lallocated) then deallocate(sp%factsml,sp%factvml,sp%eps,sp%epstop,sp%enn1,sp%elonn1,sp%eon,sp%eontop) @@ -229,11 +230,17 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) ldafft=50000+4*sp%imax ! ldafft=256+imax would be sufficient at GMAO. allocate( sp%afft(ldafft)) allocate( sp%clat(sp%jb:sp%je) ) - allocate( sp%slat(sp%jb:sp%je) ) - allocate( sp%wlat(sp%jb:sp%je) ) + allocate( sp%slat(sp%jb:sp%je) ) + allocate( sp%wlat(sp%jb:sp%je) ) call spwget(sp%iromb,sp%jcap,sp%eps,sp%epstop,sp%enn1, & sp%elonn1,sp%eon,sp%eontop) - call spffte(sp%imax,(sp%imax+2)/2,sp%imax,2,0.,0.,0,sp%afft) +! Allocate dummy_w and dummy_g arrays for spffte (unused, but required by spffte) + allocate(dummy_w((sp%imax+2)/2,2), dummy_g(sp%imax,2)) + dummy_w=zero + dummy_g=zero + call spffte(sp%imax,(sp%imax+2)/2,sp%imax,2,dummy_w,dummy_g,0,sp%afft) + if(allocated(dummy_w)) deallocate(dummy_w) + if(allocated(dummy_g)) deallocate(dummy_g) call splat(sp%idrt,sp%jmax,slatx,wlatx) jhe=(sp%jmax+1)/2 if(jhe > sp%jmax/2)wlatx(jhe)=wlatx(jhe)/2 @@ -248,14 +255,14 @@ subroutine general_init_spec_vars(sp,jcap,jcap_test,nlat_a,nlon_a,eqspace) allocate( sp%plntop(sp%jcap+1,sp%jb:sp%je) ) do j=sp%jb,sp%je call splegend(sp%iromb,sp%jcap,sp%slat(j),sp%clat(j),sp%eps, & - sp%epstop,sp%pln(1,j),sp%plntop(1,j)) + sp%epstop,sp%pln(:,j),sp%plntop(:,j)) end do else sp%precalc_pln=.false. allocate( sp%pln(sp%ncd2,1) ) allocate( sp%plntop(sp%jcap+1,1) ) end if - + ! obtain rlats and rlons half_pi=half*pi two_pi=two*pi diff --git a/GEOSaana_GridComp/GSI_GridComp/general_transform.f90 b/GEOSaana_GridComp/GSI_GridComp/general_transform.f90 index 8f1af20..2a64a75 100644 --- a/GEOSaana_GridComp/GSI_GridComp/general_transform.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/general_transform.f90 @@ -16,9 +16,9 @@ subroutine general_sptez_s(sp,wave,grid,idir) ! subprogram can be called from a multiprocessing environment. ! ! This routine differs from splib routine sptez in that -! 1) the calling list only contains the in/out arrays and +! 1) the calling list only contains the in/out arrays and ! flag for the direction in which to transform -! 2) it calls a version of sptranf that does not invoke +! 2) it calls a version of sptranf that does not invoke ! initialization routines on each entry ! 3) some generality built into the splib version is ! removed in the code below @@ -142,9 +142,9 @@ subroutine general_sptranf_s(sp_a,wave,grid,idir) ! ! subprograms called: ! sptranf1 sptranf spectral transform -! -! remarks: -! This routine assumes that splib routine sptranf0 has been +! +! remarks: +! This routine assumes that splib routine sptranf0 has been ! previously called. sptranf0 initializes arrays needed in ! the transforms. ! @@ -181,25 +181,27 @@ subroutine general_sptranf_s(sp_a,wave,grid,idir) real(r_kind),dimension(sp_a%imax,2):: g real(r_kind),dimension(sp_a%imax+2,2):: f real(r_kind),dimension(50000+4*sp_a%imax):: tmpafft + integer(i_kind), dimension(1) :: a_mp ! Initialize local variables mp=0 + a_mp=0 do i=1,2*(sp_a%jcap+1) wtop(i)=zero end do ! Transform wave to grid -! ***NOTE*** +! ***NOTE*** ! The FFT used in the transform below has been generalized to -! allow for projection of spectral coefficients onto double -! the desired number of longitudinal grid points. This -! approach is needed when transforming high wavenumber spectral -! coefficients to a coarser resoultion grid. For example, using -! splib to directly transform T878 spectral coefficients to an +! allow for projection of spectral coefficients onto double +! the desired number of longitudinal grid points. This +! approach is needed when transforming high wavenumber spectral +! coefficients to a coarser resoultion grid. For example, using +! splib to directly transform T878 spectral coefficients to an ! 1152 x 576 grid does not use Fourier modes above wavenumber 576. -! Joe Sela insightfully suggested doubling the number of points -! in the FFT and using every other point in the output grid. +! Joe Sela insightfully suggested doubling the number of points +! in the FFT and using every other point in the output grid. ! Mark Iredell coded up Joe's idea below. tmpafft(:)=sp_a%afft(:) @@ -209,7 +211,7 @@ subroutine general_sptranf_s(sp_a,wave,grid,idir) if(idir>0) then do j=sp_a%jb,sp_a%je call spsynth(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,wave,wtop,f) + sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,wave,wtop,f) call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & @@ -231,8 +233,8 @@ subroutine general_sptranf_s(sp_a,wave,grid,idir) ! high spectral representation fields to coarse physical space ! grids. The code below should not be used to transform coarse ! resolution grids to high spectral representation. Since this -! functionality is not yet needed in the GSI, the prudent action -! to take here is to print an ERROR message and terminate program +! functionality is not yet needed in the GSI, the prudent action +! to take here is to print an ERROR message and terminate program ! execution if such a transform is requested. else @@ -248,7 +250,7 @@ subroutine general_sptranf_s(sp_a,wave,grid,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop) + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,f,wave,wtop) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & ! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), & @@ -309,9 +311,9 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) ! ! subprograms called: ! sptranf1 sptranf spectral transform -! -! remarks: -! This routine assumes that splib routine sptranf0 has been +! +! remarks: +! This routine assumes that splib routine sptranf0 has been ! previously called. sptranf0 initializes arrays needed in ! the transforms. ! @@ -345,6 +347,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) ! Declare local variables integer(i_kind) i,j,ii,jj,ijn,ijs,mp,ifact,kw,kwtop,imaxp2 + integer(i_kind), dimension(1) :: a_mp real(r_kind),dimension(2*(sp_b%jcap+1)):: wtop real(r_kind),dimension(sp_b%imax,2):: g real(r_kind),dimension(sp_b%imax+2,2):: f @@ -354,6 +357,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) ! Initialize local variables mp=0 + a_mp=0 ifact = sp_b%imax/sp_a%imax do i=1,2*(sp_b%jcap+1) wtop(i)=zero @@ -363,16 +367,16 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) imaxp2=sp_b%imax+2 ! Transform wave to grid -! ***NOTE*** +! ***NOTE*** ! The FFT used in the transform below has been generalized to -! allow for projection of spectral coefficients onto double -! the desired number of longitudinal grid points. This -! approach is needed when transforming high wavenumber spectral -! coefficients to a coarser resoultion grid. For example, using -! splib to directly transform T878 spectral coefficients to an +! allow for projection of spectral coefficients onto double +! the desired number of longitudinal grid points. This +! approach is needed when transforming high wavenumber spectral +! coefficients to a coarser resoultion grid. For example, using +! splib to directly transform T878 spectral coefficients to an ! 1152 x 576 grid does not use Fourier modes above wavenumber 576. -! Joe Sela insightfully suggested doubling the number of points -! in the FFT and using every other point in the output grid. +! Joe Sela insightfully suggested doubling the number of points +! in the FFT and using every other point in the output grid. ! Mark Iredell coded up Joe's idea below. if(idir>0) then @@ -381,7 +385,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) do j=sp_a%jb,sp_a%je tmpafft(:)=sp_b%afft(:) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),0,wave,wtop,f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,wave,wtop,f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -404,7 +408,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,& sp_b%epstop,tmppln,tmpplntop) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,0,wave,wtop,f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,wave,wtop,f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -428,8 +432,8 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) ! high spectral representation fields to coarse physical space ! grids. The code below should not be used to transform coarse ! resolution grids to high spectral representation. Since this -! functionality is not yet needed in the GSI, the prudent action -! to take here is to print an ERROR message and terminate program +! functionality is not yet needed in the GSI, the prudent action +! to take here is to print an ERROR message and terminate program ! execution if such a transform is requested. else @@ -451,7 +455,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop) + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,f,wave,wtop) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & ! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), & @@ -473,7 +477,7 @@ subroutine general_sptranf_s_b(sp_a,sp_b,wave,grid,idir) sp_a%epstop,tmppln,tmpplntop) call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),0,f,wave,wtop) + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,f,wave,wtop) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & ! tmpafft,sp_a%clat(j),sp_a%slat(j),sp_a%wlat(j), & @@ -516,7 +520,7 @@ subroutine general_sptez_v(sp,waved,wavez,gridu,gridv,idir) ! 1996-02-29 iredell ! 2004-08-23 treadon - adapt splib routine sptezv for gsi use ! 2007-04-25 errico - replace use of duplicate arguments in sptranf_v -! 2008-04-03 safford - rm unused vars +! 2008-04-03 safford - rm unused vars ! 2010-02-18 parrish - copy to general_sptez_v, and pass specmod vars through ! input variable sp of type(spec_vars) ! @@ -649,7 +653,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) ! spdz2uv compute winds from divergence and vorticity ! spuv2dz compute divergence and vorticity from winds ! -! remarks: +! remarks: ! This routine assumes that splib routine sptranf0 has been ! previously called. sptranf0 initializes arrays needed in ! the transforms. @@ -685,7 +689,8 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) ! Declare local variables integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2 - integer(i_kind),dimension(2):: mp + integer(i_kind) :: mp + integer(i_kind), dimension(1) :: a_mp real(r_kind),dimension(sp_b%ncd2*2,2):: w real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop real(r_kind),dimension(sp_b%imax,2):: g @@ -697,6 +702,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) ! Set parameters mp=1 + a_mp=1 ifact = sp_b%imax/sp_a%imax kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2) kwtop=2*(sp_b%jcap+1) @@ -724,7 +730,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) do j=sp_a%jb,sp_a%je tmpafft(:)=sp_b%afft(:) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -740,7 +746,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) gridu(ijs+i)=g(ii,2) enddo call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -760,7 +766,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) call splegend(sp_b%iromb,sp_b%jcap,sp_b%slat(j),sp_b%clat(j),sp_b%eps,& sp_b%epstop,tmppln,tmpplntop) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -776,7 +782,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) gridu(ijs+i)=g(ii,2) enddo call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -821,7 +827,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),mp,f, & + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,f, & w(1,1),wtop(1,1)) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & @@ -834,7 +840,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,j),sp_a%plntop(1,j),mp,f, & + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,j),sp_a%plntop(:,j),a_mp,f, & w(1,2),wtop(1,2)) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & @@ -847,7 +853,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) do j=sp_a%jb,sp_a%je if(sp_a%wlat(j)>zero) then call splegend(sp_a%iromb,sp_a%jcap,sp_a%slat(j),sp_a%clat(j),sp_a%eps,& - sp_a%epstop,sp_a%pln(1,1),sp_a%plntop(1,1)) + sp_a%epstop,sp_a%pln(:,1),sp_a%plntop(:,1)) jj = j-sp_a%jb ijn = jj*sp_a%jn ijs = jj*sp_a%js + sp_a%ioffset @@ -857,7 +863,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,1),sp_a%plntop(1,1),mp,f, & + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,1),sp_a%plntop(:,1),a_mp,f, & w(1,1),wtop(1,1)) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & @@ -870,7 +876,7 @@ subroutine general_sptranf_v(sp_a,sp_b,waved,wavez,gridu,gridv,idir) enddo call spffte(sp_a%imax,imaxp2/2,sp_a%imax,2,f,g,-1,tmpafft) call spanaly(sp_a%iromb,sp_a%jcap,sp_a%imax,imaxp2,kw,kwtop,1, & - sp_a%wlat(j),sp_a%clat(j),sp_a%pln(1,1),sp_a%plntop(1,1),mp,f, & + sp_a%wlat(j),sp_a%clat(j),sp_a%pln(:,1),sp_a%plntop(:,1),a_mp,f, & w(1,2),wtop(1,2)) ! call sptranf1(sp_a%iromb,sp_a%jcap,sp_a%idrt,sp_a%imax,sp_a%jmax,j,j, & ! sp_a%eps,sp_a%epstop,sp_a%enn1,sp_a%elonn1,sp_a%eon,sp_a%eontop, & @@ -945,7 +951,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) ! spdz2uv compute winds from divergence and vorticity ! spuv2dz compute divergence and vorticity from winds ! -! remarks: +! remarks: ! This routine assumes that splib routine sptranf0 has been ! previously called. sptranf0 initializes arrays needed in ! the transforms. @@ -980,7 +986,8 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) ! Declare local variables integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2 - integer(i_kind),dimension(2):: mp + integer(i_kind) :: mp + integer(i_kind), dimension(1) :: a_mp real(r_kind),dimension(sp_b%ncd2*2,2):: w real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop real(r_kind),dimension(sp_b%imax,2):: g @@ -991,6 +998,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) ! Set parameters mp=1 + a_mp=1 ifact = sp_b%imax/sp_a%imax kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2) kwtop=2*(sp_b%jcap+1) @@ -1017,7 +1025,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) do j=sp_a%jb,sp_a%je tmpafft(:)=sp_b%afft(:) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -1034,7 +1042,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) enddo if(j == sp_a%jb)then call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & @@ -1056,7 +1064,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) sp_b%epstop,tmppln,tmpplntop) tmpafft(:)=sp_b%afft(:) call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -1073,7 +1081,7 @@ subroutine general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) enddo if(j == sp_a%jb)then call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & @@ -1145,7 +1153,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) ! spdz2uv compute winds from divergence and vorticity ! spuv2dz compute divergence and vorticity from winds ! -! remarks: +! remarks: ! This routine assumes that splib routine sptranf0 has been ! previously called. sptranf0 initializes arrays needed in ! the transforms. @@ -1180,7 +1188,8 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) ! Declare local variables integer(i_kind) i,j,ii,jj,ijn,ijs,ifact,kw,kwtop,imaxp2 - integer(i_kind),dimension(2):: mp + integer(i_kind) :: mp + integer(i_kind), dimension(1) :: a_mp real(r_kind),dimension(sp_b%ncd2*2,2):: w real(r_kind),dimension(2*(sp_b%jcap+1),2):: wtop real(r_kind),dimension(sp_b%imax,2):: g @@ -1191,6 +1200,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) ! Set parameters mp=1 + a_mp=1 ifact = sp_b%imax/sp_a%imax kw=(sp_b%jcap+1)*((sp_b%iromb+1)*sp_b%jcap+2) kwtop=2*(sp_b%jcap+1) @@ -1220,7 +1230,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) tmpafft(:)=sp_b%afft(:) if(j == sp_a%jb)then call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & @@ -1235,7 +1245,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) enddo end if call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),sp_b%pln(1,j),sp_b%plntop(1,j),mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),sp_b%pln(:,j),sp_b%plntop(:,j),a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & @@ -1260,7 +1270,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) tmpafft(:)=sp_b%afft(:) if(j == sp_a%jb)then call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,1),wtop(1,1),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,1),wtop(1,1),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & @@ -1275,7 +1285,7 @@ subroutine general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) enddo end if call spsynth(sp_b%iromb,sp_b%jcap,sp_b%imax,imaxp2,kw,kwtop,1, & - sp_a%clat(j),tmppln,tmpplntop,mp,w(1,2),wtop(1,2),f) + sp_a%clat(j),tmppln,tmpplntop,a_mp,w(1,2),wtop(1,2),f) call spffte(sp_b%imax,imaxp2/2,sp_b%imax,2,f,g,1,tmpafft) ! call sptranf1(sp_b%iromb,sp_b%jcap,sp_b%idrt,sp_b%imax,sp_a%jmax,j,j, & ! sp_b%eps,sp_b%epstop,sp_b%enn1,sp_b%elonn1,sp_b%eon,sp_b%eontop, & @@ -1311,9 +1321,9 @@ subroutine general_sptez_s_b(sp_a,sp_b,wave,grid,idir) ! subprogram can be called from a multiprocessing environment. ! ! This routine differs from splib routine sptez in that -! 1) the calling list only contains the in/out arrays and +! 1) the calling list only contains the in/out arrays and ! flag for the direction in which to transform -! 2) it calls a version of sptranf that does not invoke +! 2) it calls a version of sptranf that does not invoke ! initialization routines on each entry ! 3) some generality built into the splib version is ! removed in the code below @@ -1415,7 +1425,7 @@ subroutine general_sptez_v_b(sp_a,sp_b,waved,wavez,gridu,gridv,idir,iuvflag) ! 1996-02-29 iredell ! 2004-08-23 treadon - adapt splib routine sptezv for gsi use ! 2007-04-25 errico - replace use of duplicate arguments in sptranf_v -! 2008-04-03 safford - rm unused vars +! 2008-04-03 safford - rm unused vars ! 2010-02-18 parrish - copy to general_sptez_v, and pass specmod vars through ! input variable sp of type(spec_vars) ! @@ -1493,7 +1503,7 @@ subroutine general_sptez_v_b(sp_a,sp_b,waved,wavez,gridu,gridv,idir,iuvflag) if(iuvflag > 0)then ! Call spectral <--> grid transform u only (and polar v) call general_sptranf_v_u(sp_a,sp_b,waved,wavez,gridu,gridv) - else + else ! Call spectral <--> grid transform v only (and polar u) call general_sptranf_v_v(sp_a,sp_b,waved,wavez,gridu,gridv) end if diff --git a/GEOSaana_GridComp/GSI_GridComp/setupw.f90 b/GEOSaana_GridComp/GSI_GridComp/setupw.f90 index 0246155..5925111 100644 --- a/GEOSaana_GridComp/GSI_GridComp/setupw.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/setupw.f90 @@ -423,7 +423,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ikx=nint(data(ikxx,i)) if(ikx < 1 .or. ikx > nconvtype) then num_bad_ikx=num_bad_ikx+1 - if(num_bad_ikx<=10) write(6,*)' in setupw ',ikx,i,nconvtype,mype + if(verbose_hires_raob.and.num_bad_ikx<=10) write(6,*)' in setupw ',ikx,i,nconvtype,mype end if end do ! If HD raobs available move prepbufr version to monitor diff --git a/GEOSaana_GridComp/GSI_GridComp/strong_fast_global_mod.f90 b/GEOSaana_GridComp/GSI_GridComp/strong_fast_global_mod.f90 index 9e0834d..a98dd02 100644 --- a/GEOSaana_GridComp/GSI_GridComp/strong_fast_global_mod.f90 +++ b/GEOSaana_GridComp/GSI_GridComp/strong_fast_global_mod.f90 @@ -79,7 +79,7 @@ module strong_fast_global_mod integer(i_kind),allocatable:: mode_list(:,:) ! mode_list(1,j) = lat index for ew strip j ! mode_list(2,j) = vert mode number for ew strip j - ! mode_list(3,j) = pe of this lat/vert mode strip + ! mode_list(3,j) = pe of this lat/vert mode strip integer(i_kind),allocatable:: mmode_list(:,:) ! mmode_list(1,j) = m1 (zonal wave number 1) for ns strip ! mmode_list(2,j) = m2 (zonal wave number 2) for ns strip @@ -179,7 +179,7 @@ subroutine strong_bal_correction_fast_global(u_t,v_t,t_t,ps_t,mype,psi,chi,t,ps, ! bal_diagnostic - if true, then compute bal diagnostic, a measure of amplitude ! of balanced gravity mode tendencies ! fullfield - if true, full field diagnostics -! if false, incremental +! if false, incremental ! update - if false, then do not update u,v,t,ps with balance increment ! ! output argument list: @@ -300,7 +300,7 @@ subroutine strong_bal_correction_fast_global(u_t,v_t,t_t,ps_t,mype,psi,chi,t,ps, deldivhat(2,n)=deldivhat(2,n)*del2inv end do end if - + i=0 do n=m,sp_a%jcap i=i+1 @@ -377,7 +377,7 @@ subroutine strong_bal_correction_fast_global(u_t,v_t,t_t,ps_t,mype,psi,chi,t,ps, do i=1,nvmodes_keep rmstend_all_uf=rmstend_all_uf+rmstend_uf(i) rmstend_all_g_uf=rmstend_all_g_uf+rmstend_g_uf(i) - diffi = rmstend_uf(i)-rmstend_g_uf(i) + diffi = rmstend_uf(i)-rmstend_g_uf(i) diff1 = rmstend_uf(1)-rmstend_g_uf(1) if(abs(diffi)