diff --git a/src/gsi/lag_traj.f90 b/src/gsi/lag_traj.f90 index ca21a141ca..80160b3f50 100644 --- a/src/gsi/lag_traj.f90 +++ b/src/gsi/lag_traj.f90 @@ -7,10 +7,10 @@ module lag_traj ! abstract: This module contain, the trajectory model used to forecast ! the position of superpressure balloon when used in data ! assimilation. -! - A model based on a Runge-Kutta (4th order) +! - A model based on a runge-kutta (4th order) ! algorithm is fully implemented (non-linear, linear and ! adjoint). -! - A model based on a Runge-Kutta (2nd order) +! - A model based on a runge-kutta (2nd order) ! algorithm is fully implemented (non-linear, linear and ! adjoint). ! @@ -61,11 +61,11 @@ module lag_traj public:: lag_rk2iter_nl,lag_rk2iter_tl,lag_rk2iter_ad ! Number of parameters needed by the tangent linear model, for one step - ! of the Runge-Kutta (4th order) trajectory model + ! of the runge-kutta (4th order) trajectory model integer(i_kind),parameter::lag_rk4stepnpara_r=57 ! real numbers integer(i_kind),parameter::lag_rk4stepnpara_i=32 ! integer numbers ! Number of parameters needed by the tangent linear model, for one step - ! of the Runge-Kutta (2nd order) trajectory model + ! of the runge-kutta (2nd order) trajectory model integer(i_kind),parameter::lag_rk2stepnpara_r=29 ! real numbers integer(i_kind),parameter::lag_rk2stepnpara_i=16 ! integer numbers @@ -82,16 +82,16 @@ module lag_traj ! Number of parameters needed by the tangent linear model, for one iteration - ! of the Runge-Kutta (4th order) trajectory model + ! of the runge-kutta (4th order) trajectory model integer(i_kind)::lag_rk4itenpara_r,lag_rk4itenpara_i ! Number of parameters needed by the tangent linear model, for one iteration - ! of the Runge-Kutta (2nd order) trajectory model + ! of the runge-kutta (2nd order) trajectory model integer(i_kind)::lag_rk2itenpara_r,lag_rk2itenpara_i ! Number of steps needed to achieve one iteration integer(i_kind)::lag_nstepiter - ! Array indexes for the storage of the NL RK4 model parameters + ! Array indexes for the storage of the nl rk4 model parameters integer(i_kind),parameter::irk4_loc_0=1 integer(i_kind),parameter::irk4_loc_half_p =9 integer(i_kind),parameter::irk4_loc_half_pp=17 @@ -118,7 +118,7 @@ module lag_traj integer(i_kind),parameter::irk4_lon_half_pp=54 integer(i_kind),parameter::irk4_lon_1_p=56 - ! Array indexes for the storage of the NL RK2 model parameters + ! Array indexes for the storage of the nl rk2 model parameters integer(i_kind),parameter::irk2_loc_0=1 integer(i_kind),parameter::irk2_loc_1_p=9 @@ -298,7 +298,7 @@ function m2lon_nl(rv_orig_lat,rv_dx,lspec_r) real(r_kind) ,intent(in ) :: rv_orig_lat ! distance in meters along a longitude line real(r_kind) ,intent(in ) :: rv_dx - ! parameters needed for the TL and adjoint (optional) + ! parameters needed for the tl and adjoint (optional) real(r_kind),dimension(2),optional,intent( out) :: lspec_r ! equivalent distance in radians real(r_kind)::m2lon_nl @@ -425,7 +425,7 @@ function lag_haversin(rv_theta) real(r_kind),intent(in ) :: rv_theta real(r_kind)::lag_haversin - lag_haversin=sin(rv_theta/2)**2 + lag_haversin=sin(rv_theta/two)**2 end function lag_haversin ! ------------------------------------------------------------------------ ! Inverse-haversine function @@ -586,12 +586,12 @@ end function lag_d_lon ! ------------------------------------------------------------------------ - ! RK4 MODEL IMPLEMENTATION FOR ONE TIME STEP ----------------------------- + ! rk4 model implementation for one time step ----------------------------- ! ------------------------------------------------------------------------ ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm : Non linear + ! Implementation of the runge-kutta algorithm : non linear ! ------------------------------------------------------------------------ subroutine lag_rk4_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) !$$$ subprogram documentation block @@ -627,7 +627,7 @@ subroutine lag_rk4_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) real(r_kind) ,dimension(:,:) ,intent(in ) :: ufield,vfield ! Time step (seconds) real(r_kind) ,intent(in ) :: tstep - ! Parameters for the TL and Adjoint (optional) + ! Parameters for the tl and adjoint (optional) integer(i_kind),dimension(lag_rk4stepnpara_i),optional,intent( out) :: lspec_i real(r_kind) ,dimension(lag_rk4stepnpara_r),optional,intent( out) :: lspec_r @@ -755,7 +755,7 @@ subroutine lag_rk4_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) end if ! Final step : calculate the final position of the balloon using a - ! Simpson's rules corrector + ! simpson's rules corrector if (lv_spec) then rv_intu_1_p=lag_int3d_nl(ufield,& pos_1_p(1),pos_1_p(2),pos_1_p(3),& @@ -813,7 +813,7 @@ end subroutine lag_rk4_nl ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm-> TLM + ! Implementation of the runge-kutta algorithm-> tlm ! ------------------------------------------------------------------------ subroutine lag_rk4_tl(lspec_i,lspec_r,lon,lat,p,ufield,vfield) !$$$ subprogram documentation block @@ -937,7 +937,7 @@ end subroutine lag_rk4_tl ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm-> ADJOINT + ! Implementation of the runge-kutta algorithm-> adjoint ! ------------------------------------------------------------------------ subroutine lag_rk4_ad(lspec_i,lspec_r,& ad_lon,ad_lat,ad_p,ad_ufield,ad_vfield) @@ -1006,7 +1006,7 @@ subroutine lag_rk4_ad(lspec_i,lspec_r,& ad_lon=zero; ad_lat=zero; ad_p=zero; ! Final step : calculate the final position of the balloon using a - ! Simpson's rules corretor + ! simpson's rules corretor ad_p=ad_p+ad_pos_1_pp(3) ad_lon=ad_lon+ad_pos_1_pp(1) @@ -1117,13 +1117,13 @@ subroutine lag_rk4_ad(lspec_i,lspec_r,& end subroutine lag_rk4_ad ! ------------------------------------------------------------------------ - ! RK2 MODEL IMPLEMENTATION FOR ONE TIME STEP ----------------------------- + ! rk2 model implementation for one time step ----------------------------- ! ------------------------------------------------------------------------ ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm (2nd order) : Non linear - ! (Heun's Method) + ! Implementation of the runge-kutta algorithm (2nd order) : non linear + ! (heun's method) ! ------------------------------------------------------------------------ subroutine lag_rk2_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) !$$$ subprogram documentation block @@ -1159,7 +1159,7 @@ subroutine lag_rk2_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) real(r_kind) ,dimension(:,:) ,intent(in ) :: ufield,vfield ! Time step (seconds) real(r_kind) ,intent(in ) :: tstep - ! Parameters for the TL and Adjoint (optional) + ! Parameters for the tl and adjoint (optional) integer(i_kind),dimension(lag_rk2stepnpara_i),optional,intent( out) :: lspec_i real(r_kind) ,dimension(lag_rk2stepnpara_r),optional,intent( out) :: lspec_r @@ -1264,8 +1264,8 @@ end subroutine lag_rk2_nl ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm (2nd order) : Tangent linear - ! (Heun's Method) + ! Implementation of the runge-kutta algorithm (2nd order) : tangent linear + ! (heun's method) ! ------------------------------------------------------------------------ subroutine lag_rk2_tl(lspec_i,lspec_r,lon,lat,p,ufield,vfield) !$$$ subprogram documentation block @@ -1294,7 +1294,7 @@ subroutine lag_rk2_tl(lspec_i,lspec_r,lon,lat,p,ufield,vfield) !$$$ end documentation block implicit none - ! Parameters for the TL and Adjoint + ! Parameters for the tl and adjoint integer(i_kind),dimension(lag_rk2stepnpara_i),intent(in ) :: lspec_i real(r_kind) ,dimension(lag_rk2stepnpara_r),intent(in ) :: lspec_r ! longitude, latitude, pressure of the balloon @@ -1357,8 +1357,8 @@ end subroutine lag_rk2_tl ! ------------------------------------------------------------------------ - ! Implementation of the Runge-Kutta algorithm (2nd order) : Adjoint - ! (Heun's Method) + ! Implementation of the runge-kutta algorithm (2nd order) : adjoint + ! (heun's method) ! ------------------------------------------------------------------------ subroutine lag_rk2_ad(lspec_i,lspec_r,ad_lon,ad_lat,ad_p,& ad_ufield,ad_vfield) @@ -1389,7 +1389,7 @@ subroutine lag_rk2_ad(lspec_i,lspec_r,ad_lon,ad_lat,ad_p,& !$$$ end documentation block implicit none - ! Parameters for the TL and Adjoint + ! Parameters for the tl and adjoint integer(i_kind),dimension(lag_rk2stepnpara_i),intent(in ) :: lspec_i real(r_kind) ,dimension(lag_rk2stepnpara_r),intent(in ) :: lspec_r ! longitude, latitude, pressure of the balloon @@ -1481,12 +1481,12 @@ end subroutine lag_rk2_ad ! -------------------------------------------------------------------------- - ! RK2 MODEL IMPLEMENTATION FOR ONE ITERATION (discomposed in sev. time steps) + ! rk2 model implementation for one iteration (discomposed in sev. time steps) ! -------------------------------------------------------------------------- ! ------------------------------------------------------------------------ - ! Implementation for the RK2 non linear model + ! Implementation for the rk2 non linear model ! ------------------------------------------------------------------------ subroutine lag_rk2iter_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) !$$$ subprogram documentation block @@ -1522,7 +1522,7 @@ subroutine lag_rk2iter_nl(lon,lat,p,ufield,vfield,tstep,lspec_i,lspec_r) real(r_kind) ,dimension(:,:) ,intent(in ) :: ufield,vfield ! Time step (seconds) real(r_kind) ,intent(in ) :: tstep - ! Parameters for the TL and Adjoint (optional) + ! Parameters for the tl and adjoint (optional) integer(i_kind),dimension(:),optional,intent( out) :: lspec_i real(r_kind) ,dimension(:),optional,intent( out) :: lspec_r @@ -1598,7 +1598,7 @@ end subroutine lag_rk2iter_nl ! ------------------------------------------------------------------------ - ! Implementation for the RK2 tangent-linear model + ! Implementation for the rk2 tangent-linear model ! ------------------------------------------------------------------------ subroutine lag_rk2iter_tl(lspec_i,lspec_r,lon,lat,p,ufield,vfield) !$$$ subprogram documentation block @@ -1654,7 +1654,7 @@ end subroutine lag_rk2iter_tl ! ------------------------------------------------------------------------ - ! Implementation for the RK2 adjoint model + ! Implementation for the rk2 adjoint model ! ------------------------------------------------------------------------ subroutine lag_rk2iter_ad(lspec_i,lspec_r,ad_lon,ad_lat,ad_p,ad_ufield,ad_vfield) !$$$ subprogram documentation block diff --git a/src/gsi/lanczos.F90 b/src/gsi/lanczos.F90 index d6d7d5bb92..8ca2db183b 100644 --- a/src/gsi/lanczos.F90 +++ b/src/gsi/lanczos.F90 @@ -5,16 +5,16 @@ module lanczos ! prgmmr: tremolet ! ! abstract: Contains variables and routines for preconditioned -! Lanczos minimizer following Mike Fisher's algorithm. +! lanczos minimizer following mike fisher's algorithm. ! ! program history log: ! 2007-05-16 tremolet ! 2007-07-11 tremolet - increment sensitivity to obs ! 2007-11-23 todling - add timers ! 2009-01-18 todling - minimal changes to interface w/ quad-based evaljgrad -! NOTE: no attempt made to reproduce across pe's yet +! Note: no attempt made to reproduce across pe's yet ! 2009-08-18 lueken - update documentation -! 2010-03-10 treadon - add ESSL interface +! 2010-03-10 treadon - add essl interface ! 2010-03-17 todling - add analysis error estimate (congrad_siga) ! 2010-08-19 lueken - add only to module use ! 2010-09-06 todling - add ltcost parameter to allow writing out true cost @@ -31,23 +31,23 @@ module lanczos ! lanczos_precond - Preconditioner itself (called from congrad, internal) ! ! Variable Definitions: -! LMPCGL : .T. ====> precondition conjugate-gradient minimization -! R_MAX_CNUM_PC : Maximum allowed condition number for the preconditioner -! NPCVECS : number of vectors which make up the preconditioner +! lmpcgl : .t. ====> precondition conjugate-gradient minimization +! r_max_cnum_pc : Maximum allowed condition number for the preconditioner +! npcvecs : number of vectors which make up the preconditioner ! -! YVCGLPC: eigenvectors (from an earlier minimization) +! yvcglpc: eigenvectors (from an earlier minimization) ! that are used to construct the preconditioner. -! RCGLPC : eigenvalues (from an earlier minimization) +! rcglpc : eigenvalues (from an earlier minimization) ! that are used to construct the preconditioner. -! NVCGLPC: the number of eigenpairs used to form the preconditioner. +! nvcglpc: the number of eigenpairs used to form the preconditioner. ! -! YVCGLEV: eigenvectors for the current minimization. -! RCGLEV : eigenvalues for the current minimization. -! NVCGLEV: the number of eigenpairs for the current minimization. -! LTCOST : .T. to calculate true cost function (unscalled), this adds +! yvcglev: eigenvectors for the current minimization. +! rcglev : eigenvalues for the current minimization. +! nvcglev: the number of eigenpairs for the current minimization. +! ltcost : .t. to calculate true cost function (unscalled), this adds ! considerable computation cost; only used in test mode ! -! YVCGLWK: work array of eigenvectors +! yvcglwk: work array of eigenvectors ! ! attributes: ! language: f90 @@ -71,36 +71,36 @@ module lanczos ! ------------------------------------------------------------------------------ -logical :: LTCOST_= .false. -logical :: LMPCGL = .false. -logical :: LCONVERT = .false. !if true, convert the preconditioner vectors for the next outer loop. -logical :: LDECOMP = .false. !if true, carry lanczos decomposition at each iteration -real(r_kind) :: R_MAX_CNUM_PC = 10.0_r_kind +logical :: ltcost_= .false. +logical :: lmpcgl = .false. +logical :: lconvert = .false. !if true, convert the preconditioner vectors for the next outer loop. +logical :: ldecomp = .false. !if true, carry lanczos decomposition at each iteration +real(r_kind) :: r_max_cnum_pc = 10.0_r_kind real(r_kind) :: xmin_ritz = one real(r_kind) :: pkappa = one_tenth -integer(i_kind) :: NPCVECS, NVCGLPC, NVCGLEV, NWRVECS -REAL(r_kind), ALLOCATABLE :: RCGLPC(:) -REAL(r_kind), ALLOCATABLE :: RCGLEV(:) +integer(i_kind) :: npcvecs, nvcglpc, nvcglev, nwrvecs +real(r_kind), allocatable :: rcglpc(:) +real(r_kind), allocatable :: rcglev(:) integer(i_kind) :: mype,nprt,jiter,maxiter logical :: l4dvar,lanczosave -REAL(r_kind), allocatable :: zlancs(:,:) +real(r_kind), allocatable :: zlancs(:,:) -TYPE(control_vector), ALLOCATABLE, DIMENSION(:) :: YVCGLPC -TYPE(control_vector), ALLOCATABLE, DIMENSION(:) :: YVCGLEV -TYPE(control_vector), ALLOCATABLE, DIMENSION(:) :: YVCGLWK +type(control_vector), allocatable, dimension(:) :: yvcglpc +type(control_vector), allocatable, dimension(:) :: yvcglev +type(control_vector), allocatable, dimension(:) :: yvcglwk type(control_vector), allocatable, dimension(:) :: cglwork ! -------------------------------------- -integer(i_kind), PARAMETER :: N_DEFAULT_REAL_KIND = r_single -integer(i_kind), PARAMETER :: N_DOUBLE_KIND = r_double +integer(i_kind), parameter :: n_default_real_kind = r_single +integer(i_kind), parameter :: n_double_kind = r_double ! -------------------------------------- ! ------------------------------------------------------------------------------ contains ! ------------------------------------------------------------------------------ -! CONGRAD +! congrad ! ------------------------------------------------------------------------------ subroutine setup_congrad(kpe,kprt,kiter,kiterstart,kmaxit,kwrvecs, & ld4dvar,ldsave,ltcost) @@ -147,15 +147,15 @@ subroutine setup_congrad(kpe,kprt,kiter,kiterstart,kmaxit,kwrvecs, & zlancs=zero allocate(cglwork(maxiter+1)) -DO ii=1,kmaxit+1 - CALL allocate_cv(cglwork(ii)) +do ii=1,kmaxit+1 + call allocate_cv(cglwork(ii)) cglwork(ii)=zero -ENDDO +enddo if (jiter==kiterstart) then - NPCVECS=0 - NVCGLPC=0 - NVCGLEV=0 + npcvecs=0 + nvcglpc=0 + nvcglev=0 endif if (jiter>1) call setup_precond() @@ -196,7 +196,7 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) ! !$$$ end documentation block -IMPLICIT NONE +implicit none type(control_vector), intent(inout) :: xhat real(r_kind) , intent( out) :: pcost @@ -253,23 +253,23 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) !--- change of variable to account for preconditioning -if (LMPCGL) call lanczos_precond(xhat,+2) +if (lmpcgl) call lanczos_precond(xhat,+2) -zgnorm = SQRT( DOT_PRODUCT (gradx,gradx)) +zgnorm = sqrt( dot_product (gradx,gradx)) if (mype==0) write (6,*)'grepmin Starting point: Estimated gradient norm=',zgnorm -if (LMPCGL) call lanczos_precond(gradx,-2) +if (lmpcgl) call lanczos_precond(gradx,-2) cglwork(1) = gradx -znorm2l1 = DOT_PRODUCT(cglwork(1),cglwork(1)) -cglwork(1)%values = cglwork(1)%values / SQRT(znorm2l1) +znorm2l1 = dot_product(cglwork(1),cglwork(1)) +cglwork(1)%values = cglwork(1)%values / sqrt(znorm2l1) !--- save initial control vector and gradient grad0 = gradx -zqg0(1) = DOT_PRODUCT(cglwork(1),grad0) +zqg0(1) = dot_product(cglwork(1),grad0) if(nprt>=1.and.ltcost_) then if (mype==0) then @@ -282,21 +282,21 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) ingood = 0 iter = 1 -Lanczos_loop : DO +lanczos_loop : do -!--- evaluate the Hessian applied to the latest Lanczos vector +!--- evaluate the hessian applied to the latest lanczos vector do jj=1,zww%lencv zww%values(jj) = xhat%values(jj) + cglwork(iter)%values(jj) enddo - if (LMPCGL) call lanczos_precond(zww,-2) + if (lmpcgl) call lanczos_precond(zww,-2) lsavinc=.false. call evaljgrad(zww,pcostq,gradx,lsavinc,iprt,myname) pcost=pcostq - if (LMPCGL) call lanczos_precond(gradx,-2) + if (lmpcgl) call lanczos_precond(gradx,-2) do jj=1,gradx%lencv gradx%values(jj) = gradx%values(jj) - grad0%values(jj) @@ -304,15 +304,15 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) !--- calculate zdelta - zdelta(iter) = DOT_PRODUCT(cglwork(iter),gradx) + zdelta(iter) = dot_product(cglwork(iter),gradx) if (zdelta(iter)<=zero) then if (mype==0) write(6,*)'congrad stopping: J" not positive definite',zdelta(iter) iter = iter-1 - EXIT Lanczos_loop + exit lanczos_loop endif -!--- Calculate the new Lanczos vector (This is the Lanczos recurrence) +!--- Calculate the new lanczos vector (this is the lanczos recurrence) do jj=1,gradx%lencv gradx%values(jj) = gradx%values(jj) - zdelta(iter) * cglwork(iter)%values(jj) @@ -326,19 +326,19 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) !--- orthonormalize gradient against previous gradients do jm=iter,1,-1 - zdla = DOT_PRODUCT(gradx,cglwork(jm)) + zdla = dot_product(gradx,cglwork(jm)) do jj=1,gradx%lencv gradx%values(jj) = gradx%values(jj) - zdla*cglwork(jm)%values(jj) enddo enddo - zbeta(iter+1) = SQRT(DOT_PRODUCT(gradx,gradx)) + zbeta(iter+1) = sqrt(dot_product(gradx,gradx)) do jj=1,gradx%lencv cglwork(iter+1)%values(jj) = gradx%values(jj) / zbeta(iter+1) enddo - zqg0(iter+1) = DOT_PRODUCT(cglwork(iter+1),grad0) + zqg0(iter+1) = dot_product(cglwork(iter+1),grad0) !--- calculate the reduction in the gradient norm and cost @@ -359,16 +359,16 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) enddo enddo - if (LMPCGL) call lanczos_precond(zww,+2) + if (lmpcgl) call lanczos_precond(zww,+2) - preduc_norm = SQRT(DOT_PRODUCT(zww,zww)) + preduc_norm = sqrt(dot_product(zww,zww)) preduc = preduc_norm/zgnorm if (mype==0) write (6,'(2(1X,A,ES25.18))') & 'Estimated gradient norm=',preduc_norm,' reduction = ',preduc !--- determine eigenvalues and eigenvectors of the tri-diagonal problem - if((LDECOMP .or. (iter==kmaxit)) .and. lsavevecs) then + if((ldecomp .or. (iter==kmaxit)) .and. lsavevecs) then zlancs(1:iter ,4) = zdelta(1:iter) zlancs(1:iter-1,1) = zbeta (2:iter) @@ -389,9 +389,9 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) zbnds(1:iter) = abs(zbeta(iter+1)*zv(iter,1:iter)) if (mype==0) write (6,*)'congrad: error bounds are: ',zbnds(1:iter) -!--- Check for exploding or negative Ritz values +!--- Check for exploding or negative ritz values - if (ANY(zritz(1:iter)0) then @@ -442,26 +442,26 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) if (mype==0) write(6,*)'congrad: End of iteration: ',iter if (mype==0) write(6,'(/)') -! count how many eigenpairs have converged to PKAPPA precision and have +! count how many eigenpairs have converged to pkappa precision and have ! eigenvalue > xmin_ritz (which is 1 by default) ! (For the analysis, all eigenvalues should be >1. For the singular vector calculation, ! we are not interested in decaying modes.) -! However, when SVs are computed with projection operators, 1 may not +! However, when svs are computed with projection operators, 1 may not ! be an appropriate choice for xmin_ritz - imaxevecs = COUNT(zbnds(1:iter)/zritz(1:iter)<=pkappa .AND. zritz(1:iter)>xmin_ritz) + imaxevecs = count(zbnds(1:iter)/zritz(1:iter)<=pkappa .and. zritz(1:iter)>xmin_ritz) if (imaxevecs >= kmaxevecs) then if (mype==0) write(6,*)imaxevecs,' eigenpairs converged to precision ',pkappa if (mype==0) write(6,'(/)') - EXIT Lanczos_loop + exit lanczos_loop endif end if ! Tests for end of iterations if (iter >= kmaxit .or. (preduc <= zreqrd .and. iter >= kminit)) & - EXIT Lanczos_loop + exit lanczos_loop @@ -484,7 +484,7 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) call allocate_cv(gradf) gradf=zero - if (LMPCGL) then + if (lmpcgl) then call lanczos_precond(xiter,-2) endif call evaljgrad(xiter,pcostq,gradf,lsavinc,nprt,myname) @@ -513,7 +513,7 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) xiter%values(ii) = xiter%values(ii) + cglwork(jj)%values(ii)*zlancs(jj,3) enddo enddo - if (LMPCGL) call lanczos_precond(xiter,-2) + if (lmpcgl) call lanczos_precond(xiter,-2) xsens=xiter ! Compute observation impact @@ -531,7 +531,7 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) iter = iter+1 if (ingood>0) itheta1 = itheta1+1 -ENDDO Lanczos_loop +enddo lanczos_loop !--- end of Lanczos iteration @@ -572,7 +572,7 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) !--- transform control variable and gradient back to unpreconditioned space -if (LMPCGL) then +if (lmpcgl) then call lanczos_precond(xhat ,-2) call lanczos_precond(gradx,+2) endif @@ -593,15 +593,15 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) if (lanczosave) then do jj=1,iter clfile='lanczvec.XXX.YYYY' - WRITE(clfile(10:12),'(I3.3)') jiter - WRITE(clfile(14:17),'(I4.4)') jj + write(clfile(10:12),'(I3.3)') jiter + write(clfile(14:17),'(I4.4)') jj call write_cv(cglwork(jj),clfile) - ENDDO + enddo if (mype==0) then iunit=get_lun() clfile='zlanczos.XXX' - WRITE(clfile(10:12),'(I3.3)') jiter + write(clfile(10:12),'(I3.3)') jiter write(6,*)'Writing Lanczos coef. to file ',clfile open(iunit,file=trim(clfile),form='unformatted') @@ -616,54 +616,54 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) if (l4dvar .and. lsavevecs) then zbnds(1:iter) = zbnds(1:iter)/zritz(1:iter) - NVCGLEV = 0 + nvcglev = 0 do jk=iter,1,-1 if (zbnds(jk) <= pkappa .AND. zritz(jk) > xmin_ritz) then - NVCGLEV=NVCGLEV+1 + nvcglev=nvcglev+1 endif - ENDDO + enddo if (mype==0) write(6,*) & - 'Number of eigenpairs converged to requested accuracy NVCGLEV=',NVCGLEV + 'Number of eigenpairs converged to requested accuracy NVCGLEV=',nvcglev - NVCGLEV= min(NVCGLEV,NWRVECS) + nvcglev= min(nvcglev,nwrvecs) if (mype==0) write(6,*) & - 'Number of eigenevctors to be calculated is NVCGLEV=',NVCGLEV + 'Number of eigenevctors to be calculated is NVCGLEV=',nvcglev - ALLOCATE(RCGLEV(NVCGLEV)) - ALLOCATE (YVCGLEV(NVCGLEV)) - DO ii=1,NVCGLEV - CALL allocate_cv(YVCGLEV(ii)) - ENDDO + allocate(rcglev(nvcglev)) + allocate(yvcglev(nvcglev)) + do ii=1,nvcglev + call allocate_cv(yvcglev(ii)) + enddo ii=0 do jk=iter,1,-1 - if (zbnds(jk) <= pkappa .AND. zritz(jk) > xmin_ritz .AND. ii xmin_ritz .and. ii0) then + if (mype==0.and.nvcglev>0) then write(6,'(/)') write(6,*)'Calculated eigenvectors for the following eigenvalues:' - write(6,*)'RCGLEV=',RCGLEV(1:NVCGLEV) + write(6,*)'RCGLEV=',rcglev(1:nvcglev) write(6,'(/)') endif @@ -686,8 +686,8 @@ subroutine congrad(xhat,pcost,gradx,preduc,kmaxit,iobsconv,lsavevecs) !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- -! STEQR - Simplified interface to LAPACK routines SSTEQR/DSTEQR -! and to ESSL routines SGEEV/DGEEV +! steqr - Simplified interface to lapack routines ssteqr/dsteqr +! and to essl routines sgeev/dgeev !----------------------------------------------------------------------- subroutine steqr !$$$ subprogram documentation block @@ -699,7 +699,7 @@ subroutine steqr ! ! program history log: ! 2009-08-05 lueken - added subprogram doc block -! 2010-03-10 treadon - add ESSL interface +! 2010-03-10 treadon - add essl interface ! ! input argument list: ! @@ -722,7 +722,7 @@ subroutine steqr real(r_kind),allocatable:: a(:,:),aux(:),w_order(:) complex(r_kind),allocatable:: w(:),z(:,:) -! Use ESSL +! Use essl iopt=1 n=iter lda=n @@ -730,7 +730,7 @@ subroutine steqr naux=2*n allocate(select(n),indx(n),a(lda,n),w(n),z(ldz,n),aux(naux),w_order(n)) -! Load matrix A. +! Load matrix a. a=zero do i=1,n-1 a(i, i)=zlancs(i,4) ! load diagonal @@ -745,25 +745,25 @@ subroutine steqr z=zero aux=zero -! Call ESSL routines - if (r_kind == N_DEFAULT_REAL_KIND) then - call SGEEV(iopt, a, lda, w, z, ldz, select, n, aux, naux) +! Call essl routines + if (r_kind == n_default_real_kind) then + call sgeev(iopt, a, lda, w, z, ldz, select, n, aux, naux) do i=1,n w_order(i)=real(w(i),r_kind) end do - call SSORTX(w_order,1,n,indx) ! sort eigenvalues into ascending order - ELSEIF (r_kind == N_DOUBLE_KIND) then - call DGEEV(iopt, a, lda, w, z, ldz, select, n, aux, naux) + call ssortx(w_order,1,n,indx) ! sort eigenvalues into ascending order + elseif (r_kind == n_double_kind) then + call dgeev(iopt, a, lda, w, z, ldz, select, n, aux, naux) do i=1,n w_order(i)=real(w(i),r_kind) end do - call DSORTX(w_order,1,n,indx) ! sort eigenvalues into ascending order + call dsortx(w_order,1,n,indx) ! sort eigenvalues into ascending order else write(6,*)'STEQR: r_kind is neither default real nor double precision' call stop2(319) endif -! Load ESSL eigenvalues and eigenvectors into output arrays +! Load essl eigenvalues and eigenvectors into output arrays do j=1,n zlancs(j,4)=w_order(j) ! eigenvalues jj=indx(j) @@ -777,11 +777,11 @@ subroutine steqr #else -! Use LAPACK - if (r_kind == N_DEFAULT_REAL_KIND) then - call SSTEQR ('I',iter,zlancs(1,4),zlancs,zv,kmaxit+1,zsstwrk,info) - ELSEIF (r_kind == N_DOUBLE_KIND) then - call DSTEQR ('I',iter,zlancs(1,4),zlancs,zv,kmaxit+1,zsstwrk,info) +! Use lapack + if (r_kind == n_default_real_kind) then + call ssteqr ('I',iter,zlancs(1,4),zlancs,zv,kmaxit+1,zsstwrk,info) + elseif (r_kind == n_double_kind) then + call dsteqr ('I',iter,zlancs(1,4),zlancs,zv,kmaxit+1,zsstwrk,info) else write(6,*)'STEQR: r_kind is neither default real nor double precision' call stop2(319) @@ -797,8 +797,8 @@ subroutine steqr end subroutine steqr !----------------------------------------------------------------------- -! PTSV - Simplified interface to LAPACK routines SPTSV/DPTSV -! and to ESSL routines SPTF,S/DPTF,S +! ptsv - Simplified interface to lapack routines sptsv/dptsv +! and to essl routines sptf,s/dptf,s !----------------------------------------------------------------------- subroutine ptsv !$$$ subprogram documentation block @@ -810,7 +810,7 @@ subroutine ptsv ! ! program history log: ! 2009-08-05 lueken - added subprogram doc block -! 2010-03-10 treadon- add ESSL interface +! 2010-03-10 treadon- add essl interface ! ! input argument list: ! @@ -829,7 +829,7 @@ subroutine ptsv integer(i_kind) :: i,n,iopt real(r_kind),allocatable:: c(:),d(:),bx(:) -! Use ESSL +! Use essl iopt=0 n=iter allocate(c(n),d(n),bx(n)) @@ -837,18 +837,18 @@ subroutine ptsv ! Load matrices c=zero do i=1,n-1 - c(i+1) = zlancs(i+1,2) ! lower subdiagonal of A - d(i) = zlancs(i,1) ! main diagonal of A - bx(i)=zlancs(i,3) ! right hand side B + c(i+1) = zlancs(i+1,2) ! lower subdiagonal of a + d(i) = zlancs(i,1) ! main diagonal of a + bx(i)=zlancs(i,3) ! right hand side b end do d(n) =zlancs(n,1) bx(n)=zlancs(n,3) -! Factorize and solve system of equations using ESSL routines - if (r_kind == N_DEFAULT_REAL_KIND) then +! Factorize and solve system of equations using essl routines + if (r_kind == n_default_real_kind) then call sptf(n,c,d,iopt) ! factorize call spts(n,c,d,bx) ! solve - ELSEIF (r_kind == N_DOUBLE_KIND) then + elseif (r_kind == n_double_kind) then call dptf(n,c,d,iopt) ! factorize call dpts(n,c,d,bx) ! solve else @@ -856,7 +856,7 @@ subroutine ptsv call stop2(321) endif -! Load ESSL result into output arrays +! Load essl result into output arrays do i=1,n zlancs(i,3)=bx(i) end do @@ -866,11 +866,11 @@ subroutine ptsv #else -! Use LAPACK - if (r_kind == N_DEFAULT_REAL_KIND) then - call SPTSV (iter,1,zlancs(1,1),zlancs(2,2),zlancs(1,3),kmaxit+1,info) - ELSEIF (r_kind == N_DOUBLE_KIND) then - call DPTSV (iter,1,zlancs(1,1),zlancs(2,2),zlancs(1,3),kmaxit+1,info) +! Use lapack + if (r_kind == n_default_real_kind) then + call sptsv (iter,1,zlancs(1,1),zlancs(2,2),zlancs(1,3),kmaxit+1,info) + elseif (r_kind == n_double_kind) then + call dptsv (iter,1,zlancs(1,1),zlancs(2,2),zlancs(1,3),kmaxit+1,info) else write(6,*) 'r_kind is neither default real nor double precision' call stop2(321) @@ -896,7 +896,7 @@ subroutine congrad_ad(xsens,kiter) ! subprogram: congrad_ad ! prgmmr: ! -! abstract: Apply product of adjoint of estimated Hessian to a vector. +! abstract: Apply product of adjoint of estimated hessian to a vector. ! ! program history log: ! 2009-08-05 lueken - added subprogram doc block @@ -928,7 +928,7 @@ subroutine congrad_ad(xsens,kiter) zzz=dot_product(xsens,xsens) if (mype==0) write(6,888)'congrad_ad: Norm input=',sqrt(zzz) -if (LMPCGL) call lanczos_precond(xsens,-2) +if (lmpcgl) call lanczos_precond(xsens,-2) zaa=zero do jj=1,kiter @@ -948,7 +948,7 @@ subroutine congrad_ad(xsens,kiter) enddo enddo -if (LMPCGL) call lanczos_precond(xsens,-2) +if (lmpcgl) call lanczos_precond(xsens,-2) zzz=dot_product(xsens,xsens) if (mype==0) write(6,888)'congrad_ad: Norm output=',sqrt(zzz) @@ -973,7 +973,7 @@ subroutine congrad_siga(siga,ivecs,rc) ! program history log: ! 2010-03-17 todling - initia code ! 2010-05-16 todling - update to use gsi_bundle -! 2013-01-26 parrish - WCOSS debug compile flagged type mismatch error for +! 2013-01-26 parrish - wcoss debug compile flagged type mismatch error for ! "call bkg_stddev(aux,mval(ii))". ! I changed to ! "call bkg_stddev(aux%step(ii),mval(ii))". @@ -1009,12 +1009,12 @@ subroutine congrad_siga(siga,ivecs,rc) real(r_kind) :: zz rc=0 -NPCVECS = NVCGLEV -ivecs=MIN(npcvecs,nwrvecs) +npcvecs = nvcglev +ivecs=min(npcvecs,nwrvecs) if (ivecs<1) then - if (mype==0) write(6,*)'save_precond: cannot get siga, ivecs=', ivecs - rc=1 - return + if (mype==0) write(6,*)'save_precond: cannot get siga, ivecs=', ivecs + rc=1 + return endif call allocate_preds(sbias) @@ -1023,22 +1023,22 @@ subroutine congrad_siga(siga,ivecs,rc) end do call allocate_cv(aux) -!-- calculate increment on analysis error covariance diag(Delta P) +!-- calculate increment on analysis error covariance diag(delta p) siga=zero -DO jj=1,ivecs - zz=sqrt(one-one/sqrt(RCGLEV(jj))) - aux%values = zz * YVCGLEV(jj)%values - call control2model(aux,mval,sbias) - do ii=1,nsubwin - call gsi_bundlehadamard(siga,mval(ii),mval(ii)) - enddo -ENDDO +do jj=1,ivecs + zz=sqrt(one-one/sqrt(rcglev(jj))) + aux%values = zz * yvcglev(jj)%values + call control2model(aux,mval,sbias) + do ii=1,nsubwin + call gsi_bundlehadamard(siga,mval(ii),mval(ii)) + enddo +enddo do ii=1,nsubwin -!-- get B standard deviations +!-- get b standard deviations call bkg_stddev(aux%step(ii),mval(ii)) -!-- calculate diag(Pa) = diag(B) - diag(Delta P) -! i.e., add diag(B) as rank-1 update to diag(Delta P) +!-- calculate diag(pa) = diag(b) - diag(delta p) +! i.e., add diag(b) as rank-1 update to diag(delta p) call gsi_bundlehadamard(siga,mval(ii),mval(ii)) enddo @@ -1052,7 +1052,7 @@ end subroutine congrad_siga ! ------------------------------------------------------------------------------ ! ------------------------------------------------------------------------------ -! SAVE_PRECOND - Save eigenvectors from CONGRAD for next minimization +! save_precond - Save eigenvectors from congrad for next minimization ! ------------------------------------------------------------------------------ subroutine save_precond(ldsave) !$$$ subprogram documentation block @@ -1076,90 +1076,90 @@ subroutine save_precond(ldsave) ! !$$$ end documentation block -IMPLICIT NONE +implicit none logical, intent(in ) :: ldsave -REAL(r_kind), ALLOCATABLE :: zmat(:,:) -INTEGER(i_kind) :: ii,jj, info, iunit, ivecs -REAL(r_kind) :: zz -CHARACTER(LEN=13) :: clfile +real(r_kind), allocatable :: zmat(:,:) +integer(i_kind) :: ii,jj, info, iunit, ivecs +real(r_kind) :: zz +character(len=13) :: clfile if (ldsave) then !--- read eigenvalues of the preconditioner - NPCVECS = NVCGLEV+NVCGLPC + npcvecs = nvcglev+nvcglpc if (mype==0) write(6,*)'save_precond: NVCGLEV,NVCGLPC,NPCVECS=', & - & NVCGLEV,NVCGLPC,NPCVECS + nvcglev,nvcglpc,npcvecs - ALLOCATE(YVCGLWK(npcvecs)) + allocate(yvcglwk(npcvecs)) ii=0 - if(.not.lCONVERT) then - DO jj=1,NVCGLEV + if(.not.lconvert) then + do jj=1,nvcglev ii=ii+1 - ! zz=sqrt(RCGLPC(jj)-one) - CALL allocate_cv(YVCGLWK(ii)) - YVCGLWK(ii)%values = YVCGLEV(jj)%values - CALL deallocate_cv(YVCGLEV(jj)) - ENDDO - IF (ALLOCATED(YVCGLEV)) DEALLOCATE(YVCGLEV) + ! zz=sqrt(rcglpc(jj)-one) + call allocate_cv(yvcglwk(ii)) + yvcglwk(ii)%values = yvcglev(jj)%values + call deallocate_cv(yvcglev(jj)) + enddo + if (allocated(yvcglev)) deallocate(yvcglev) - NVCGLEV=0 + nvcglev=0 else !--- copy preconditioner vectors to work file - if (mype==0.and.NVCGLPC>0) write(6,*)'save_precond: RCGLPC=',RCGLPC - DO jj=1,NVCGLPC + if (mype==0.and.nvcglpc>0) write(6,*)'save_precond: RCGLPC=',rcglpc + do jj=1,nvcglpc ii=ii+1 - zz=sqrt(RCGLPC(jj)-one) - CALL allocate_cv(YVCGLWK(ii)) - YVCGLWK(ii)%values = zz * YVCGLPC(jj)%values - CALL deallocate_cv(YVCGLPC(jj)) - ENDDO - IF (ALLOCATED(YVCGLPC)) DEALLOCATE(YVCGLPC) -! IF (ALLOCATED( RCGLPC)) deallocate( RCGLPC) - NVCGLPC=0 - -!--- copy and transform eigenvectors of preconditioned Hessian - - if (mype==0.and.NVCGLEV>0) write(6,*)'save_precond: RCGLEV=',RCGLEV - DO jj=1,NVCGLEV + zz=sqrt(rcglpc(jj)-one) + call allocate_cv(yvcglwk(ii)) + yvcglwk(ii)%values = zz * yvcglpc(jj)%values + call deallocate_cv(yvcglpc(jj)) + enddo + if (allocated(yvcglpc)) deallocate(yvcglpc) +! if (allocated( rcglpc)) deallocate( rcglpc) + nvcglpc=0 + +!--- copy and transform eigenvectors of preconditioned hessian + + if (mype==0.and.nvcglev>0) write(6,*)'save_precond: RCGLEV=',rcglev + do jj=1,nvcglev ii=ii+1 - zz=sqrt(RCGLEV(jj)-one) - CALL allocate_cv(YVCGLWK(ii)) - YVCGLWK(ii)%values = zz * YVCGLEV(jj)%values - CALL deallocate_cv(YVCGLEV(jj)) - ENDDO - IF (ALLOCATED(YVCGLEV)) DEALLOCATE(YVCGLEV) - IF (ALLOCATED( RCGLEV)) deallocate( RCGLEV) - NVCGLEV=0 + zz=sqrt(rcglev(jj)-one) + call allocate_cv(yvcglwk(ii)) + yvcglwk(ii)%values = zz * yvcglev(jj)%values + call deallocate_cv(yvcglev(jj)) + enddo + if (allocated(yvcglev)) deallocate(yvcglev) + if (allocated( rcglev)) deallocate( rcglev) + nvcglev=0 if (mype==0) write(6,*)'save_precond: NVCGLPC,NVCGLEV,npcvecs,ii=', & - NVCGLPC,NVCGLEV,npcvecs,ii + nvcglpc,nvcglev,npcvecs,ii if (ii/=npcvecs) then write(6,*)'save_precond: error number of vectors',ii,npcvecs call stop2(139) end if -!--- form the inner matrix for the Shermann-Morrison-Woodbury inversion +!--- form the inner matrix for the shermann-morrison-woodbury inversion - ALLOCATE(zmat(npcvecs,npcvecs)) + allocate(zmat(npcvecs,npcvecs)) do jj=1,npcvecs do ii=jj,npcvecs - zmat(ii,jj) = DOT_PRODUCT (YVCGLWK(jj),YVCGLWK(ii)) - ENDDO + zmat(ii,jj) = dot_product (yvcglwk(jj),yvcglwk(ii)) + enddo zmat(jj,jj) = zmat(jj,jj) + one - ENDDO + enddo !--- Cholesky decompose if (mype==0) write(6,*)'save_precond: call dpotrf npcvecs=',npcvecs - if (r_kind==N_DEFAULT_REAL_KIND) then - call SPOTRF('L',npcvecs,zmat,npcvecs,info) - ELSEIF (r_kind==N_DOUBLE_KIND) then - call DPOTRF('L',npcvecs,zmat,npcvecs,info) + if (r_kind==n_default_real_kind) then + call spotrf('L',npcvecs,zmat,npcvecs,info) + elseif (r_kind==n_double_kind) then + call dpotrf('L',npcvecs,zmat,npcvecs,info) else write(6,*)'save_precond: r_kind is neither default real nor double precision' call stop2(323) @@ -1175,54 +1175,54 @@ subroutine save_precond(ldsave) do jj=1,npcvecs do ii=1,jj-1 - YVCGLWK(jj)%values = YVCGLWK(jj)%values - zmat(jj,ii)*YVCGLWK(ii)%values + yvcglwk(jj)%values = yvcglwk(jj)%values - zmat(jj,ii)*yvcglwk(ii)%values enddo - YVCGLWK(jj)%values = YVCGLWK(jj)%values / zmat(jj,jj) - ENDDO + yvcglwk(jj)%values = yvcglwk(jj)%values / zmat(jj,jj) + enddo endif !--- Save the eigenvectors if (l4dvar) then - ivecs=MIN(npcvecs,nwrvecs) - DO jj=1,ivecs + ivecs=min(npcvecs,nwrvecs) + do jj=1,ivecs clfile='evec.XXX.YYYY' - WRITE(clfile(6:8) ,'(I3.3)') jiter - WRITE(clfile(10:13),'(I4.4)') jj - call write_cv(YVCGLWK(jj),clfile) - ENDDO + write(clfile(6:8) ,'(I3.3)') jiter + write(clfile(10:13),'(I4.4)') jj + call write_cv(yvcglwk(jj),clfile) + enddo if (mype==0) then iunit=78 clfile='eval.XXX' - WRITE(clfile(6:8),'(I3.3)') jiter + write(clfile(6:8),'(I3.3)') jiter open(iunit,file=clfile) - write(iunit,*)RCGLEV + write(iunit,*)rcglev close(iunit) endif if (mype==0) then iunit=78 clfile='numpcvecs.XXX' - WRITE(clfile(11:13),'(I3.3)') jiter + write(clfile(11:13),'(I3.3)') jiter open(iunit,file=clfile) write(iunit,*)ivecs close(iunit) endif - DO ii=1,npcvecs - CALL deallocate_cv(YVCGLWK(ii)) - ENDDO - DEALLOCATE(YVCGLWK) + do ii=1,npcvecs + call deallocate_cv(yvcglwk(ii)) + enddo + deallocate(yvcglwk) else do ii=nwrvecs+1,npcvecs - CALL deallocate_cv(YVCGLWK(ii)) + call deallocate_cv(yvcglwk(ii)) enddo - npcvecs=MIN(npcvecs,nwrvecs) + npcvecs=min(npcvecs,nwrvecs) endif - if (ALLOCATED(zmat)) DEALLOCATE(zmat) + if (allocated(zmat)) deallocate(zmat) endif do ii=1,maxiter+1 @@ -1234,7 +1234,7 @@ subroutine save_precond(ldsave) end subroutine save_precond ! ------------------------------------------------------------------------------ -! SETUP_PRECOND - Calculates the preconditioner for congrad +! setup_precond - Calculates the preconditioner for congrad ! ------------------------------------------------------------------------------ subroutine setup_precond() !$$$ subprogram documentation block @@ -1246,7 +1246,7 @@ subroutine setup_precond() ! ! program history log: ! 2009-08-05 lueken - added subprogram doc block -! 2010-03-10 treadon - add ESSL interface +! 2010-03-10 treadon - add essl interface ! ! input argument list: ! @@ -1258,315 +1258,315 @@ subroutine setup_precond() ! !$$$ end documentation block - IMPLICIT NONE +implicit none - INTEGER(i_kind), allocatable :: indarr(:) - REAL(r_kind), allocatable :: zq(:),zlam(:),zU(:,:),zUUT(:,:),zwork(:),zzz(:) - INTEGER(i_kind) :: info,ik,inpcv,ji,jj,jk,ii,iunit - REAL(r_kind) :: za, zps - CHARACTER(LEN=13) :: clfile +integer(i_kind), allocatable :: indarr(:) +real(r_kind), allocatable :: zq(:),zlam(:),zu(:,:),zuut(:,:),zwork(:),zzz(:) +integer(i_kind) :: info,ik,inpcv,ji,jj,jk,ii,iunit +real(r_kind) :: za, zps +character(len=13) :: clfile #ifdef ibm_sp -! Declare variables and Work arrays used by ESSL - integer(i_kind):: iopt, ldz, n, naux - real(r_kind), allocatable :: w(:),z(:,:),ap(:),aux(:) +! Declare variables and work arrays used by essl +integer(i_kind):: iopt, ldz, n, naux +real(r_kind), allocatable :: w(:),z(:,:),ap(:),aux(:) #endif !--- read vectors, apply change of variable and copy to work file - if (l4dvar) then - iunit=78 - clfile='numpcvecs.XXX' - WRITE(clfile(11:13),'(I3.3)') jiter-1 - open(iunit,file=clfile) - read(iunit,*)npcvecs - close(iunit) - - if (npcvecs<1) then - write(6,*)'SETUP_PRECOND: no vectors for preconditioner',npcvecs - call stop2(140) - end if +if (l4dvar) then + iunit=78 + clfile='numpcvecs.XXX' + write(clfile(11:13),'(I3.3)') jiter-1 + open(iunit,file=clfile) + read(iunit,*)npcvecs + close(iunit) + + if (npcvecs<1) then + write(6,*)'SETUP_PRECOND: no vectors for preconditioner',npcvecs + call stop2(140) + end if - ALLOCATE(YVCGLWK(npcvecs)) - DO ii=1,npcvecs - CALL allocate_cv(YVCGLWK(ii)) - ENDDO + allocate(yvcglwk(npcvecs)) + do ii=1,npcvecs + call allocate_cv(yvcglwk(ii)) + enddo - do jj=1,npcvecs - clfile='evec.XXX.YYYY' - WRITE(clfile(6:8) ,'(I3.3)') jiter-1 - WRITE(clfile(10:13),'(I4.4)') jj - call read_cv(yvcglwk(jj),clfile) - ENDDO - endif - if(.not. LCONVERT) then - NVCGLPC=npcvecs - if(NVCGLPC > 0) then - ALLOCATE (RCGLPC(NVCGLPC)) + do jj=1,npcvecs + clfile='evec.XXX.YYYY' + write(clfile(6:8) ,'(I3.3)') jiter-1 + write(clfile(10:13),'(I4.4)') jj + call read_cv(yvcglwk(jj),clfile) + enddo +endif +if(.not. lconvert) then + nvcglpc=npcvecs + if(nvcglpc > 0) then + allocate (rcglpc(nvcglpc)) - iunit=78 - clfile='eval.XXX' - WRITE(clfile(6:8),'(I3.3)') jiter-1 - open(iunit,file=clfile) - read(iunit,*)RCGLPC - close(iunit) - do ii=1,NVCGLPC - RCGLPC(ii) = MIN(R_MAX_CNUM_PC,RCGLPC(ii)) - end do - + iunit=78 + clfile='eval.XXX' + write(clfile(6:8),'(I3.3)') jiter-1 + open(iunit,file=clfile) + read(iunit,*)rcglpc + close(iunit) + do ii=1,nvcglpc + rcglpc(ii) = min(r_max_cnum_pc,rcglpc(ii)) + end do + - ALLOCATE (YVCGLPC(NVCGLPC)) - DO jj=1,NVCGLPC - CALL ALLOCATE_CV(YVCGLPC(jj)) - ENDDO - DO jj=1,NVCGLPC - YVCGLPC(jj) = zero - do jk=1,YVCGLPC(jj)%lencv - YVCGLPC(jj)%values(jk) = yvcglwk(jj)%values(jk) - enddo - ENDDO + allocate (yvcglpc(nvcglpc)) + do jj=1,nvcglpc + call allocate_cv(yvcglpc(jj)) + enddo + do jj=1,nvcglpc + yvcglpc(jj) = zero + do jk=1,yvcglpc(jj)%lencv + yvcglpc(jj)%values(jk) = yvcglwk(jj)%values(jk) + enddo + enddo - LMPCGL = .true. - else - NVCGLPC = 0 - LMPCGL = .false. - endif + lmpcgl = .true. else - if (mype==0) write(6,*)'allocate arrays with npcvecs=',npcvecs - allocate(indarr(npcvecs)) - allocate(zq(npcvecs),zlam(npcvecs),zU(npcvecs,npcvecs)) - allocate(zUUT(npcvecs,npcvecs),zwork(3*npcvecs),zzz(npcvecs)) + nvcglpc = 0 + lmpcgl = .false. + endif +else + if (mype==0) write(6,*)'allocate arrays with npcvecs=',npcvecs + allocate(indarr(npcvecs)) + allocate(zq(npcvecs),zlam(npcvecs),zu(npcvecs,npcvecs)) + allocate(zuut(npcvecs,npcvecs),zwork(3*npcvecs),zzz(npcvecs)) -!--- Perform Householder transformations to reduce the matrix of vectors +!--- Perform householder transformations to reduce the matrix of vectors !--- to upper triangular - do jj=1,npcvecs - CALL ALLGATHER_CVSECTION(yvcglwk(jj),zq(1:jj),1,jj) + do jj=1,npcvecs + call allgather_cvsection(yvcglwk(jj),zq(1:jj),1,jj) - zps = DOT_PRODUCT(yvcglwk(jj),yvcglwk(jj)) - DOT_PRODUCT(zq(1:jj),zq(1:jj)) + zps = dot_product(yvcglwk(jj),yvcglwk(jj)) - dot_product(zq(1:jj),zq(1:jj)) - if (zq(jj) < zero) then - zU(jj,jj) = -sqrt(zps+zq(jj)*zq(jj)) - else - zU(jj,jj) = sqrt(zps+zq(jj)*zq(jj)) - endif + if (zq(jj) < zero) then + zu(jj,jj) = -sqrt(zps+zq(jj)*zq(jj)) + else + zu(jj,jj) = sqrt(zps+zq(jj)*zq(jj)) + endif - zq(jj) = zq(jj) - zU(jj,jj) + zq(jj) = zq(jj) - zu(jj,jj) - do jk=1,jj-1 - zU(jk,jj) = zq(jk) - ENDDO + do jk=1,jj-1 + zu(jk,jj) = zq(jk) + enddo - zps = zps + zq(jj)*zq(jj) + zps = zps + zq(jj)*zq(jj) - zzz(1:jj-1)=zero - zzz(jj)=zq(jj) - CALL SET_CVSECTION(zzz(1:jj),yvcglwk(jj),1,jj) + zzz(1:jj-1)=zero + zzz(jj)=zq(jj) + call set_cvsection(zzz(1:jj),yvcglwk(jj),1,jj) - do jk=1,yvcglwk(jj)%lencv - yvcglwk(jj)%values(jk) = yvcglwk(jj)%values(jk) * sqrt(two/zps) - enddo + do jk=1,yvcglwk(jj)%lencv + yvcglwk(jj)%values(jk) = yvcglwk(jj)%values(jk) * sqrt(two/zps) + enddo -!--- we now have the Householder vector in yvcglwk(jj), and the non-zero -!--- elements of the transformed vector in ZU. Now apply the Householder +!--- we now have the householder vector in yvcglwk(jj), and the non-zero +!--- elements of the transformed vector in zu. Now apply the householder !--- transformations to the remaining vectors. - do ji=jj+1,npcvecs - zps = DOT_PRODUCT (yvcglwk(jj),yvcglwk(ji)) - do jk=1,yvcglwk(ji)%lencv - yvcglwk(ji)%values(jk) = yvcglwk(ji)%values(jk) - zps*yvcglwk(jj)%values(jk) - enddo - ENDDO - ENDDO + do ji=jj+1,npcvecs + zps = dot_product (yvcglwk(jj),yvcglwk(ji)) + do jk=1,yvcglwk(ji)%lencv + yvcglwk(ji)%values(jk) = yvcglwk(ji)%values(jk) - zps*yvcglwk(jj)%values(jk) + enddo + enddo + enddo !--- Multiply the upper triangle by its transpose and find eigenvectors !--- and eigenvalues - do jj=1,npcvecs - do ji=jj+1,npcvecs - zU(ji,jj) = zero - enddo + do jj=1,npcvecs + do ji=jj+1,npcvecs + zu(ji,jj) = zero enddo + enddo - do jj=1,npcvecs - do ji=jj,npcvecs - zUUT(ji,jj) = zero - do jk=ji,npcvecs - zUUT(ji,jj) = zUUT(ji,jj) + zU(ji,jk)*zU(jj,jk) - ENDDO - ENDDO - ENDDO + do jj=1,npcvecs + do ji=jj,npcvecs + zuut(ji,jj) = zero + do jk=ji,npcvecs + zuut(ji,jj) = zuut(ji,jj) + zu(ji,jk)*zu(jj,jk) + enddo + enddo + enddo #ifdef ibm_sp -! USE ESSL - iopt=1 - n=npcvecs - ldz=n - naux=3*n - - allocate(w(n),z(n,n),aux(naux),ap(n*n)) - w=zero - z=zero - aux=zero - -! Load zuut in ESSL lower-packed storage mode - ap=zero - jk=0 - do jj=1,n - do ii=jj,n - jk=jk+1 - ap(jk)=zuut(ii,jj) - end do - end do - -! Call ESSL routines - if (r_kind==N_DEFAULT_REAL_KIND) then - call sspev(iopt,ap,w,z,ldz,n,aux,naux) - ELSEIF (r_kind==N_DOUBLE_KIND) then - call dspev(iopt,ap,w,z,ldz,n,aux,naux) - else - write(6,*)'SETUP_PRECOND: r_kind is neither default real nor double precision' - call stop2(325) - endif - -! Load ESSL results into output arrays - do jj=1,n - zlam(jj)=w(jj) - do ii=1,n - zuut(ii,jj)=z(ii,jj) - end do - end do +! Use essl + iopt=1 + n=npcvecs + ldz=n + naux=3*n + + allocate(w(n),z(n,n),aux(naux),ap(n*n)) + w=zero + z=zero + aux=zero + +! Load zuut in essl lower-packed storage mode + ap=zero + jk=0 + do jj=1,n + do ii=jj,n + jk=jk+1 + ap(jk)=zuut(ii,jj) + end do + end do + +! Call essl routines + if (r_kind==n_default_real_kind) then + call sspev(iopt,ap,w,z,ldz,n,aux,naux) + elseif (r_kind==n_double_kind) then + call dspev(iopt,ap,w,z,ldz,n,aux,naux) + else + write(6,*)'SETUP_PRECOND: r_kind is neither default real nor double precision' + call stop2(325) + endif + +! Load essl results into output arrays + do jj=1,n + zlam(jj)=w(jj) + do ii=1,n + zuut(ii,jj)=z(ii,jj) + end do + end do -! Deallocate work arrays - deallocate(w,z,aux,ap) +! Deallocate work arrays + deallocate(w,z,aux,ap) #else -! Use LAPACK routines - if (r_kind==N_DEFAULT_REAL_KIND) then - call SSYEV('V','L',npcvecs,zUUT,npcvecs,zlam,zwork,SIZE(zwork),info) - ELSEIF (r_kind==N_DOUBLE_KIND) then - call DSYEV('V','L',npcvecs,zUUT,npcvecs,zlam,zwork,SIZE(zwork),info) - else - write(6,*)'SETUP_PRECOND: r_kind is neither default real nor double precision' - call stop2(325) - endif - if (info/=0) then - write(6,*)'SETUP_PRECOND: SSYEV/DSYEV returned with info=',info - write(6,*)'SETUP_PRECOND: SSYEV/DSYEV returned non-zero return code' - call stop2(326) - endif - +! Use lapack routines + if (r_kind==n_default_real_kind) then + call ssyev('V','L',npcvecs,zuut,npcvecs,zlam,zwork,size(zwork),info) + elseif (r_kind==n_double_kind) then + call dsyev('V','L',npcvecs,zuut,npcvecs,zlam,zwork,size(zwork),info) + else + write(6,*)'SETUP_PRECOND: r_kind is neither default real nor double precision' + call stop2(325) + endif + if (info/=0) then + write(6,*)'SETUP_PRECOND: SSYEV/DSYEV returned with info=',info + write(6,*)'SETUP_PRECOND: SSYEV/DSYEV returned non-zero return code' + call stop2(326) + endif + #endif !--- convert to eigenvalues of the preconditioner - do jk=1,npcvecs - zlam(jk) = one / (one - zlam(jk)) - ENDDO - - if (mype==0) write(6,*)'SETUP_PRECOND: eigenvalues found are: ',(zlam(ji),ji=1,npcvecs) + do jk=1,npcvecs + zlam(jk) = one / (one - zlam(jk)) + enddo + + if (mype==0) write(6,*)'SETUP_PRECOND: eigenvalues found are: ',(zlam(ji),ji=1,npcvecs) !--- sort eigenvalues with eigenvalues larger than 1 after eigenvalues !--- smaller than 1 and with eigenvalues larger than 1 sorted in decreasing !--- order - do ji=1,npcvecs - indarr(ji) = ji - ENDDO - -!--- straight insertion sort courtesy of Numerical Recipies - - do jj=2,npcvecs - za = zlam(jj) - ik = indarr(jj) - do ji=jj-1,1,-1 - if (zlam(ji)>one .and. (zlam(ji)>=za .or. za<=one)) then - ii=ji - exit - else - ii=0 - endif - zlam(ji+1) = zlam(ji) - indarr(ji+1) = indarr(ji) - ENDDO - zlam(ii+1) = za - indarr(ii+1) = ik - ENDDO - - inpcv = npcvecs - - do while (zlam(inpcv) <= zero) - if (mype==0) write(6,*)'Warning - eigenvalue less than 1: ',zlam(inpcv) - inpcv = inpcv-1 - if (inpcv == 0) then - if (mype==0) write(6,*)'SETUP_PRECOND: cannot form preconditioner - '//& - 'no positive eigenvalues.' - if (mype==0) write(6,*)'SETUP_PRECOND: minimisation will not be preconditioned.' - EXIT - endif - enddo - - IF (inpcv>0) THEN - if (mype==0) write(6,*)'Number of preconditioning vectors selected is ',inpcv - if (mype==0) write(6,*)'SETUP_PRECOND: selected eigenvalues are: ',(zlam(ji),ji=1,inpcv) - - IF (ALLOCATED(YVCGLPC)) THEN - DO jj=1,NVCGLPC - CALL DEALLOCATE_CV(YVCGLPC(jj)) - ENDDO - DEALLOCATE(YVCGLPC) - NVCGLPC=0 - ENDIF - IF (ALLOCATED(RCGLPC)) DEALLOCATE(RCGLPC) - - !--- Save eigenvalues - NVCGLPC = inpcv - ALLOCATE (RCGLPC(NVCGLPC)) - RCGLPC(:) = MIN(R_MAX_CNUM_PC,zlam(1:NVCGLPC)) + do ji=1,npcvecs + indarr(ji) = ji + enddo + +!--- straight insertion sort courtesy of numerical recipies + + do jj=2,npcvecs + za = zlam(jj) + ik = indarr(jj) + do ji=jj-1,1,-1 + if (zlam(ji)>one .and. (zlam(ji)>=za .or. za<=one)) then + ii=ji + exit + else + ii=0 + endif + zlam(ji+1) = zlam(ji) + indarr(ji+1) = indarr(ji) + enddo + zlam(ii+1) = za + indarr(ii+1) = ik + enddo + + inpcv = npcvecs + + do while (zlam(inpcv) <= zero) + if (mype==0) write(6,*)'Warning - eigenvalue less than 1: ',zlam(inpcv) + inpcv = inpcv-1 + if (inpcv == 0) then + if (mype==0) write(6,*)'SETUP_PRECOND: cannot form preconditioner - '//& + 'no positive eigenvalues.' + if (mype==0) write(6,*)'SETUP_PRECOND: minimisation will not be preconditioned.' + exit + endif + enddo + + if (inpcv>0) then + if (mype==0) write(6,*)'Number of preconditioning vectors selected is ',inpcv + if (mype==0) write(6,*)'SETUP_PRECOND: selected eigenvalues are: ',(zlam(ji),ji=1,inpcv) - ALLOCATE (YVCGLPC(NVCGLPC)) - DO jj=1,NVCGLPC - CALL ALLOCATE_CV(YVCGLPC(jj)) - ENDDO + if (allocated(yvcglpc)) then + do jj=1,nvcglpc + call deallocate_cv(yvcglpc(jj)) + enddo + deallocate(yvcglpc) + nvcglpc=0 + endif + if (allocated(rcglpc)) deallocate(rcglpc) -!--- apply Householder transformations to the eigenvectors to get the + !--- Save eigenvalues + nvcglpc = inpcv + allocate (rcglpc(nvcglpc)) + rcglpc(:) = min(r_max_cnum_pc,zlam(1:nvcglpc)) + + allocate (yvcglpc(nvcglpc)) + do jj=1,nvcglpc + call allocate_cv(yvcglpc(jj)) + enddo + +!--- apply householder transformations to the eigenvectors to get the !--- eigenvectors of the preconditioner - DO jj=1,NVCGLPC - YVCGLPC(jj) = zero - CALL SET_CVSECTION(zuut(1:npcvecs,indarr(jj)),YVCGLPC(jj),1,npcvecs) + do jj=1,nvcglpc + yvcglpc(jj) = zero + call set_cvsection(zuut(1:npcvecs,indarr(jj)),yvcglpc(jj),1,npcvecs) - do ji=npcvecs,1,-1 - zps = DOT_PRODUCT (yvcglwk(ji),YVCGLPC(jj)) - do jk=1,YVCGLPC(jj)%lencv - YVCGLPC(jj)%values(jk) = YVCGLPC(jj)%values(jk) - zps*yvcglwk(ji)%values(jk) - enddo - ENDDO - ENDDO - LMPCGL = .true. - ELSE - NVCGLPC = 0 - LMPCGL = .false. - ENDIF - - deallocate(indarr) - deallocate(zq,zlam,zU,zUUT,zwork,zzz) - endif - NPCVECS = 0 - DO jj=1,npcvecs - CALL DEALLOCATE_CV(YVCGLWK(jj)) - ENDDO - DEALLOCATE(YVCGLWK) + do ji=npcvecs,1,-1 + zps = dot_product (yvcglwk(ji),yvcglpc(jj)) + do jk=1,yvcglpc(jj)%lencv + yvcglpc(jj)%values(jk) = yvcglpc(jj)%values(jk) - zps*yvcglwk(ji)%values(jk) + enddo + enddo + enddo + lmpcgl = .true. + else + nvcglpc = 0 + lmpcgl = .false. + endif + + deallocate(indarr) + deallocate(zq,zlam,zu,zuut,zwork,zzz) +endif +npcvecs = 0 +do jj=1,npcvecs + call deallocate_cv(yvcglwk(jj)) +enddo +deallocate(yvcglwk) - return +return end subroutine setup_precond ! ------------------------------------------------------------------------------ -! PRECOND - Preconditioner for minimization +! precond - Preconditioner for minimization ! ------------------------------------------------------------------------------ subroutine lanczos_precond(ycvx,kmat) !$$$ subprogram documentation block @@ -1593,37 +1593,37 @@ subroutine lanczos_precond(ycvx,kmat) ! !$$$ end documentation block -IMPLICIT NONE +implicit none -TYPE(CONTROL_VECTOR),INTENT(INOUT) :: ycvx -INTEGER(i_kind) ,INTENT(IN ) :: kmat +type(control_vector),intent(inout) :: ycvx +integer(i_kind) ,intent(in ) :: kmat -REAL(r_kind) :: zevals(NVCGLPC),zdp(NVCGLPC) -INTEGER(i_kind) :: jk, ji +real(r_kind) :: zevals(nvcglpc),zdp(nvcglpc) +integer(i_kind) :: jk, ji if (kmat== 1 ) then - zevals(:) = RCGLPC(:) -ELSEIF (kmat==-1 ) then - zevals(:) = one/RCGLPC(:) -ELSEIF (kmat== 2) then - zevals(1:NVCGLPC) = sqrt(RCGLPC(:)) -ELSEIF (kmat==-2) then - zevals(1:NVCGLPC) = one/sqrt(RCGLPC(:)) + zevals(:) = rcglpc(:) +elseif (kmat==-1 ) then + zevals(:) = one/rcglpc(:) +elseif (kmat== 2) then + zevals(1:nvcglpc) = sqrt(rcglpc(:)) +elseif (kmat==-2) then + zevals(1:nvcglpc) = one/sqrt(rcglpc(:)) else write(6,*)'Error: invalid value for kmat in precond: ',kmat write(6,*)'PRECOND: invalid value for kmat' call stop2(327) endif -do jk=1,NVCGLPC - zdp(jk) = (zevals(jk)-one)*DOT_PRODUCT(ycvx,YVCGLPC(jk)) +do jk=1,nvcglpc + zdp(jk) = (zevals(jk)-one)*dot_product(ycvx,yvcglpc(jk)) enddo -DO jk=1,NVCGLPC - DO ji=1,ycvx%lencv - ycvx%values(ji) = ycvx%values(ji) + YVCGLPC(jk)%values(ji) * zdp(jk) - ENDDO -ENDDO +do jk=1,nvcglpc + do ji=1,ycvx%lencv + ycvx%values(ji) = ycvx%values(ji) + yvcglpc(jk)%values(ji) * zdp(jk) + enddo +enddo return end subroutine lanczos_precond @@ -1651,7 +1651,7 @@ subroutine read_lanczos(kmaxit) ! !$$$ end documentation block -IMPLICIT NONE +implicit none integer(i_kind) , intent(inout) :: kmaxit @@ -1673,7 +1673,7 @@ subroutine read_lanczos(kmaxit) if (mype==0) then iunit=get_lun() clfile='zlanczos.XXX' - WRITE(clfile(10:12),'(I3.3)') jiter + write(clfile(10:12),'(I3.3)') jiter write(6,*)'Reading Lanczos coef. from file ',clfile open(iunit,file=trim(clfile),form='unformatted') diff --git a/src/gsi/lightinfo.f90 b/src/gsi/lightinfo.f90 index b0ebcdacfd..7c1c63732c 100755 --- a/src/gsi/lightinfo.f90 +++ b/src/gsi/lightinfo.f90 @@ -5,12 +5,12 @@ module lightinfo ! prgmmr: Apodaca org: CSU/CIRA date: 2015-08-11 ! abstract: This module contains variables related to the ! direct assimilation of lightning observations -! (e.g. GOES-16 GLM). +! (e.g. goes-16 glm). ! ! program history log: ! ! . . . . . . . . -! 2016-05-03 Apodaca - updates regarding GLM observation errors +! 2016-05-03 Apodaca - updates regarding glm observation errors ! ! Subroutines Included: ! sub init_light - initialize lightning related variables to defaults @@ -27,8 +27,8 @@ module lightinfo ! def gross_light - gross error for lightning obs ! def glermax - gross error parameter - max error ! def glermin - gross error parameter - min error -! def b_light - b value for variational QC -! def pg_light - pg value for variational QC +! def b_light - b value for variational qc +! def pg_light - pg value for variational qc ! def nulight - satellite/instrument ! def iuse_light - use to turn off lightning data ! @@ -70,7 +70,7 @@ subroutine init_light ! assimilation routines ! ! program history log: -! 2016-05-03 apodaca, updates regarding GLM observation errors +! 2016-05-03 apodaca, updates regarding glm observation errors ! 2016-06-13 apodaca, documentation ! ! input argument list: @@ -91,7 +91,7 @@ subroutine init_light deltiml = r1200 ! model timestep diag_light =.true. ! flag to toggle creation of lightning diagnostic file mype_light = 0 ! task to print light info. Note that mype_light - ! MUST equal mype_rad (see radinfo.f90) in order for + ! must equal mype_rad (see radinfo.f90) in order for ! statspcp.f90 to print out the correct information end subroutine init_light @@ -106,7 +106,7 @@ subroutine lightinfo_read ! ! program history log: ! 2015-08-14 apodaca - original code based on pcpinfo and aeroinfo -! 2016-05-03 apodaca - updates regarding GLM observation errors +! 2016-05-03 apodaca - updates regarding glm observation errors ! input argument list: ! mype - mpi task id ! @@ -180,7 +180,7 @@ subroutine lightinfo_read rewind(lunin) !---------------------------------------------------------- -! READ INFO FILE +! Read info file !---------------------------------------------------------- j=0 do k=1,nlines diff --git a/src/gsi/m_aeroNode.F90 b/src/gsi/m_aeroNode.F90 deleted file mode 100644 index b1dfe70050..0000000000 --- a/src/gsi/m_aeroNode.F90 +++ /dev/null @@ -1,329 +0,0 @@ -module m_aeroNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_aeroNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type aeroNode (aerosal AOD) -! -! program history log: -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! -! input argument list: see Fortran 90 style document below -! -! output argument list: see Fortran 90 style document below -! -! attributes: -! language: Fortran 90 and/or above -! machine: -! -!$$$ end subprogram documentation block - -! module interface: - use m_obsdiagNode, only: obs_diag,aofp_obs_diag => fptr_obsdiagNode - use m_obsdiagNode, only: obs_diags - - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: aeroNode - - type,extends(obsNode):: aeroNode - !type(aero_ob_type),pointer :: llpoint => NULL() - type(aofp_obs_diag), dimension(:), pointer :: diags => NULL() - real(r_kind),dimension(:),pointer :: res => NULL() ! aerosol property residual - real(r_kind),dimension(:),pointer :: err2 => NULL() ! aerosol property error squared - real(r_kind),dimension(:),pointer :: raterr2 => NULL() ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time in sec - real(r_kind) :: wij(4) =0._r_kind ! horizontal interpolation weights - real(r_kind),dimension(:,:),pointer :: daod_dvar => NULL()! jacobians_aero (nsig*n_aerosols,nchan) - real(r_kind),dimension(:),pointer :: prs => NULL() ! pressure levels - integer(i_kind),dimension(:),pointer :: ipos => NULL() - integer(i_kind),dimension(:),pointer :: icx => NULL() - integer(i_kind) :: ij(4) =0 ! horizontal locations - integer(i_kind) :: nlaero =0 ! number of channels - !logical :: luse ! flag indicating if ob is used in pen. - !integer(i_kind) :: idv,iob ! device id and obs index for sorting - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution - integer(i_kind),dimension(:),pointer :: ich => NULL() - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: getTLDdp => getTLDdp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - procedure:: clean => obsNode_clean_ - end type aeroNode - - public:: aeroNode_typecast - public:: aeroNode_nextcast - interface aeroNode_typecast; module procedure typecast_ ; end interface - interface aeroNode_nextcast; module procedure nextcast_ ; end interface - - public:: aeroNode_appendto - interface aeroNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_aeroNode" - -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(aeroNode) - use m_obsNode, only: obsNode - implicit none - type(aeroNode),pointer:: ptr_ - class(obsNode),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer. - select type(aNode) - type is(aeroNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(aeroNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type(aeroNode),pointer:: ptr_ - class(obsNode),target ,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) -!-- append aNode to linked-list oLL - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(aeroNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[aeroNode]" -end function mytype - -subroutine obsNode_clean_(aNode) - implicit none - class(aeroNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_clean_' -_ENTRY_(myname_) -!_TRACEV_(myname_,'%mytype() =',aNode%mytype()) - if(associated(aNode%daod_dvar)) deallocate(aNode%daod_dvar) - if(associated(aNode%diags )) deallocate(aNode%diags ) - if(associated(aNode%res )) deallocate(aNode%res ) - if(associated(aNode%err2 )) deallocate(aNode%err2 ) - if(associated(aNode%raterr2)) deallocate(aNode%raterr2) - if(associated(aNode%prs )) deallocate(aNode%prs ) - if(associated(aNode%ipos )) deallocate(aNode%ipos ) - if(associated(aNode%icx )) deallocate(aNode%icx ) -_EXIT_(myname_) -return -end subroutine obsNode_clean_ - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use aeroinfo, only: nsigaerojac - use m_obsdiagNode, only: obsdiagLookup_locate - implicit none - class(aeroNode), intent(inout):: aNode - integer(i_kind), intent(in ):: iunit - integer(i_kind), intent( out):: istat - type(obs_diags), intent(in ):: diagLookup - logical,optional,intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xread_' - integer(i_kind),allocatable,dimension(:):: ich - integer(i_kind):: nlaero,nlevp,k - logical:: skip_ -_ENTRY_(myname_) - - skip_=.false. - if(present(skip)) skip_=skip - - if(skip_) then - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%nlaero), iostat =',istat) - _EXIT_(myname_) - return - endif - - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(ich,...)), iostat =',istat) - _EXIT_(myname_) - return - endif - - else - read(iunit,iostat=istat) aNode%nlaero - if (istat/=0) then - call perr(myname_,'read(%nlaero), iostat =',istat) - _EXIT_(myname_) - return - end if - - if(associated(aNode%daod_dvar)) deallocate(aNode%daod_dvar) - if(associated(aNode%diags )) deallocate(aNode%diags ) - if(associated(aNode%res )) deallocate(aNode%res ) - if(associated(aNode%err2 )) deallocate(aNode%err2 ) - if(associated(aNode%raterr2)) deallocate(aNode%raterr2) - if(associated(aNode%prs )) deallocate(aNode%prs ) - if(associated(aNode%ipos )) deallocate(aNode%ipos ) - if(associated(aNode%icx )) deallocate(aNode%icx ) - - nlaero=aNode%nlaero - nlevp = max(nlaero,1) - - allocate( anode%daod_dvar(nsigaerojac,nlaero), & - aNode%diags(nlaero), & - aNode%res (nlaero), & - aNode%err2 (nlaero), & - aNode%raterr2(nlaero), & - aNode%prs (nlevp ), & - aNode%ipos (nlaero), & - aNode%icx (nlaero) ) - - allocate( ich (nlaero) ) - - read(iunit,iostat=istat) ich , & - aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%prs , & - aNode%ipos , & - aNode%icx , & - aNode%daod_dvar, & - aNode%wij , & - aNode%ij - if (istat/=0) then - call perr(myname_,'read(ich,%(...)), iostat =',istat) - _EXIT_(myname_) - return - end if - - do k=1,aNode%nlaero - aNode%diags(k)%ptr => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,ich(k)) - if(.not. associated(aNode%diags(k)%ptr)) then - call perr(myname_,'obsdiagLookup_locate(k), k =',k) - call perr(myname_,' %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call perr(myname_,' ich =',ich(k)) - call die(myname_) - endif - enddo - deallocate(ich) - - endif ! if(skip_); else -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(aeroNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xwrite_' - integer(i_kind):: k -_ENTRY_(myname_) - - write(junit,iostat=jstat) aNode%nlaero - if(jstat/=0) then - call perr(myname_,'write(%nlaero), iostat =',jstat) - _EXIT_(myname_) - return - endif - - write(junit,iostat=jstat) (/(aNode%diags(k)%ptr%ich,k=1,aNode%nlaero)/), & !(nlaero) - aNode%res , & !(nlaero) - aNode%err2 , & !(nlaero) - aNode%raterr2, & !(nlaero) - aNode%prs , & !(nlevp ) - aNode%ipos , & !(nlaero) - aNode%icx , & !(nlaero) - aNode%daod_dvar, & !(nsigaerojac,nlaero) - aNode%wij , & - aNode%ij - if (jstat/=0) then - call perr(myname_,'write(ich,%(...)), iostat =',jstat) - _EXIT_(myname_) - return - end if -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - use m_cvgridLookup, only: cvgridLookup_getiw - implicit none - class(aeroNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,aNode%wij) -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(aeroNode),intent(in):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' - integer(i_kind):: k -_ENTRY_(myname_) - isvalid_= all ( (/ (associated(aNode%diags(k)%ptr), k=1,aNode%nlaero) /) ) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine getTLDdp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(aeroNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - integer(kind=i_kind):: k - do k=1,aNode%nlaero - tlddp = tlddp + aNode%diags(k)%ptr%tldepart(jiter)*aNode%diags(k)%ptr%tldepart(jiter) - enddo - if(present(nob)) nob=nob+aNode%nlaero -return -end subroutine getTLDdp_ - -end module m_aeroNode diff --git a/src/gsi/m_aerolNode.F90 b/src/gsi/m_aerolNode.F90 deleted file mode 100644 index 8ad15bd026..0000000000 --- a/src/gsi/m_aerolNode.F90 +++ /dev/null @@ -1,239 +0,0 @@ -module m_aerolNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_aerolNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type aerolNode (unfinished leveled aerosol data?) -! -! program history log: -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! -! input argument list: see Fortran 90 style document below -! -! output argument list: see Fortran 90 style document below -! -! attributes: -! language: Fortran 90 and/or above -! machine: -! -!$$$ end subprogram documentation block - -! module interface: - - use m_obsdiagNode, only: obs_diag - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: aerolNode - - type,extends(obsNode):: aerolNode - !type(aerol_ob_type),pointer :: llpoint => NULL() - type(obs_diag), pointer :: diags => NULL() - real(r_kind) :: res =0._r_kind ! aerosol residual - real(r_kind) :: err2 =0._r_kind ! aerosol obs error squared - real(r_kind) :: raterr2=0._r_kind ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time - real(r_kind) :: b =0._r_kind ! variational quality control parameter - real(r_kind) :: pg =0._r_kind ! variational quality control parameter - real(r_kind) :: wij(8) =0._r_kind ! horizontal interpolation weights - integer(i_kind) :: ij(8) =0_i_kind ! horizontal locations - !logical :: luse ! flag indicating if ob is used in pen. - - !integer(i_kind) :: idv,iob ! device id and obs index for sorting - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution - real (r_kind) :: dlev =0._r_kind ! reference to the vertical grid - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: getTLDdp => getTLDdp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - ! procedure:: clean => obsNode_clean_ - end type aerolNode - - public:: aerolNode_typecast - public:: aerolNode_nextcast - interface aerolNode_typecast; module procedure typecast_ ; end interface - interface aerolNode_nextcast; module procedure nextcast_ ; end interface - - public:: aerolNode_appendto - interface aerolNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_aerolNode" - -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(aerolNode) - use m_obsNode, only: obsNode - implicit none - type(aerolNode),pointer:: ptr_ - class(obsNode ),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer. - select type(aNode) - type is(aerolNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(aerolNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type(aerolNode),pointer:: ptr_ - class(obsNode ),target ,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) -!-- append aNode to linked-list oLL - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(aerolNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[aerolNode]" -end function mytype - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use m_obsdiagNode, only: obsdiagLookup_locate - use m_obsdiagNode, only: obs_diags - implicit none - class(aerolNode), intent(inout):: aNode - integer(i_kind), intent(in ):: iunit - integer(i_kind), intent( out):: istat - type(obs_diags), intent(in ):: diagLookup - logical,optional,intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xread_' - logical:: skip_ -_ENTRY_(myname_) - - skip_=.false. - if(present(skip)) skip_=skip - if(skip_) then - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(res,...)), istat =',istat) - _EXIT_(myname_) - return - endif - else - read(iunit,iostat=istat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%dlev , & - aNode%wij , & - aNode%ij - if(istat/=0) then - call perr(myname_,'read(%(res,...)), istat =',istat) - _EXIT_(myname_) - return - endif - - aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) - if(.not. associated(aNode%diags)) then - call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call die(myname_) - endif - endif -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(aerolNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xwrite_' -_ENTRY_(myname_) - - write(junit,iostat=jstat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%dlev , & - aNode%wij , & - aNode%ij -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - use m_cvgridLookup, only: cvgridLookup_getiw - implicit none - class(aerolNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%dlev,aNode%ij,aNode%wij) -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(aerolNode),intent(in):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' -_ENTRY_(myname_) - isvalid_=associated(aNode%diags) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine getTLDdp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(aerolNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) - if(present(nob)) nob=nob+1 -return -end subroutine getTLDdp_ - -end module m_aerolNode diff --git a/src/gsi/m_aerolnode.F90 b/src/gsi/m_aerolnode.F90 new file mode 100644 index 0000000000..87573987cd --- /dev/null +++ b/src/gsi/m_aerolnode.F90 @@ -0,0 +1,239 @@ +module m_aerolnode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_aerolnode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type aerolnode (unfinished leveled aerosol data?) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + + use m_obsdiagnode, only: obs_diag + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: aerolnode + + type,extends(obsnode):: aerolnode + !type(aerol_ob_type),pointer :: llpoint => null() + type(obs_diag), pointer :: diags => null() + real(r_kind) :: res =0._r_kind ! aerosol residual + real(r_kind) :: err2 =0._r_kind ! aerosol obs error squared + real(r_kind) :: raterr2=0._r_kind ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time + real(r_kind) :: b =0._r_kind ! variational quality control parameter + real(r_kind) :: pg =0._r_kind ! variational quality control parameter + real(r_kind) :: wij(8) =0._r_kind ! horizontal interpolation weights + integer(i_kind) :: ij(8) =0_i_kind ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + real (r_kind) :: dlev =0._r_kind ! reference to the vertical grid + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + ! procedure:: clean => obsnode_clean_ + end type aerolnode + + public:: aerolnode_typecast + public:: aerolnode_nextcast + interface aerolnode_typecast; module procedure typecast_ ; end interface + interface aerolnode_nextcast; module procedure nextcast_ ; end interface + + public:: aerolnode_appendto + interface aerolnode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_aerolnode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(aerolnode) + use m_obsnode, only: obsnode + implicit none + type(aerolnode),pointer:: ptr_ + class(obsnode ),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(anode) + type is(aerolnode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(aerolnode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type(aerolnode),pointer:: ptr_ + class(obsnode ),target ,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) +!-- append anode to linked-list oll + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(aerolnode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[aerolnode]" +end function mytype + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use m_obsdiagnode, only: obsdiaglookup_locate + use m_obsdiagnode, only: obs_diags + implicit none + class(aerolnode), intent(inout):: anode + integer(i_kind), intent(in ):: iunit + integer(i_kind), intent( out):: istat + type(obs_diags), intent(in ):: diaglookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'::obsnode_xread_' + logical:: skip_ +_ENTRY_(myname_) + + skip_=.false. + if(present(skip)) skip_=skip + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,...)), istat =',istat) + _EXIT_(myname_) + return + endif + else + read(iunit,iostat=istat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%dlev , & + anode%wij , & + anode%ij + if(istat/=0) then + call perr(myname_,'read(%(res,...)), istat =',istat) + _EXIT_(myname_) + return + endif + + anode%diags => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,1) + if(.not. associated(anode%diags)) then + call perr(myname_,'obsdiaglookup_locate(), %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(aerolnode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'::obsnode_xwrite_' +_ENTRY_(myname_) + + write(junit,iostat=jstat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%dlev , & + anode%wij , & + anode%ij +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + use m_cvgridlookup, only: cvgridlookup_getiw + implicit none + class(aerolnode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_sethop_' +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%dlev,anode%ij,anode%wij) +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(aerolnode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_isvalid_' +_ENTRY_(myname_) + isvalid_=associated(anode%diags) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(aerolnode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + anode%diags%tldepart(jiter)*anode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 +return +end subroutine gettlddp_ + +end module m_aerolnode diff --git a/src/gsi/m_aeronode.F90 b/src/gsi/m_aeronode.F90 new file mode 100644 index 0000000000..9ac0b2584b --- /dev/null +++ b/src/gsi/m_aeronode.F90 @@ -0,0 +1,329 @@ +module m_aeronode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_aeronode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type aeronode (aerosal aod) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + use m_obsdiagnode, only: obs_diag,aofp_obs_diag => fptr_obsdiagnode + use m_obsdiagnode, only: obs_diags + + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: aeronode + + type,extends(obsnode):: aeronode + !type(aero_ob_type),pointer :: llpoint => null() + type(aofp_obs_diag), dimension(:), pointer :: diags => null() + real(r_kind),dimension(:),pointer :: res => null() ! aerosol property residual + real(r_kind),dimension(:),pointer :: err2 => null() ! aerosol property error squared + real(r_kind),dimension(:),pointer :: raterr2 => null() ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: wij(4) =0._r_kind ! horizontal interpolation weights + real(r_kind),dimension(:,:),pointer :: daod_dvar => null()! jacobians_aero (nsig*n_aerosols,nchan) + real(r_kind),dimension(:),pointer :: prs => null() ! pressure levels + integer(i_kind),dimension(:),pointer :: ipos => null() + integer(i_kind),dimension(:),pointer :: icx => null() + integer(i_kind) :: ij(4) =0 ! horizontal locations + integer(i_kind) :: nlaero =0 ! number of channels + !logical :: luse ! flag indicating if ob is used in pen. + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + integer(i_kind),dimension(:),pointer :: ich => null() + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + procedure:: clean => obsnode_clean_ + end type aeronode + + public:: aeronode_typecast + public:: aeronode_nextcast + interface aeronode_typecast; module procedure typecast_ ; end interface + interface aeronode_nextcast; module procedure nextcast_ ; end interface + + public:: aeronode_appendto + interface aeronode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_aeronode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(aeronode) + use m_obsnode, only: obsnode + implicit none + type(aeronode),pointer:: ptr_ + class(obsnode),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(anode) + type is(aeronode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(aeronode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type(aeronode),pointer:: ptr_ + class(obsnode),target ,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) +!-- append anode to linked-list oll + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(aeronode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[aeronode]" +end function mytype + +subroutine obsnode_clean_(anode) + implicit none + class(aeronode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_clean_' +_ENTRY_(myname_) +!_TRACEV_(myname_,'%mytype() =',anode%mytype()) + if(associated(anode%daod_dvar)) deallocate(anode%daod_dvar) + if(associated(anode%diags )) deallocate(anode%diags ) + if(associated(anode%res )) deallocate(anode%res ) + if(associated(anode%err2 )) deallocate(anode%err2 ) + if(associated(anode%raterr2)) deallocate(anode%raterr2) + if(associated(anode%prs )) deallocate(anode%prs ) + if(associated(anode%ipos )) deallocate(anode%ipos ) + if(associated(anode%icx )) deallocate(anode%icx ) +_EXIT_(myname_) + return +end subroutine obsnode_clean_ + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use aeroinfo, only: nsigaerojac + use m_obsdiagnode, only: obsdiaglookup_locate + implicit none + class(aeronode), intent(inout):: anode + integer(i_kind), intent(in ):: iunit + integer(i_kind), intent( out):: istat + type(obs_diags), intent(in ):: diaglookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'::obsnode_xread_' + integer(i_kind),allocatable,dimension(:):: ich + integer(i_kind):: nlaero,nlevp,k + logical:: skip_ +_ENTRY_(myname_) + + skip_=.false. + if(present(skip)) skip_=skip + + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%nlaero), iostat =',istat) + _EXIT_(myname_) + return + endif + + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(ich,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) anode%nlaero + if (istat/=0) then + call perr(myname_,'read(%nlaero), iostat =',istat) + _EXIT_(myname_) + return + end if + + if(associated(anode%daod_dvar)) deallocate(anode%daod_dvar) + if(associated(anode%diags )) deallocate(anode%diags ) + if(associated(anode%res )) deallocate(anode%res ) + if(associated(anode%err2 )) deallocate(anode%err2 ) + if(associated(anode%raterr2)) deallocate(anode%raterr2) + if(associated(anode%prs )) deallocate(anode%prs ) + if(associated(anode%ipos )) deallocate(anode%ipos ) + if(associated(anode%icx )) deallocate(anode%icx ) + + nlaero=anode%nlaero + nlevp = max(nlaero,1) + + allocate( anode%daod_dvar(nsigaerojac,nlaero), & + anode%diags(nlaero), & + anode%res (nlaero), & + anode%err2 (nlaero), & + anode%raterr2(nlaero), & + anode%prs (nlevp ), & + anode%ipos (nlaero), & + anode%icx (nlaero) ) + + allocate( ich (nlaero) ) + + read(iunit,iostat=istat) ich , & + anode%res , & + anode%err2 , & + anode%raterr2, & + anode%prs , & + anode%ipos , & + anode%icx , & + anode%daod_dvar, & + anode%wij , & + anode%ij + if (istat/=0) then + call perr(myname_,'read(ich,%(...)), iostat =',istat) + _EXIT_(myname_) + return + end if + + do k=1,anode%nlaero + anode%diags(k)%ptr => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,ich(k)) + if(.not. associated(anode%diags(k)%ptr)) then + call perr(myname_,'obsdiaglookup_locate(k), k =',k) + call perr(myname_,' %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call perr(myname_,' ich =',ich(k)) + call die(myname_) + endif + enddo + deallocate(ich) + + endif ! if(skip_); else +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(aeronode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'::obsnode_xwrite_' + integer(i_kind):: k +_ENTRY_(myname_) + + write(junit,iostat=jstat) anode%nlaero + if(jstat/=0) then + call perr(myname_,'write(%nlaero), iostat =',jstat) + _EXIT_(myname_) + return + endif + + write(junit,iostat=jstat) (/(anode%diags(k)%ptr%ich,k=1,anode%nlaero)/), & !(nlaero) + anode%res , & !(nlaero) + anode%err2 , & !(nlaero) + anode%raterr2, & !(nlaero) + anode%prs , & !(nlevp ) + anode%ipos , & !(nlaero) + anode%icx , & !(nlaero) + anode%daod_dvar, & !(nsigaerojac,nlaero) + anode%wij , & + anode%ij + if (jstat/=0) then + call perr(myname_,'write(ich,%(...)), iostat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + use m_cvgridlookup, only: cvgridlookup_getiw + implicit none + class(aeronode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_sethop_' +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%ij,anode%wij) +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(aeronode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_isvalid_' + integer(i_kind):: k +_ENTRY_(myname_) + isvalid_= all ( (/ (associated(anode%diags(k)%ptr), k=1,anode%nlaero) /) ) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(aeronode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + integer(kind=i_kind):: k + do k=1,anode%nlaero + tlddp = tlddp + anode%diags(k)%ptr%tldepart(jiter)*anode%diags(k)%ptr%tldepart(jiter) + enddo + if(present(nob)) nob=nob+anode%nlaero + return +end subroutine gettlddp_ + +end module m_aeronode diff --git a/src/gsi/m_berror_stats.f90 b/src/gsi/m_berror_stats.f90 index 82d4cf0d6c..5c15aff20f 100644 --- a/src/gsi/m_berror_stats.f90 +++ b/src/gsi/m_berror_stats.f90 @@ -129,7 +129,7 @@ end subroutine get_dims ! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS ! !BOP ------------------------------------------------------------------- ! -! !IROUTINE: lset - set (logical) parameter options internal to B +! !IROUTINE: lset - set (logical) parameter options internal to b ! ! !DESCRIPTION: ! @@ -184,7 +184,7 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit) real(r_single),dimension(nlat,nsig,nsig),intent( out) :: agvin real(r_single),dimension(nlat,nsig) ,intent( out) :: bvin,wgvin real(r_single),dimension(nlat,nsig) ,intent( out) :: pputin - integer(i_kind) ,intent(in ) :: mype ! "my" processor ID + integer(i_kind) ,intent(in ) :: mype ! "my" processor id integer(i_kind),optional ,intent(in ) :: lunit ! an alternative unit ! !REVISION HISTORY: @@ -194,7 +194,7 @@ subroutine read_bal(agvin,bvin,wgvin,pputin,fut2ps,mype,lunit) ! 25Feb10 - Zhu ! - change the structure of background error file ! - read in agvin,wgvin,bvin only -! 09Oct12 - Gu add fut2ps to project unbalanced temp to surface pressure in static B modeling +! 09Oct12 - Gu add fut2ps to project unbalanced temp to surface pressure in static b modeling !EOP ___________________________________________________________________ character(len=*),parameter :: myname_=myname//'::read_bal' @@ -261,7 +261,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt use kinds,only : r_single,r_kind use gridmod,only : nlat,nlon,nsig - use radiance_mod, only: n_clouds_fwd,cloud_names_fwd + use radiance_mod, only: n_clouds_fwd,cloud_names_fwd implicit none @@ -280,7 +280,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt integer(i_kind) ,intent(in ) :: qoption integer(i_kind) ,intent(in ) :: cwoption - integer(i_kind) ,intent(in ) :: mype ! "my" processor ID + integer(i_kind) ,intent(in ) :: mype ! "my" processor id integer(i_kind),optional ,intent(in ) :: lunit ! an alternative unit ! !REVISION HISTORY: @@ -470,8 +470,8 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt enddo if (cwcoveqqcov_) then if(iq>0.and.icw>0) then - hwll(:,:,icw)=hwll(:,:,iq) - vz (:,:,icw)=vz (:,:,iq) + hwll(:,:,icw)=hwll(:,:,iq) + vz (:,:,icw)=vz (:,:,iq) end if end if if (cwoption==1 .or. cwoption==3) then @@ -508,7 +508,7 @@ subroutine read_wgt(corz,corp,hwll,hwllp,vz,corsst,hsst,varq,qoption,varcw,cwopt deallocate(found3d,found2d) - return + return end subroutine read_wgt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -531,7 +531,7 @@ subroutine setcoroz_(coroz,mype) use gridmod,only: lon1,lat1 use guess_grids,only: ntguessig - use guess_grids,only: ges_prsi ! interface pressures (kPa) + use guess_grids,only: ges_prsi ! interface pressures (kpa) use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_metguess_mod,only: gsi_metguess_bundle @@ -539,18 +539,18 @@ subroutine setcoroz_(coroz,mype) implicit none real(r_single),dimension(nlat,nsig),intent( out) :: coroz ! of ozone - integer(i_kind) ,intent(in ) :: mype ! ID of this processor + integer(i_kind) ,intent(in ) :: mype ! id of this processor ! !REVISION HISTORY: ! 31Jul08 - Jing Guo -! - adopted from PREWGT of previous version +! - adopted from prewgt of previous version ! 2013-10-19 oz guess field in metguess now !EOP ___________________________________________________________________ character(len=*),parameter :: myname_=myname//'::setcoroz_' real(r_kind),parameter :: r25 = one/25.0_r_kind - real(r_kind),dimension(:,:,:),pointer :: ges_oz=>NULL() + real(r_kind),dimension(:,:,:),pointer :: ges_oz=>null() real(r_kind),dimension(nsig+1,npe) :: work_oz,work_oz1 real(r_kind),dimension(nsig) :: ozmz @@ -590,7 +590,7 @@ subroutine setcoroz_(coroz,mype) enddo enddo enddo - work_oz(nsig+1,mm1)=float(lon1*lat1) + work_oz(nsig+1,mm1)=real(lon1*lat1,r_kind) call mpi_allreduce(work_oz,work_oz1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) @@ -641,12 +641,12 @@ subroutine sethwlloz_(hwlloz,mype) use kinds, only: r_single,r_kind use mpimod, only: levs_id use gridmod, only: nnnn1o,nsig,nlon,nlat - use constants,only: two,three,pi,rearth_equator + use constants,only: two,three,four,pi,rearth_equator implicit none real(r_single),dimension(nlat,nsig),intent( out) :: hwlloz - integer(i_kind) ,intent(in ) :: mype ! ID of this processor + integer(i_kind) ,intent(in ) :: mype ! id of this processor ! !REVISION HISTORY: ! 31Jul08 - Jing Guo @@ -671,10 +671,10 @@ subroutine sethwlloz_(hwlloz,mype) if ( k1>0 ) then write(6,*) myname_,'(PREWGT): mype = ',mype, k1 if ( k1<=nsig*3/4 ) then - ! fact=1./hwl - fact=r40000/(r400*nlon) + ! fact=1._r_kind/hwl + fact=r40000/(r400*nlon) else - fact=r40000/(nlon*(r800-r400*(nsig-k1)/(nsig-nsig*3/4))) + fact=r40000/(nlon*(r800-r400*(nsig-k1)/(nsig-nsig*three/four))) endif fact=fact*three hwlloz(:,k1)=s2u/fact @@ -741,7 +741,7 @@ subroutine setcorchem_(cname,corchem,rc) use gridmod,only: lon1,lat1 use guess_grids,only: ntguessig - use guess_grids,only: ges_prsi ! interface pressures (kPa) + use guess_grids,only: ges_prsi ! interface pressures (kpa) use gsi_chemguess_mod,only: gsi_chemguess_bundle use gsi_bundlemod, only: gsi_bundlegetpointer @@ -774,10 +774,10 @@ subroutine setcorchem_(cname,corchem,rc) ! sanity check if ( mype==0 ) write(6,*) myname_,'(PREWGT): mype = ',mype - ! Get information for how to use CO2 + ! Get information for how to use co2 iptr=-1 if ( size(gsi_chemguess_bundle)>0 ) then ! check to see if bundle's allocated - call gsi_bundlegetpointer(gsi_chemguess_bundle(1),cname,iptr,ierror) + call gsi_bundlegetpointer(gsi_chemguess_bundle(1),cname,iptr,ierror) if ( ierror/=0 ) then rc=-2 ! field not found return @@ -808,12 +808,12 @@ subroutine setcorchem_(cname,corchem,rc) do i=2,lat1+1 work_chem(k,mm1) = work_chem(k,mm1) + gsi_chemguess_bundle(ntguessig)%r3(iptr)%q(i,j,k)* & (ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) - !_RT not sure yet how to handle scaling factor (rozcon) in general - !_RT rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) + !not sure yet how to handle scaling factor (rozcon) in general + ! rozcon*(ges_prsi(i,j,k,ntguessig)-ges_prsi(i,j,k+1,ntguessig)) enddo enddo enddo - work_chem(nsig+1,mm1)=float(lon1*lat1) + work_chem(nsig+1,mm1)=real(lon1*lat1,r_kind) call mpi_allreduce(work_chem,work_chem1,(nsig+1)*npe,mpi_rtype,mpi_sum,& mpi_comm_world,ierror) diff --git a/src/gsi/m_berror_stats_reg.f90 b/src/gsi/m_berror_stats_reg.f90 index 896a43da4a..5b92d7fded 100644 --- a/src/gsi/m_berror_stats_reg.f90 +++ b/src/gsi/m_berror_stats_reg.f90 @@ -8,25 +8,25 @@ ! ! !INTERFACE: - module m_berror_stats_reg - use kinds,only : i_kind,r_kind - use constants, only: zero,one,max_varname_length - use gridmod, only: nsig - use chemmod, only : berror_chem,upper2lower,lower2upper - use m_berror_stats, only: berror_stats +module m_berror_stats_reg + use kinds,only : i_kind,r_kind + use constants, only: zero,one,max_varname_length + use gridmod, only: nsig + use chemmod, only : berror_chem,upper2lower,lower2upper + use m_berror_stats, only: berror_stats - implicit none + implicit none - private ! except + private ! except - ! reconfigurable parameters, via NAMELIST/setup/ - public :: berror_stats ! reconfigurable filename + ! reconfigurable parameters, via namelist/setup/ + public :: berror_stats ! reconfigurable filename ! interfaces to file berror_stats. - public :: berror_set_reg ! set internal parameters - public :: berror_get_dims_reg ! get dimensions, jfunc::createj_func() - public :: berror_read_bal_reg ! get cross-cov.stats., balmod::prebal() - public :: berror_read_wgt_reg ! get auto-cov.stats., prewgt() + public :: berror_set_reg ! set internal parameters + public :: berror_get_dims_reg ! get dimensions, jfunc::createj_func() + public :: berror_read_bal_reg ! get cross-cov.stats., balmod::prebal() + public :: berror_read_wgt_reg ! get auto-cov.stats., prewgt() ! !REVISION HISTORY: ! 25Feb10 - Zhu - adopt code format from m_berror_stats @@ -50,7 +50,7 @@ module m_berror_stats_reg ! be defined for this module in the /setup/ namelist. integer(i_kind),parameter :: default_unit_ = 22 - integer(i_kind),parameter :: ERRCODE=2 + integer(i_kind),parameter :: errcode=2 logical,save:: cwcoveqqcov_ integer(i_kind),allocatable,dimension(:):: lsig @@ -67,13 +67,13 @@ module m_berror_stats_reg ! ! !INTERFACE: - subroutine berror_get_dims_reg(msig,mlat,unit) +subroutine berror_get_dims_reg(msig,mlat,unit) - implicit none + implicit none - integer(i_kind) ,intent( out) :: msig ! dimension of levels - integer(i_kind) ,intent( out) :: mlat ! dimension of latitudes - integer(i_kind),optional,intent(in ) :: unit ! logical unit [22] + integer(i_kind) ,intent( out) :: msig ! dimension of levels + integer(i_kind) ,intent( out) :: mlat ! dimension of latitudes + integer(i_kind),optional,intent(in ) :: unit ! logical unit [22] ! !REVISION HISTORY: !EOP ___________________________________________________________________ @@ -94,17 +94,17 @@ end subroutine berror_get_dims_reg ! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS ! !BOP ------------------------------------------------------------------- ! -! !IROUTINE: berror_set_reg - set (logical) parameter options internal to B +! !IROUTINE: berror_set_reg - set (logical) parameter options internal to b ! ! !DESCRIPTION: ! ! !INTERFACE: - subroutine berror_set_reg(opt,value) +subroutine berror_set_reg(opt,value) - implicit none - character(len=*),intent(in) :: opt - logical(i_kind), intent(in) :: value + implicit none + character(len=*),intent(in) :: opt + logical(i_kind), intent(in) :: value ! !REVISION HISTORY: ! 2014-10-15 - Zhu - adopted from m_berror_stat to make code structure similar @@ -134,18 +134,18 @@ end subroutine berror_set_reg ! ! !INTERFACE: - subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) - use kinds,only : r_single - use guess_grids, only: ges_psfcavg,ges_prslavg +subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) + use kinds,only : r_single + use guess_grids, only: ges_psfcavg,ges_prslavg - implicit none + implicit none - integer(i_kind) ,intent(in ) :: msig,mlat - real(r_kind),dimension(0:mlat+1,nsig,nsig),intent( out) :: agvi - real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: wgvi - real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: bvi - integer(i_kind) ,intent(in ) :: mype ! "my" processor ID - integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit + integer(i_kind) ,intent(in ) :: msig,mlat + real(r_kind),dimension(0:mlat+1,nsig,nsig),intent( out) :: agvi + real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: wgvi + real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: bvi + integer(i_kind) ,intent(in ) :: mype ! "my" processor id + integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit ! !REVISION HISTORY: ! 25Feb10 - Zhu - extracted from rdgstat_reg @@ -167,102 +167,102 @@ subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit) real(r_single),dimension(:,:,:),allocatable:: agv_avn real(r_kind),dimension(:),allocatable:: rlsigo -! Open background error statistics file - inerr=default_unit_ - if(present(unit)) inerr=unit - open(inerr,file=berror_stats,form='unformatted',status='old') - -! Read header. - rewind inerr - read(inerr) nsigstat,nlatstat - - if(mype==0) then - write(6,*) myname_,'(PREBAL_REG): get balance variables', & - '"',trim(berror_stats),'". ', & - 'mype,nsigstat,nlatstat =', & - mype,nsigstat,nlatstat - end if - - allocate ( clat_avn(mlat) ) - allocate ( sigma_avn(1:msig) ) - allocate ( rlsigo(1:msig) ) - allocate ( agv_avn(0:mlat+1,1:msig,1:msig) ) - allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) ) - -! Read background error file to get balance variables - read(inerr)clat_avn,(sigma_avn(k),k=1,msig) - read(inerr)agv_avn,bv_avn,wgv_avn - close(inerr) - -! compute vertical(pressure) interpolation index and weight - do k=1,nsig - rlsig(k)=log(ges_prslavg(k)/ges_psfcavg) - enddo - do k=1,msig - rlsigo(k)=log(sigma_avn(k)) - enddo - - allocate(lsig(nsig),coef1(nsig),coef2(nsig)) - do k=1,nsig - if(rlsig(k)>=rlsigo(1))then - m=1 - m1=2 - lsig(k)=1 - coef1(k)=one - coef2(k)=zero - else if(rlsig(k)<=rlsigo(msig))then - m=msig-1 - m1=msig - lsig(k)=msig-1 - coef1(k)=zero - coef2(k)=one - else - m_loop: do m=1,msig-1 - m1=m+1 - if((rlsig(k)<=rlsigo(m)) .and. & - (rlsig(k)>rlsigo(m1)) )then - lsig(k)=m - exit m_loop - end if - enddo m_loop - coef1(k)=(rlsigo(m1)-rlsig(k))/(rlsigo(m1)-rlsigo(m)) - coef2(k)=one-coef1(k) - if(lsig(k)==msig)then - lsig(k)=msig-1 - coef2(k)=one - coef1(k)=zero - endif - endif - end do - -! Load agv wgv bv - do k=1,nsig - m=lsig(k) - m1=m+1 - do i=1,mlat - wgvi(i,k)=wgv_avn(i,m)*coef1(k)+wgv_avn(i,m1)*coef2(k) - bvi (i,k)=bv_avn (i,m)*coef1(k)+bv_avn (i,m1)*coef2(k) - enddo - - do j=1,nsig - l=lsig(j) - l1=l+1 - do i=1,mlat - agvi(i,j,k)=(agv_avn(i,l,m)*coef1(j)+agv_avn(i,l1,m)*coef2(j))*coef1(k) & - +(agv_avn(i,l,m1)*coef1(j)+agv_avn(i,l1,m1)*coef2(j))*coef2(k) - enddo - enddo - enddo - - agvi(0,:,:)=agvi(1,:,:) - wgvi(0,:)=wgvi(1,:) - bvi(0,:)=bvi(1,:) - agvi(mlat+1,:,:)=agvi(mlat,:,:) - wgvi(mlat+1,:)=wgvi(mlat,:) - bvi(mlat+1,:)=bvi(mlat,:) +! Open background error statistics file + inerr=default_unit_ + if(present(unit)) inerr=unit + open(inerr,file=berror_stats,form='unformatted',status='old') + +! Read header. + rewind inerr + read(inerr) nsigstat,nlatstat + + if(mype==0) then + write(6,*) myname_,'(PREBAL_REG): get balance variables', & + '"',trim(berror_stats),'". ', & + 'mype,nsigstat,nlatstat =', & + mype,nsigstat,nlatstat + end if + + allocate ( clat_avn(mlat) ) + allocate ( sigma_avn(1:msig) ) + allocate ( rlsigo(1:msig) ) + allocate ( agv_avn(0:mlat+1,1:msig,1:msig) ) + allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) ) + +! Read background error file to get balance variables + read(inerr)clat_avn,(sigma_avn(k),k=1,msig) + read(inerr)agv_avn,bv_avn,wgv_avn + close(inerr) + +! compute vertical(pressure) interpolation index and weight + do k=1,nsig + rlsig(k)=log(ges_prslavg(k)/ges_psfcavg) + enddo + do k=1,msig + rlsigo(k)=log(sigma_avn(k)) + enddo + + allocate(lsig(nsig),coef1(nsig),coef2(nsig)) + do k=1,nsig + if(rlsig(k)>=rlsigo(1))then + m=1 + m1=2 + lsig(k)=1 + coef1(k)=one + coef2(k)=zero + else if(rlsig(k)<=rlsigo(msig))then + m=msig-1 + m1=msig + lsig(k)=msig-1 + coef1(k)=zero + coef2(k)=one + else + m_loop: do m=1,msig-1 + m1=m+1 + if((rlsig(k)<=rlsigo(m)) .and. & + (rlsig(k)>rlsigo(m1)))then + lsig(k)=m + exit m_loop + end if + enddo m_loop + coef1(k)=(rlsigo(m1)-rlsig(k))/(rlsigo(m1)-rlsigo(m)) + coef2(k)=one-coef1(k) + if(lsig(k)==msig)then + lsig(k)=msig-1 + coef2(k)=one + coef1(k)=zero + endif + endif + end do + +! Load agv wgv bv + do k=1,nsig + m=lsig(k) + m1=m+1 + do i=1,mlat + wgvi(i,k)=wgv_avn(i,m)*coef1(k)+wgv_avn(i,m1)*coef2(k) + bvi (i,k)=bv_avn (i,m)*coef1(k)+bv_avn (i,m1)*coef2(k) + enddo + + do j=1,nsig + l=lsig(j) + l1=l+1 + do i=1,mlat + agvi(i,j,k)=(agv_avn(i,l,m)*coef1(j)+agv_avn(i,l1,m)*coef2(j))*coef1(k) & + +(agv_avn(i,l,m1)*coef1(j)+agv_avn(i,l1,m1)*coef2(j))*coef2(k) + enddo + enddo + enddo + + agvi(0,:,:)=agvi(1,:,:) + wgvi(0,:)=wgvi(1,:) + bvi(0,:)=bvi(1,:) + agvi(mlat+1,:,:)=agvi(mlat,:,:) + wgvi(mlat+1,:)=wgvi(mlat,:) + bvi(mlat+1,:)=bvi(mlat,:) - deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo) - return + deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo) + return end subroutine berror_read_bal_reg !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS ! @@ -274,37 +274,37 @@ end subroutine berror_read_bal_reg ! ! !INTERFACE: - subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qoption,varcw,cwoption,mype,unit) +subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qoption,varcw,cwoption,mype,unit) - use kinds,only : r_single,r_kind - use gridmod,only : nsig - use gsi_io, only : verbose - use control_vectors,only: nrf,nc2d,nc3d,mvars,nvars - use control_vectors,only: cvars => nrf_var - use control_vectors,only: cvars2d,cvars3d,cvarsmd - use guess_grids, only: ges_psfcavg,ges_prslavg - use constants, only: zero,one,ten,three - use mpeu_util,only: getindex - use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd + use kinds,only : r_single,r_kind + use gridmod,only : nsig + use gsi_io, only : verbose + use control_vectors,only: nrf,nc2d,nc3d,mvars,nvars + use control_vectors,only: cvars => nrf_var + use control_vectors,only: cvars2d,cvars3d,cvarsmd + use guess_grids, only: ges_psfcavg,ges_prslavg + use constants, only: zero,one,ten,three + use mpeu_util,only: getindex + use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd - implicit none + implicit none - integer(i_kind) ,intent(in ) :: qoption - integer(i_kind) ,intent(in ) :: cwoption - integer(i_kind) ,intent(in ) :: msig,mlat - integer(i_kind) ,intent(in ) :: mype ! "my" processor ID - integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit + integer(i_kind) ,intent(in ) :: qoption + integer(i_kind) ,intent(in ) :: cwoption + integer(i_kind) ,intent(in ) :: msig,mlat + integer(i_kind) ,intent(in ) :: mype ! "my" processor id + integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit - real(r_kind),dimension(:,:,:),intent(inout):: corz - real(r_kind),dimension(:,:) ,intent(inout) :: corp + real(r_kind),dimension(:,:,:),intent(inout):: corz + real(r_kind),dimension(:,:) ,intent(inout) :: corp - real(r_kind),dimension(0:mlat+1,1:nsig,1:nc3d), intent(inout):: hwll - real(r_kind),dimension(0:mlat+1,nvars-nc3d) , intent(inout):: hwllp - real(r_kind),dimension(nsig,0:mlat+1,1:nc3d),intent(inout):: vz - real(r_kind),dimension(mlat,nsig),intent(inout)::varq - real(r_kind),dimension(mlat,nsig),intent(inout)::varcw + real(r_kind),dimension(0:mlat+1,1:nsig,1:nc3d), intent(inout):: hwll + real(r_kind),dimension(0:mlat+1,nvars-nc3d) , intent(inout):: hwllp + real(r_kind),dimension(nsig,0:mlat+1,1:nc3d),intent(inout):: vz + real(r_kind),dimension(mlat,nsig),intent(inout)::varq + real(r_kind),dimension(mlat,nsig),intent(inout)::varcw - real(r_kind),dimension(nsig),intent(out):: rlsig + real(r_kind),dimension(nsig),intent(out):: rlsig ! !REVISION HISTORY: ! 25Feb10 - Zhu - extract from rdgstat_reg @@ -490,17 +490,17 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt end if if (isig==1) then - do n=1,nc2d - if (nrf2_loc(n)==loc) then - do i=1,mlat + do n=1,nc2d + if (nrf2_loc(n)==loc) then + do i=1,mlat corp(i,n)=corz_avn(i,1) hwllp(i,n)=hwll_avn(i,1) - end do - hwllp(0,n)=hwll_avn(0,1) - hwllp(mlat+1,n)=hwll_avn(mlat+1,1) - exit - end if - end do + end do + hwllp(0,n)=hwll_avn(0,1) + hwllp(mlat+1,n)=hwll_avn(mlat+1,1) + exit + end if + end do end if deallocate ( corz_avn ) @@ -583,13 +583,13 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt endif if( nrf3_dbz>0 )then - if(.not. nrf3_t>0) then - write(6,*)'not as expect,stop' - stop - endif - corz(:,:,nrf3_dbz)=10.0_r_kind - hwll(:,:,nrf3_dbz)=hwll(:,:,nrf3_t) - vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t) + if(.not. nrf3_t>0) then + write(6,*)'not as expect,stop' + stop + endif + corz(:,:,nrf3_dbz)=10.0_r_kind + hwll(:,:,nrf3_dbz)=hwll(:,:,nrf3_t) + vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t) endif if (nrf3_oz>0) then @@ -651,15 +651,15 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt if (nrf3_sfwter>0) then - if(mype==0) write(6,*)'Replace default with appropriate statistics for variable sfwter' - corz(1:mlat,1:nsig,nrf3_sfwter)=corz(1:mlat,1:nsig,nrf3_sf) - hwll(0:mlat+1,1:nsig,nrf3_sfwter)=hwll(0:mlat+1,1:nsig,nrf3_sf) + if(mype==0) write(6,*)'Replace default with appropriate statistics for variable sfwter' + corz(1:mlat,1:nsig,nrf3_sfwter)=corz(1:mlat,1:nsig,nrf3_sf) + hwll(0:mlat+1,1:nsig,nrf3_sfwter)=hwll(0:mlat+1,1:nsig,nrf3_sf) endif if (nrf3_vpwter>0) then - if(mype==0) write(6,*)'Replace default with appropriate statistics for variable vpwter' - corz(1:mlat,1:nsig,nrf3_vpwter)=corz(1:mlat,1:nsig,nrf3_vp) - hwll(0:mlat+1,1:nsig,nrf3_vpwter)=hwll(0:mlat+1,1:nsig,nrf3_vp) + if(mype==0) write(6,*)'Replace default with appropriate statistics for variable vpwter' + corz(1:mlat,1:nsig,nrf3_vpwter)=corz(1:mlat,1:nsig,nrf3_vp) + hwll(0:mlat+1,1:nsig,nrf3_vpwter)=hwll(0:mlat+1,1:nsig,nrf3_vp) endif ! 2d variable @@ -742,17 +742,17 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt end do end if if (n==nrf2_howv) then - call read_howv_stats(mlat,1,2,cov_dum) - do i=1,mlat - corp(i,n)=cov_dum(i,1,1) !#ww3 - hwllp(i,n) = cov_dum(i,1,2) - end do - hwllp(0,n) = hwllp(1,n) - hwllp(mlat+1,n) = hwllp(mlat,n) - - if (mype==0) print*, 'corp(i,n) = ', corp(:,n) - if (mype==0) print*, ' hwllp(i,n) = ', hwllp(:,n) -! corp(:,n)=cov_dum(:,1) + call read_howv_stats(mlat,1,2,cov_dum) + do i=1,mlat + corp(i,n)=cov_dum(i,1,1) !#ww3 + hwllp(i,n) = cov_dum(i,1,2) + end do + hwllp(0,n) = hwllp(1,n) + hwllp(mlat+1,n) = hwllp(mlat,n) + + if (mype==0) print*, 'corp(i,n) = ', corp(:,n) + if (mype==0) print*, ' hwllp(i,n) = ', hwllp(:,n) +! corp(:,n)=cov_dum(:,1) !do i=1,mlat ! corp(i,n)=0.4_r_kind !#ww3 !end do @@ -917,12 +917,12 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! ! abstract: improrts the bkerror stats for howv (variance, correlation lengths, ! etc). Each quantity has been calculated externaly, it is interpolated at -! the grid of URMA and it is saved as binary file. The exact format +! the grid of urma and it is saved as binary file. The exact format ! is shown in the code. ! The imported arrays can be 2d or 1d, aka can be lon/lat dependent or ! two one of the two (most probably on lat) ! -! For cross reference the MATLAB code is provided +! For cross reference the matlab code is provided ! Export :: ! fid = fopen(['URMA_variance_Q',num2str(i1),'.bin'],'w'); ! dum =squeeze(variance(i1,:,:)); @@ -943,7 +943,7 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) ! ! program history log: ! 2016-08-03 stelios -! 2016-08-26 stelios : Compatible with GSI. +! 2016-08-26 stelios : Compatible with gsi. ! input argument list: ! filename - The name of the file ! output argument list: @@ -967,7 +967,6 @@ subroutine read_howv_stats(nlat,nlon,npar,arrout) logical :: file_exists integer(i_kind), parameter :: lun34=34 character(len=256),dimension(npar) :: filename - integer(i_kind), parameter :: dp1 = kind(1D0) ! filename(1) = 'howv_var_berr.bin' filename(2) = 'howv_lng_berr.bin' diff --git a/src/gsi/m_cldchNode.F90 b/src/gsi/m_cldchNode.F90 deleted file mode 100644 index 068595529f..0000000000 --- a/src/gsi/m_cldchNode.F90 +++ /dev/null @@ -1,243 +0,0 @@ -module m_cldchNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_cldchNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type cldchNode (surface cldch(ibility)) -! -! program history log: -! 2015-07-10 pondeca - add cloud ceiling height (cldch) -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! -! input argument list: see Fortran 90 style document below -! -! output argument list: see Fortran 90 style document below -! -! attributes: -! language: Fortran 90 and/or above -! machine: -! -!$$$ end subprogram documentation block - -! module interface: - use m_obsdiagNode, only: obs_diag - use m_obsdiagNode, only: obs_diags - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: cldchNode - - type,extends(obsNode):: cldchNode - !type(cldch_ob_type),pointer :: llpoint => NULL() - type(obs_diag), pointer :: diags => NULL() - real(r_kind) :: res ! cldch residual - real(r_kind) :: err2 ! cldch error squared - real(r_kind) :: raterr2 ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time in sec - real(r_kind) :: b ! variational quality control parameter - real(r_kind) :: pg ! variational quality control parameter - real(r_kind) :: wij(4) ! horizontal interpolation weights - integer(i_kind) :: ij(4) ! horizontal locations - !logical :: luse ! flag indicating if ob is used in pen. - - !integer(i_kind) :: idv,iob ! device id and obs index for sorting - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: getTLDdp => getTLDdp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - ! procedure:: clean => obsNode_clean_ - end type cldchNode - - public:: cldchNode_typecast - public:: cldchNode_nextcast - interface cldchNode_typecast; module procedure typecast_ ; end interface - interface cldchNode_nextcast; module procedure nextcast_ ; end interface - - public:: cldchNode_appendto - interface cldchNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_cldchNode" - -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(cldchNode) - use m_obsNode, only: obsNode - implicit none - type(cldchNode),pointer:: ptr_ - class(obsNode ),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer. - select type(aNode) - type is(cldchNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(cldchNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type(cldchNode),pointer:: ptr_ - class(obsNode ),target,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) -!-- append aNode to linked-list oLL - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(cldchNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[cldchNode]" -end function mytype - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use m_obsdiagNode, only: obsdiagLookup_locate - implicit none - class(cldchNode), intent(inout):: aNode - integer(i_kind), intent(in ):: iunit - integer(i_kind), intent( out):: istat - type(obs_diags), intent(in ):: diagLookup - logical,optional,intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xread_' - logical:: skip_ -_ENTRY_(myname_) - - skip_=.false. - if(present(skip)) skip_=skip - - if(skip_) then - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(res,...)), iostat =',istat) - _EXIT_(myname_) - return - endif - - else - read(iunit,iostat=istat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%wij , & - aNode%ij - if (istat/=0) then - call perr(myname_,'read(%(res,...)), iostat =',istat) - _EXIT_(myname_) - return - end if - - aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) - if(.not. associated(aNode%diags)) then - call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call die(myname_) - endif - endif ! if(skip_); else -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(cldchNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xwrite_' -_ENTRY_(myname_) - - write(junit,iostat=jstat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%wij , & - aNode%ij - if (jstat/=0) then - call perr(myname_,'write(%(res,...)), iostat =',jstat) - _EXIT_(myname_) - return - end if -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - use m_cvgridLookup, only: cvgridLookup_getiw - implicit none - class(cldchNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,aNode%wij) -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(cldchNode),intent(in):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' -_ENTRY_(myname_) - isvalid_= associated(aNode%diags) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine getTLDdp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(cldchNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) - if(present(nob)) nob=nob+1 -return -end subroutine getTLDdp_ - -end module m_cldchNode diff --git a/src/gsi/m_cldchnode.F90 b/src/gsi/m_cldchnode.F90 new file mode 100644 index 0000000000..59403c80ec --- /dev/null +++ b/src/gsi/m_cldchnode.F90 @@ -0,0 +1,243 @@ +module m_cldchnode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_cldchnode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type cldchnode (surface cldch(ibility)) +! +! program history log: +! 2015-07-10 pondeca - add cloud ceiling height (cldch) +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! input argument list: see Fortran 90 style document below +! +! output argument list: see Fortran 90 style document below +! +! attributes: +! language: Fortran 90 and/or above +! machine: +! +!$$$ end subprogram documentation block + +! module interface: + use m_obsdiagnode, only: obs_diag + use m_obsdiagnode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: cldchnode + + type,extends(obsnode):: cldchnode + !type(cldch_ob_type),pointer :: llpoint => null() + type(obs_diag), pointer :: diags => null() + real(r_kind) :: res ! cldch residual + real(r_kind) :: err2 ! cldch error squared + real(r_kind) :: raterr2 ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b ! variational quality control parameter + real(r_kind) :: pg ! variational quality control parameter + real(r_kind) :: wij(4) ! horizontal interpolation weights + integer(i_kind) :: ij(4) ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + ! procedure:: clean => obsnode_clean_ + end type cldchnode + + public:: cldchnode_typecast + public:: cldchnode_nextcast + interface cldchnode_typecast; module procedure typecast_ ; end interface + interface cldchnode_nextcast; module procedure nextcast_ ; end interface + + public:: cldchnode_appendto + interface cldchnode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_cldchnode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(cldchnode) + use m_obsnode, only: obsnode + implicit none + type(cldchnode),pointer:: ptr_ + class(obsnode ),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(anode) + type is(cldchnode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(cldchnode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type(cldchnode),pointer:: ptr_ + class(obsnode ),target,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) +!-- append anode to linked-list oLL + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(cldchnode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[cldchnode]" +end function mytype + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use m_obsdiagnode, only: obsdiaglookup_locate + implicit none + class(cldchnode), intent(inout):: anode + integer(i_kind), intent(in ):: iunit + integer(i_kind), intent( out):: istat + type(obs_diags), intent(in ):: diaglookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'::obsnode_xread_' + logical:: skip_ +_ENTRY_(myname_) + + skip_=.false. + if(present(skip)) skip_=skip + + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%wij , & + anode%ij + if (istat/=0) then + call perr(myname_,'read(%(res,...)), iostat =',istat) + _EXIT_(myname_) + return + end if + + anode%diags => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,1) + if(.not. associated(anode%diags)) then + call perr(myname_,'obsdiaglookup_locate(), %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call die(myname_) + endif + endif ! if(skip_); else +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(cldchnode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'::obsnode_xwrite_' +_ENTRY_(myname_) + + write(junit,iostat=jstat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%wij , & + anode%ij + if (jstat/=0) then + call perr(myname_,'write(%(res,...)), iostat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + use m_cvgridlookup, only: cvgridlookup_getiw + implicit none + class(cldchnode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_sethop_' +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%ij,anode%wij) +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(cldchnode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_isvalid_' +_ENTRY_(myname_) + isvalid_= associated(anode%diags) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(cldchnode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + anode%diags%tldepart(jiter)*anode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 + return +end subroutine gettlddp_ + +end module m_cldchnode