From 81acc83adae2ad8ecba1db116b3e79dfec6e5def Mon Sep 17 00:00:00 2001 From: Michael Lueken Date: Fri, 19 Jun 2020 17:02:59 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#13. Continuing to clear through coding standard issues in the master. Finished through src/gsi/intlwcp.f90. --- src/gsi/inc2guess.f90 | 48 +- src/gsi/insitu_info.f90 | 67 +-- src/gsi/intall.f90 | 110 ++-- src/gsi/intaod.f90 | 28 +- src/gsi/intcldch.f90 | 28 +- src/gsi/intco.f90 | 208 ++++---- src/gsi/intdbz.f90 | 164 +++--- src/gsi/intdw.f90 | 42 +- src/gsi/intgps.f90 | 82 +-- src/gsi/intgust.f90 | 32 +- src/gsi/inthowv.f90 | 30 +- src/gsi/intjcmod.f90 | 202 ++++---- src/gsi/intjo.f90 | 222 ++++---- src/gsi/intlag.f90 | 54 +- src/gsi/intlcbas.f90 | 28 +- src/gsi/intlight.f90 | 1096 +++++++++++++++++++-------------------- src/gsi/intlwcp.f90 | 162 +++--- 17 files changed, 1304 insertions(+), 1299 deletions(-) diff --git a/src/gsi/inc2guess.f90 b/src/gsi/inc2guess.f90 index 0e5aa0fe9d..d7d0aed317 100644 --- a/src/gsi/inc2guess.f90 +++ b/src/gsi/inc2guess.f90 @@ -59,26 +59,26 @@ subroutine inc2guess(sval) character(len=10),allocatable,dimension(:) :: guess integer(i_kind) i,j,k,it,ii,ic,id,ngases,nguess,ier,istatus real(r_kind) :: zt - real(r_kind),pointer,dimension(:,: ) :: ptr2dinc=>NULL() - real(r_kind),pointer,dimension(:,: ) :: ptr2dges=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ptr3dinc=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ptr3dges=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ges_div_it=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ges_vor_it=>NULL() + real(r_kind),pointer,dimension(:,: ) :: ptr2dinc=>null() + real(r_kind),pointer,dimension(:,: ) :: ptr2dges=>null() + real(r_kind),pointer,dimension(:,:,:) :: ptr3dinc=>null() + real(r_kind),pointer,dimension(:,:,:) :: ptr3dges=>null() + real(r_kind),pointer,dimension(:,:,:) :: ges_div_it=>null() + real(r_kind),pointer,dimension(:,:,:) :: ges_vor_it=>null() ! Inquire about guess fields -call gsi_metguess_get('dim',nguess,istatus) -if (nguess>0) then - allocate(guess(nguess)) - call gsi_metguess_get('gsinames',guess,istatus) -endif + call gsi_metguess_get('dim',nguess,istatus) + if (nguess>0) then + allocate(guess(nguess)) + call gsi_metguess_get('gsinames',guess,istatus) + endif ! Inquire about chemistry fields -call gsi_chemguess_get('dim',ngases,istatus) -if (ngases>0) then - allocate(gases(ngases)) - call gsi_chemguess_get('gsinames',gases,istatus) -endif + call gsi_chemguess_get('dim',ngases,istatus) + if (ngases>0) then + allocate(gases(ngases)) + call gsi_chemguess_get('gsinames',gases,istatus) + endif !******************************************************************************* @@ -86,7 +86,7 @@ subroutine inc2guess(sval) do it=1,nfldsig if (nobs_bins>1) then zt = hrdifsig(it) - ii = NINT(zt/hr_obsbin)+1 + ii = nint(zt/hr_obsbin)+1 else ii = 1 endif @@ -113,7 +113,7 @@ subroutine inc2guess(sval) call gsi_bundlegetpointer (sval(ii), guess(ic),ptr3dinc,istatus) call gsi_bundlegetpointer (gsi_metguess_bundle(it),guess(ic),ptr3dges,istatus) if (trim(guess(ic))=='oz'.or.trim(guess(ic))=='q') then - call copy_positive_fldr3_(ptr3dges,ptr3dinc) + call copy_positive_fldr3_(ptr3dges,ptr3dinc) else ptr3dges = ptr3dinc endif @@ -146,7 +146,7 @@ subroutine inc2guess(sval) end do if(ngases>0)then - deallocate(gases) + deallocate(gases) endif call gsi_bundlegetpointer (sval(ii),'sst',ptr2dinc,istatus) @@ -164,8 +164,9 @@ subroutine inc2guess(sval) return contains subroutine copy_positive_fldr2_(ges,xinc) - real(r_kind),pointer :: ges(:,:) - real(r_kind),pointer :: xinc(:,:) + implicit none + real(r_kind),pointer,intent(inout) :: ges(:,:) + real(r_kind),pointer,intent(in ) :: xinc(:,:) real(r_kind) ana do j=1,lon2 do i=1,lat2 @@ -175,8 +176,9 @@ subroutine copy_positive_fldr2_(ges,xinc) end do end subroutine copy_positive_fldr2_ subroutine copy_positive_fldr3_(ges,xinc) - real(r_kind),pointer :: ges(:,:,:) - real(r_kind),pointer :: xinc(:,:,:) + implicit none + real(r_kind),pointer,intent(inout) :: ges(:,:,:) + real(r_kind),pointer,intent(in ) :: xinc(:,:,:) real(r_kind) ana do k=1,nsig do j=1,lon2 diff --git a/src/gsi/insitu_info.f90 b/src/gsi/insitu_info.f90 index de1f2c6407..ccb9f79313 100644 --- a/src/gsi/insitu_info.f90 +++ b/src/gsi/insitu_info.f90 @@ -50,8 +50,9 @@ module insitu_info subroutine mbuoy_info(mype) !************************************************************************** ! -! assign the depth dependent moored buoy station ID +! assign the depth dependent moored buoy station id ! + implicit none integer(i_kind), intent(in) :: mype allocate(cid_mbuoy(n_3mdiscus)) @@ -60,13 +61,13 @@ subroutine mbuoy_info(mype) ! cid_mbuoy = ' ' ! -! COMPS moored buoy (depth = 1.2m) +! comps moored buoy (depth = 1.2m) ! cid_mbuoy( 1) = '42022' cid_mbuoy( 2) = '42023' cid_mbuoy( 3) = '42024' ! -! SCRIPPS moored buoy (depth = 0.45m) +! scripps moored buoy (depth = 0.45m) ! cid_mbuoy( 4) = '31201' cid_mbuoy( 5) = '41112' @@ -106,7 +107,7 @@ subroutine mbuoy_info(mype) cid_mbuoy(39) = '51203' cid_mbuoy(40) = '52200' ! -! TRITON buoys (depth = 1.5m) +! triton buoys (depth = 1.5m) ! cid_mbuoy(41) = '52071' cid_mbuoy(42) = '52072' @@ -133,7 +134,7 @@ subroutine mbuoy_info(mype) cid_mbuoy(63) = '52045' cid_mbuoy(64) = '52046' ! -! NDBC 3-meter buoy (depth = 0.6m) +! ndbc 3-meter buoy (depth = 0.6m) ! cid_mbuoy(71) = '41004' cid_mbuoy(72) = '41008' @@ -193,7 +194,7 @@ subroutine mbuoy_info(mype) cid_mbuoy(126) = '51001' cid_mbuoy(127) = '51028' ! -! Canadian 3-meter buoy (depth = 0.6m) +! canadian 3-meter buoy (depth = 0.6m) ! cid_mbuoy(128) = '44258' cid_mbuoy(129) = '45132' @@ -219,7 +220,7 @@ subroutine mbuoy_info(mype) cid_mbuoy(149) = '46207' cid_mbuoy(150) = '46208' ! -! MBARI moored buoy (depth = 0.6m) +! mbari moored buoy (depth = 0.6m) ! cid_mbuoy(151) = '46091' cid_mbuoy(152) = '46092' @@ -233,8 +234,9 @@ end subroutine mbuoy_info subroutine mbuoyb_info(mype) !************************************************************************** ! -! assign the depth dependent moored buoyb station ID +! assign the depth dependent moored buoyb station id ! + implicit none integer(i_kind), intent(in) :: mype allocate(cid_mbuoyb(n_3mdiscus)) @@ -243,13 +245,13 @@ subroutine mbuoyb_info(mype) ! cid_mbuoyb = ' ' ! -! COMPS moored buoy (depth = 1.2m) +! comps moored buoy (depth = 1.2m) ! cid_mbuoyb( 1) = '4200022' cid_mbuoyb( 2) = '4200023' cid_mbuoyb( 3) = '4200024' ! -! SCRIPPS moored buoy (depth = 0.45m) +! scripps moored buoy (depth = 0.45m) ! cid_mbuoyb( 4) = '3100201' cid_mbuoyb( 5) = '4100112' @@ -289,7 +291,7 @@ subroutine mbuoyb_info(mype) cid_mbuoyb(39) = '5100203' cid_mbuoyb(40) = '5200200' ! -! TRITON buoys (depth = 1.5m) +! triton buoys (depth = 1.5m) ! cid_mbuoyb(41) = '5200071' cid_mbuoyb(42) = '5200072' @@ -316,7 +318,7 @@ subroutine mbuoyb_info(mype) cid_mbuoyb(63) = '5200045' cid_mbuoyb(64) = '5200046' ! -! NDBC 3-meter buoy (depth = 0.6m) +! ndbc 3-meter buoy (depth = 0.6m) ! cid_mbuoyb(71) = '4100004' cid_mbuoyb(72) = '4100008' @@ -376,7 +378,7 @@ subroutine mbuoyb_info(mype) cid_mbuoyb(126) = '5100001' cid_mbuoyb(127) = '5100028' ! -! Canadian 3-meter buoy (depth = 0.6m) +! canadian 3-meter buoy (depth = 0.6m) ! cid_mbuoyb(128) = '4400258' cid_mbuoyb(129) = '4500132' @@ -402,7 +404,7 @@ subroutine mbuoyb_info(mype) cid_mbuoyb(149) = '4600207' cid_mbuoyb(150) = '4600208' ! -! MBARI moored buoy (depth = 0.6m) +! mbari moored buoy (depth = 0.6m) ! cid_mbuoyb(151) = '4600091' cid_mbuoyb(152) = '4600092' @@ -418,26 +420,27 @@ subroutine read_ship_info(mype) ! ! read ship info from an external file to determine the depth and instrument ! - integer(i_kind), intent(in) :: mype + implicit none + integer(i_kind), intent(in) :: mype - integer(i_kind) ios - logical iexist + integer(i_kind) ios + logical iexist - filename='insituinfo' - inquire(file=trim(filename),exist=iexist) - if(iexist) then - open(lunship,file=filename,form='formatted',iostat=ios) - allocate (ship%id(n_ship),ship%depth(n_ship),ship%sensor(n_ship)) - if(ios==0) then - do i = 1, n_ship - read(lunship,'(a10,f6.1,1x,a5)') ship%id(i),ship%depth(i),ship%sensor(i) - enddo - endif - else - n_ship=0 - allocate (ship%id(n_ship),ship%depth(n_ship),ship%sensor(n_ship)) - endif + filename='insituinfo' + inquire(file=trim(filename),exist=iexist) + if(iexist) then + open(lunship,file=filename,form='formatted',iostat=ios) + allocate (ship%id(n_ship),ship%depth(n_ship),ship%sensor(n_ship)) + if(ios==0) then + do i = 1, n_ship + read(lunship,'(a10,f6.1,1x,a5)') ship%id(i),ship%depth(i),ship%sensor(i) + enddo + endif + else + n_ship=0 + allocate (ship%id(n_ship),ship%depth(n_ship),ship%sensor(n_ship)) + endif - if(mype == 0) write(6,*) ' in read_ship_info, n_ship = ', n_ship + if(mype == 0) write(6,*) ' in read_ship_info, n_ship = ', n_ship end subroutine read_ship_info end module insitu_info diff --git a/src/gsi/intall.f90 b/src/gsi/intall.f90 index 43c2b9f350..0c719a95e2 100644 --- a/src/gsi/intall.f90 +++ b/src/gsi/intall.f90 @@ -14,10 +14,10 @@ module intallmod ! 2009-08-13 lueken - update documentation ! 2012-02-08 kleist - changes related to 4d-ensemble-var additions and consolidation of ! int... individual modules to one intjcmod -! 2015-09-03 guo - obsmod::yobs has been replaced with m_obsHeadBundle, +! 2015-09-03 guo - obsmod::yobs has been replaced with m_obsheadbundle, ! where yobs is created and destroyed when and where it ! is needed. -! 2018-08-10 guo - removed obsHeadBundle references. +! 2018-08-10 guo - removed obsheadbundle references. ! - replaced intjo() related implementations with a new ! polymorphic implementation of intjomod::intjo(). ! @@ -35,8 +35,8 @@ module intallmod implicit none -PRIVATE -PUBLIC intall +private +public intall contains @@ -44,36 +44,36 @@ module intallmod subroutine intall(sval,sbias,rval,rbias) !$$$ subprogram documentation block ! . . . . -! subprogram: intall calculate RHS for analysis equation +! subprogram: intall calculate rhs for analysis equation ! prgmmr: derber org: np23 date: 2003-12-18 ! -! abstract: calculate RHS for all variables (nonlinear qc version) +! abstract: calculate rhs for all variables (nonlinear qc version) ! ! A description of nonlinear qc follows: ! -! The observation penalty Jo is defined as +! The observation penalty jo is defined as ! -! Jo = - (sum over obs) 2*log(Po) +! jo = - (sum over obs) 2*log(Po) ! ! where, ! -! Po = Wnotgross*exp(-.5*(Hn(x+xb) - yo)**2 ) + Wgross +! po = wnotgross*exp(-.5*(hn(x+xb) - yo)**2 ) + wgross ! with -! Hn = the forward model (possibly non-linear) normalized by +! hn = the forward model (possibly non-linear) normalized by ! observation error ! x = the current estimate of the analysis increment ! xb = the background state ! yo = the observation normalized by observation error ! -! Note: The factor 2 in definition of Jo is present because the -! penalty Jo as used in this code is 2*(usual definition +! Note: The factor 2 in definition of jo is present because the +! penalty jo as used in this code is 2*(usual definition ! of penalty) ! -! Wgross = Pgross*cg +! wgross = pgross*cg ! -! Wnotgross = 1 - Wgross +! wnotgross = 1 - wgross ! -! Pgross = probability of gross error for observation (assumed +! pgross = probability of gross error for observation (assumed ! here to have uniform distribution over the possible ! range of values) ! @@ -97,17 +97,17 @@ subroutine intall(sval,sbias,rval,rbias) ! pg_rad=.002 ! probability of gross error for radiances ! ! -! Given the above Jo, the gradient of Jo is as follows: +! Given the above jo, the gradient of jo is as follows: ! -! T -! gradx(Jo) = - (sum over observations) 2*H (Hn(x+xb)-yo)*(Po - Wgross)/Po +! t +! gradx(jo) = - (sum over observations) 2*h (hn(x+xb)-yo)*(po - wgross)/po ! ! where, ! -! H = tangent linear model of Hn about x+xb +! h = tangent linear model of hn about x+xb ! ! -! Note that if Pgross = 0.0, then Wnotgross=1.0 and Wgross=0.0. That is, +! Note that if pgross = 0.0, then wnotgross=1.0 and wgross=0.0. That is, ! the code runs as though nonlinear quality control were not present ! (which is indeed the case since the gross error probability is 0). ! @@ -124,7 +124,7 @@ subroutine intall(sval,sbias,rval,rbias) ! for wind components into int for st,vp ! 2004-11-30 treadon - add brightness temperatures to nonlinear ! quality control -! 2004-12-03 treadon - replace mpe_iallreduce (IBM extension) with +! 2004-12-03 treadon - replace mpe_iallreduce (ibm extension) with ! standard mpi_allreduce ! 2005-01-20 okamoto - add u,v to intrad ! 2005-02-23 wu - changes related to normalized rh option @@ -134,27 +134,27 @@ subroutine intall(sval,sbias,rval,rbias) ! for 2dvar only surface analysis option ! 2005-06-03 parrish - add horizontal derivatives ! 2005-07-10 kleist - add dynamic constraint term -! 2005-09-29 kleist - expand Jc term, include time derivatives vector -! 2005-11-21 kleist - separate tendencies from Jc term, add call to calctends adjoint +! 2005-09-29 kleist - expand jc term, include time derivatives vector +! 2005-11-21 kleist - separate tendencies from jc term, add call to calctends adjoint ! 2005-12-01 cucurull - add code for GPS local bending angle, add use obsmod for ref_obs ! 2005-12-20 parrish - add arguments to call to intt to allow for option of using boundary ! layer forward tlm. ! 2006-02-03 derber - modify to increase reproducibility ! 2006-03-17 park - correct error in call to intt--rval,sval --> rvaluv,svaluv ! in order to correctly pass wind variables. -! 2006-04-06 kleist - include both Jc formulations +! 2006-04-06 kleist - include both jc formulations ! 2006-07-26 parrish - correct inconsistency in computation of space and time derivatives of q ! currently, if derivatives computed, for q it is normalized q, but ! should be mixing ratio. ! 2006-07-26 parrish - add strong constraint initialization option ! 2007-03-19 tremolet - binning of observations -! 2007-04-13 tremolet - split Jo and 3dvar components into intjo and int3dvar +! 2007-04-13 tremolet - split jo and 3dvar components into intjo and int3dvar ! 2007-10-01 todling - add timers ! 2011-10-20 todling - observation operators refer to state- not control-vec (cvars->svars) -! 2014-03-19 pondeca - Add RHS calculation for wspd10m constraint -! 2014-05-07 pondeca - Add RHS calculation for howv constraint -! 2014-06-17 carley/zhu - Add RHS calculation for lcbas constraint -! 2015-07-10 pondeca - Add RHS calculation for cldch constraint +! 2014-03-19 pondeca - Add rhs calculation for wspd10m constraint +! 2014-05-07 pondeca - Add rhs calculation for howv constraint +! 2014-06-17 carley/zhu - Add rhs calculation for lcbas constraint +! 2015-07-10 pondeca - Add rhs calculation for cldch constraint ! 2019-03-13 eliu - add precipitation component ! ! input argument list: @@ -164,7 +164,7 @@ subroutine intall(sval,sbias,rval,rbias) ! rbias ! ! output argument list: -! rval - RHS on grid +! rval - rhs on grid ! rbias ! ! remarks: @@ -219,11 +219,11 @@ subroutine intall(sval,sbias,rval,rbias) qpred=zero_quad -! Compute RHS in physical space (rval,qpred) +! Compute rhs in physical space (rval,qpred) call intjo(rval,qpred,sval,sbias) if(.not.ltlint)then -! RHS for moisture constraint +! rhs for moisture constraint if (.not.ljc4tlevs) then call intlimq(rval(ibin_anl),sval(ibin_anl),ntguessig) else @@ -258,41 +258,41 @@ subroutine intall(sval,sbias,rval,rbias) end do end if end if ! ljclimqc -! RHS for gust constraint +! rhs for gust constraint if (getindex(svars2d,'gust')>0)call intlimg(rval(1),sval(1)) -! RHS for vis constraint +! rhs for vis constraint if (getindex(svars2d,'vis')>0) call intlimv(rval(1),sval(1)) -! RHS for pblh constraint +! rhs for pblh constraint if (getindex(svars2d,'pblh')>0) call intlimp(rval(1),sval(1)) -! RHS for wspd10m constraint +! rhs for wspd10m constraint if (getindex(svars2d,'wspd10m')>0) call intlimw10m(rval(1),sval(1)) -! RHS for howv constraint +! rhs for howv constraint if (getindex(svars2d,'howv')>0) call intlimhowv(rval(1),sval(1)) -! RHS for lcbas constraint +! rhs for lcbas constraint if (getindex(svars2d,'lcbas')>0) call intliml(rval(1),sval(1)) -! RHS for cldch constraint +! rhs for cldch constraint if (getindex(svars2d,'cldch')>0) call intlimcldch(rval(1),sval(1)) end if -! RHS for dry ps constraint: part 1 +! rhs for dry ps constraint: part 1 if(ljcpdry)then - if (.not.ljc4tlevs) then - call intjcpdry1(sval(ibin_anl),1,mass) - else - call intjcpdry1(sval,nobs_bins,mass) - end if + if (.not.ljc4tlevs) then + call intjcpdry1(sval(ibin_anl),1,mass) + else + call intjcpdry1(sval,nobs_bins,mass) + end if -! Put reduces together to minimize wait time -! First, use MPI to get global mean increment - call mpl_allreduce(2*nobs_bins,qpvals=mass) +! Put reduces together to minimize wait time +! First, use MPI to get global mean increment + call mpl_allreduce(2*nobs_bins,qpvals=mass) end if @@ -301,16 +301,16 @@ subroutine intall(sval,sbias,rval,rbias) call mpl_allreduce(nrclen,qpvals=qpred) -! RHS for dry ps constraint: part 2 +! rhs for dry ps constraint: part 2 if(ljcpdry)then - if (.not.ljc4tlevs) then - call intjcpdry2(rval(ibin_anl),1,mass) - else - call intjcpdry2(rval,nobs_bins,mass) - end if + if (.not.ljc4tlevs) then + call intjcpdry2(rval(ibin_anl),1,mass) + else + call intjcpdry2(rval,nobs_bins,mass) + end if end if -! RHS for Jc DFI +! rhs for jc dfi if (ljcdfi .and. nobs_bins>1) call intjcdfi(rval,sval) if(nsclen > 0)then diff --git a/src/gsi/intaod.f90 b/src/gsi/intaod.f90 index c55e59e339..283fb513cb 100644 --- a/src/gsi/intaod.f90 +++ b/src/gsi/intaod.f90 @@ -8,7 +8,7 @@ module intaodmod ! ! program history log: ! 2010-10-20 hclin - modified from intrad for total aod -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intaod_ @@ -21,17 +21,17 @@ module intaodmod ! !$$$ end documentation block - use m_obsNode , only: obsNode - use m_aeroNode , only: aeroNode - use m_aeroNode , only: aeroNode_typecast - use m_aeroNode , only: aeroNode_nextcast + use m_obsnode , only: obsnode + use m_aeronode , only: aeronode + use m_aeronode , only: aeronode_typecast + use m_aeronode , only: aeronode_nextcast implicit none private public :: intaod interface intaod; module procedure & - intaod_ + intaod_ end interface contains @@ -73,11 +73,11 @@ subroutine intaod_(aerohead,rval,sval) use gsi_bundlemod, only: gsi_bundleputvar use gsi_chemguess_mod, only: gsi_chemguess_get use mpeu_util, only: getindex - use m_obsdiagNode, only: obsdiagNode_set + use m_obsdiagnode, only: obsdiagnode_set implicit none ! Declare passed variables - class(obsNode),pointer,intent(in) :: aerohead + class(obsnode),pointer,intent(in) :: aerohead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -88,7 +88,7 @@ subroutine intaod_(aerohead,rval,sval) real(r_kind) val real(r_kind) w1,w2,w3,w4 ! real(r_kind) cg_aero,p0,wnotgross,wgross - type(aeroNode), pointer :: aeroptr + type(aeronode), pointer :: aeroptr real(r_kind),pointer,dimension(:) :: sv_chem real(r_kind),pointer,dimension(:) :: rv_chem @@ -99,7 +99,7 @@ subroutine intaod_(aerohead,rval,sval) if ( naero <= 0 ) return !aeroptr => aerohead - aeroptr => aeroNode_typecast(aerohead) + aeroptr => aeronode_typecast(aerohead) do while (associated(aeroptr)) j1=aeroptr%ij(1) j2=aeroptr%ij(2) @@ -114,7 +114,7 @@ subroutine intaod_(aerohead,rval,sval) tval(k)=zero end do -! Begin Forward model +! Begin forward model ! calculate temperature, q, ozone, sst vector at observation location i1n(1) = j1 i2n(1) = j2 @@ -155,10 +155,10 @@ subroutine intaod_(aerohead,rval,sval) if (lsaveobsens) then val = val*aeroptr%err2(nn)*aeroptr%raterr2(nn) !-- aeroptr%diags(nn)%ptr%obssen(jiter) = val - call obsdiagNode_set(aeroptr%diags(nn)%ptr,jiter=jiter,obssen=val) + call obsdiagnode_set(aeroptr%diags(nn)%ptr,jiter=jiter,obssen=val) else !-- if (aeroptr%luse) aeroptr%diags(nn)%ptr%tldepart(jiter) = val - if (aeroptr%luse) call obsdiagNode_set(aeroptr%diags(nn)%ptr,jiter=jiter,tldepart=val) + if (aeroptr%luse) call obsdiagnode_set(aeroptr%diags(nn)%ptr,jiter=jiter,tldepart=val) endif endif @@ -214,7 +214,7 @@ subroutine intaod_(aerohead,rval,sval) endif ! < l_do_adjoint > !aeroptr => aeroptr%llpoint - aeroptr => aeroNode_nextcast(aeroptr) + aeroptr => aeronode_nextcast(aeroptr) ! call stop2(999) end do diff --git a/src/gsi/intcldch.f90 b/src/gsi/intcldch.f90 index ca75f57368..18a0617b22 100644 --- a/src/gsi/intcldch.f90 +++ b/src/gsi/intcldch.f90 @@ -8,7 +8,7 @@ module intcldchmod ! ! program history log: ! 2015-07-10 Manuel Pondeca -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intcldch @@ -21,15 +21,15 @@ module intcldchmod ! !$$$ end documentation block -use m_obsNode , only: obsNode -use m_cldchNode, only: cldchNode -use m_cldchNode, only: cldchNode_typecast -use m_cldchNode, only: cldchNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode , only: obsnode +use m_cldchnode, only: cldchnode +use m_cldchnode, only: cldchnode_typecast +use m_cldchnode, only: cldchnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intcldch +private +public intcldch contains @@ -70,7 +70,7 @@ subroutine intcldch(cldchhead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in) :: cldchhead + class(obsnode),pointer,intent(in) :: cldchhead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -83,7 +83,7 @@ subroutine intcldch(cldchhead,rval,sval) real(r_kind) cg_cldch,p0,grad,wnotgross,wgross,pg_cldch real(r_kind),pointer,dimension(:) :: scldch real(r_kind),pointer,dimension(:) :: rcldch - type(cldchNode), pointer :: cldchptr + type(cldchnode), pointer :: cldchptr ! Retrieve pointers ! Simply return if any pointer not found @@ -93,7 +93,7 @@ subroutine intcldch(cldchhead,rval,sval) if(ier/=0)return !cldchptr => cldchhead - cldchptr => cldchNode_typecast(cldchhead) + cldchptr => cldchnode_typecast(cldchhead) do while (associated(cldchptr)) j1=cldchptr%ij(1) j2=cldchptr%ij(2) @@ -112,10 +112,10 @@ subroutine intcldch(cldchhead,rval,sval) if (lsaveobsens) then grad = val*cldchptr%raterr2*cldchptr%err2 !-- cldchptr%diags%obssen(jiter) = grad - call obsdiagNode_set(cldchptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(cldchptr%diags,jiter=jiter,obssen=grad) else !-- if (cldchptr%luse) cldchptr%diags%tldepart(jiter)=val - if (cldchptr%luse) call obsdiagNode_set(cldchptr%diags,jiter=jiter,tldepart=val) + if (cldchptr%luse) call obsdiagnode_set(cldchptr%diags,jiter=jiter,tldepart=val) endif endif @@ -148,7 +148,7 @@ subroutine intcldch(cldchhead,rval,sval) endif !cldchptr => cldchptr%llpoint - cldchptr => cldchNode_nextcast(cldchptr) + cldchptr => cldchnode_nextcast(cldchptr) end do diff --git a/src/gsi/intco.f90 b/src/gsi/intco.f90 index d67a361129..6d98dd7e37 100644 --- a/src/gsi/intco.f90 +++ b/src/gsi/intco.f90 @@ -13,8 +13,8 @@ module intcomod ! 2008-11-26 Todling - remove intoz_tl; add interface back ! 2009-08-13 lueken - update documentation ! 2010-06-02 tangborn - converted intoz into intco -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - implemented obs adjoint test +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intco_ @@ -28,18 +28,18 @@ module intcomod ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_colvkNode , only: colvkNode -use m_colvkNode , only: colvkNode_typecast -use m_colvkNode , only: colvkNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_colvknode , only: colvknode +use m_colvknode , only: colvknode_typecast +use m_colvknode , only: colvknode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intco +private +public intco interface intco; module procedure & - intco_ + intco_ end interface contains @@ -61,7 +61,7 @@ subroutine intco_(colvkhead,rval,sval) ! 2010-06-02 tangborn - made version for carbon monoxide ! ! input argument list: -! colvkhead - level carbon monoxide obs type pointer to obs structure for MOPITT +! colvkhead - level carbon monoxide obs type pointer to obs structure for mopitt ! sco - carbon monoxide increment in grid space ! ! output argument list: @@ -77,7 +77,7 @@ subroutine intco_(colvkhead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in ) :: colvkhead + class(obsnode),pointer,intent(in ) :: colvkhead type(gsi_bundle),intent(in ) :: sval type(gsi_bundle),intent(inout) :: rval @@ -98,7 +98,7 @@ subroutine intcolev_(colvkhead,rval,sval) ! program history log: ! 1995-07-11 derber ! 2010-06-07 tangborn - carbon monoxide based on ozone code -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - introduced ladtest_obs ! 2014-12-03 derber - modify so that use of obsdiags can be turned off ! ! input argument list: @@ -125,7 +125,7 @@ subroutine intcolev_(colvkhead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in ) :: colvkhead + class(obsnode),pointer,intent(in ) :: colvkhead type(gsi_bundle) ,intent(in ) :: sval type(gsi_bundle) ,intent(inout) :: rval @@ -142,7 +142,7 @@ subroutine intcolev_(colvkhead,rval,sval) real(r_kind),allocatable,dimension(:) :: coak real(r_kind),allocatable,dimension(:) :: vali real(r_kind),allocatable,dimension(:) :: val_ret - type(colvkNode), pointer :: colvkptr + type(colvknode), pointer :: colvkptr ! If no co observations return if(.not. associated(colvkhead))return @@ -166,11 +166,11 @@ subroutine intcolev_(colvkhead,rval,sval) enddo enddo ! -! MOPITT CARBON MONOXIDE: LAYER CO +! mopitt carbon monoxide: layer co ! ! Loop over carbon monoxide observations. !colvkptr => colvkhead - colvkptr => colvkNode_typecast(colvkhead) + colvkptr => colvknode_typecast(colvkhead) do while (associated(colvkptr)) ! Set location @@ -185,107 +185,107 @@ subroutine intcolev_(colvkhead,rval,sval) ! with several of the terms already calculated - allocate(vali(colvkptr%nlco)) - allocate(coak(colvkptr%nlco)) - allocate(val_ret(colvkptr%nlco)) - - do k=1,colvkptr%nlco ! loop over MOPITT ave. ker. contribution levels - pob = colvkptr%prs(k) - k1=int(pob) - k2=min(k1+1,nsig) - w1=colvkptr%wij(1,k) - w2=colvkptr%wij(2,k) - w3=colvkptr%wij(3,k) - w4=colvkptr%wij(4,k) - w5=colvkptr%wij(5,k) - w6=colvkptr%wij(6,k) - w7=colvkptr%wij(7,k) - w8=colvkptr%wij(8,k) - val1= w1* sco(j1,k1)+ & - w2* sco(j2,k1)+ & - w3* sco(j3,k1)+ & - w4* sco(j4,k1)+ & - w5* sco(j1,k2)+ & - w6* sco(j2,k2)+ & - w7* sco(j3,k2)+ & - w8* sco(j4,k2) - vali(k)=val1 - enddo - -! Averaging kernel + allocate(vali(colvkptr%nlco)) + allocate(coak(colvkptr%nlco)) + allocate(val_ret(colvkptr%nlco)) + + do k=1,colvkptr%nlco ! loop over mopitt ave. ker. contribution levels + pob = colvkptr%prs(k) + k1=int(pob) + k2=min(k1+1,nsig) + w1=colvkptr%wij(1,k) + w2=colvkptr%wij(2,k) + w3=colvkptr%wij(3,k) + w4=colvkptr%wij(4,k) + w5=colvkptr%wij(5,k) + w6=colvkptr%wij(6,k) + w7=colvkptr%wij(7,k) + w8=colvkptr%wij(8,k) + val1= w1* sco(j1,k1)+ & + w2* sco(j2,k1)+ & + w3* sco(j3,k1)+ & + w4* sco(j4,k1)+ & + w5* sco(j1,k2)+ & + w6* sco(j2,k2)+ & + w7* sco(j3,k2)+ & + w8* sco(j4,k2) + vali(k)=val1 + enddo - do k=1,colvkptr%nlco ! loop over MOPITT retrieval levels - val1=zero_quad - do j=1,colvkptr%nlco ! loop over MOPITT ak levels - val1=val1+colvkptr%ak(k,j)*vali(j) - enddo +! Averaging kernel + + do k=1,colvkptr%nlco ! loop over mopitt retrieval levels + val1=zero_quad + do j=1,colvkptr%nlco ! loop over mopitt ak levels + val1=val1+colvkptr%ak(k,j)*vali(j) + enddo + + if(luse_obsdiag)then + if (lsaveobsens) then + valx=val1*colvkptr%err2(k)*colvkptr%raterr2(k) + !-- colvkptr%diags(k)%ptr%obssen(jiter)=valx + call obsdiagnode_set(colvkptr%diags(k)%ptr,jiter=jiter,obssen=real(valx,r_kind)) + else + !-- if (colvkptr%luse) colvkptr%diags(k)%ptr%tldepart(jiter)=val1 + if (colvkptr%luse) call obsdiagnode_set(colvkptr%diags(k)%ptr,tldepart=real(val1,r_kind)) + endif + endif - if(luse_obsdiag)then - if (lsaveobsens) then - valx=val1*colvkptr%err2(k)*colvkptr%raterr2(k) - !-- colvkptr%diags(k)%ptr%obssen(jiter)=valx - call obsdiagNode_set(colvkptr%diags(k)%ptr,jiter=jiter,obssen=real(valx,r_kind)) + if (l_do_adjoint) then + if (.not. lsaveobsens) then + if( ladtest_obs ) then + valx = val1 else - !-- if (colvkptr%luse) colvkptr%diags(k)%ptr%tldepart(jiter)=val1 - if (colvkptr%luse) call obsdiagNode_set(colvkptr%diags(k)%ptr,tldepart=real(val1,r_kind)) - endif - endif + val1=val1-colvkptr%res(k) - if (l_do_adjoint) then - if (.not. lsaveobsens) then - if( ladtest_obs ) then - valx = val1 - else - val1=val1-colvkptr%res(k) - - valx=val1*colvkptr%err2(k) - valx=valx*colvkptr%raterr2(k) - end if - endif - val_ret(k)=valx - endif - enddo ! k + valx=val1*colvkptr%err2(k) + valx=valx*colvkptr%raterr2(k) + end if + endif + val_ret(k)=valx + endif + enddo ! k ! Averaging kernel First - spread values to ak contribution levels - if(l_do_adjoint)then - do k=1,colvkptr%nlco !loop over ak levels - coak(k)=zero_quad - do j=1,colvkptr%nlco !loop over profile levels - coak(k)=coak(k)+colvkptr%ak(j,k)*val_ret(j) ! Contribution to kth ak level from jth retrieval level - enddo - enddo + if(l_do_adjoint)then + do k=1,colvkptr%nlco !loop over ak levels + coak(k)=zero_quad + do j=1,colvkptr%nlco !loop over profile levels + coak(k)=coak(k)+colvkptr%ak(j,k)*val_ret(j) ! Contribution to kth ak level from jth retrieval level + enddo + enddo ! Adjoint of interpolation - spreads each ave. kernel level to interpolant gridpoints - do kk=colvkptr%nlco,1,-1 !loop over averaging kernel levels - pob = colvkptr%prs(kk) - k1=int(pob) - k2=min(k1+1,nsig) - w1=colvkptr%wij(1,kk) - w2=colvkptr%wij(2,kk) - w3=colvkptr%wij(3,kk) - w4=colvkptr%wij(4,kk) - w5=colvkptr%wij(5,kk) - w6=colvkptr%wij(6,kk) - w7=colvkptr%wij(7,kk) - w8=colvkptr%wij(8,kk) - rco(j1,k1) = rco(j1,k1) + coak(kk)*w1 - rco(j2,k1) = rco(j2,k1) + coak(kk)*w2 - rco(j3,k1) = rco(j3,k1) + coak(kk)*w3 - rco(j4,k1) = rco(j4,k1) + coak(kk)*w4 - rco(j1,k2) = rco(j1,k2) + coak(kk)*w5 - rco(j2,k2) = rco(j2,k2) + coak(kk)*w6 - rco(j3,k2) = rco(j3,k2) + coak(kk)*w7 - rco(j4,k2) = rco(j4,k2) + coak(kk)*w8 - enddo ! kk + do kk=colvkptr%nlco,1,-1 !loop over averaging kernel levels + pob = colvkptr%prs(kk) + k1=int(pob) + k2=min(k1+1,nsig) + w1=colvkptr%wij(1,kk) + w2=colvkptr%wij(2,kk) + w3=colvkptr%wij(3,kk) + w4=colvkptr%wij(4,kk) + w5=colvkptr%wij(5,kk) + w6=colvkptr%wij(6,kk) + w7=colvkptr%wij(7,kk) + w8=colvkptr%wij(8,kk) + rco(j1,k1) = rco(j1,k1) + coak(kk)*w1 + rco(j2,k1) = rco(j2,k1) + coak(kk)*w2 + rco(j3,k1) = rco(j3,k1) + coak(kk)*w3 + rco(j4,k1) = rco(j4,k1) + coak(kk)*w4 + rco(j1,k2) = rco(j1,k2) + coak(kk)*w5 + rco(j2,k2) = rco(j2,k2) + coak(kk)*w6 + rco(j3,k2) = rco(j3,k2) + coak(kk)*w7 + rco(j4,k2) = rco(j4,k2) + coak(kk)*w8 + enddo ! kk deallocate(coak,vali,val_ret) - endif ! l_do_adjoint - !colvkptr => colvkptr%llpoint - colvkptr => colvkNode_nextcast(colvkptr) + endif ! l_do_adjoint + !colvkptr => colvkptr%llpoint + colvkptr => colvknode_nextcast(colvkptr) ! End loop over observations enddo diff --git a/src/gsi/intdbz.f90 b/src/gsi/intdbz.f90 index aa225aa83f..6b632d781e 100644 --- a/src/gsi/intdbz.f90 +++ b/src/gsi/intdbz.f90 @@ -8,7 +8,7 @@ module intdbzmod ! ! program history log: ! 2017-05-12 Y. Wang and X. Wang - add tangent linear of dbz operator to directly assimilate reflectivity -! for both ARW and NMMB models (Wang and Wang 2017 MWR). POC: xuguang.wang@ou.edu +! for both arw and nmmb models (wang and wang 2017 mwr). ! 2019-07-11 todling - introduced wrf_vars_mod ! ! subroutines included: @@ -22,18 +22,18 @@ module intdbzmod ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_dbzNode, only: dbzNode -use m_dbzNode, only: dbzNode_typecast -use m_dbzNode, only: dbzNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_dbznode, only: dbznode +use m_dbznode, only: dbznode_typecast +use m_dbznode, only: dbznode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intdbz +private +public intdbz interface intdbz; module procedure & - intdbz_ + intdbz_ end interface contains @@ -57,15 +57,15 @@ subroutine intdbz_(dbzhead,rval,sval) ! 2005-08-02 derber - modify for variational qc parameters for each ob ! 2005-09-28 derber - consolidate location and weight arrays ! 2006-07-28 derber - modify to use new inner loop obs data structure -! - unify NL qc +! - unify nl qc ! 2007-02-15 rancic - add foto ! 2007-03-19 tremolet - binning of observations ! 2007-06-05 tremolet - use observation diagnostics structure ! 2007-07-09 tremolet - observation sensitivity -! 2008-01-04 tremolet - Don't apply H^T if l_do_adjoint is false -! 2008-11-28 todling - turn FOTO optional; changed ptr%time handle +! 2008-01-04 tremolet - Don't apply h^t if l_do_adjoint is false +! 2008-11-28 todling - turn foto optional; changed ptr%time handle ! 2010-05-13 todlng - update to use gsi_bundle; update interface -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - introduced ladtest_obs ! 2014-12-03 derber - modify so that use of obsdiags can be turned off ! ! attributes: @@ -87,7 +87,7 @@ subroutine intdbz_(dbzhead,rval,sval) implicit none ! Declare passed variables - class(obsNode), pointer, intent(in ) :: dbzhead + class(obsnode), pointer, intent(in ) :: dbzhead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -99,35 +99,35 @@ subroutine intdbz_(dbzhead,rval,sval) real(r_kind) qrtl,qstl, qgtl real(r_kind),pointer,dimension(:) :: sqr,sqs,sqg,sdbz real(r_kind),pointer,dimension(:) :: rqr,rqs,rqg,rdbz - type(dbzNode), pointer :: dbzptr + type(dbznode), pointer :: dbzptr -! If no dbz data return +! If no dbz data return if(.not. associated(dbzhead))return ! Retrieve pointers ! Simply return if any pointer not found ier=0 if(dbz_exist)then - call gsi_bundlegetpointer(sval,'dbz',sdbz,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'dbz',rdbz,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'dbz',sdbz,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'dbz',rdbz,istatus);ier=istatus+ier else - call gsi_bundlegetpointer(sval,'qr',sqr,istatus);ier=istatus+ier - if (wrf_mass_regional) then - call gsi_bundlegetpointer(sval,'qs',sqs,istatus);ier=istatus+ier - call gsi_bundlegetpointer(sval,'qg',sqg,istatus);ier=istatus+ier - end if + call gsi_bundlegetpointer(sval,'qr',sqr,istatus);ier=istatus+ier + if (wrf_mass_regional) then + call gsi_bundlegetpointer(sval,'qs',sqs,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'qg',sqg,istatus);ier=istatus+ier + end if - call gsi_bundlegetpointer(rval,'qr',rqr,istatus);ier=istatus+ier - if (wrf_mass_regional) then - call gsi_bundlegetpointer(rval,'qs',rqs,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'qg',rqg,istatus);ier=istatus+ier - end if + call gsi_bundlegetpointer(rval,'qr',rqr,istatus);ier=istatus+ier + if (wrf_mass_regional) then + call gsi_bundlegetpointer(rval,'qs',rqs,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'qg',rqg,istatus);ier=istatus+ier + end if end if if(ier/=0)return - dbzptr => dbzNode_typecast(dbzhead) + dbzptr => dbznode_typecast(dbzhead) do while (associated(dbzptr)) j1=dbzptr%ij(1) j2=dbzptr%ij(2) @@ -149,20 +149,20 @@ subroutine intdbz_(dbzhead,rval,sval) ! Forward mode l if( dbz_exist )then - val = w1* sdbz(j1)+w2* sdbz(j2)+w3* sdbz(j3)+w4* sdbz(j4)+ & - w5* sdbz(j5)+w6* sdbz(j6)+w7* sdbz(j7)+w8* sdbz(j8) + val = w1* sdbz(j1)+w2* sdbz(j2)+w3* sdbz(j3)+w4* sdbz(j4)+ & + w5* sdbz(j5)+w6* sdbz(j6)+w7* sdbz(j7)+w8* sdbz(j8) else - qrtl = w1* sqr(j1)+w2* sqr(j2)+w3* sqr(j3)+w4* sqr(j4)+ & - w5* sqr(j5)+w6* sqr(j6)+w7* sqr(j7)+w8* sqr(j8) - if ( wrf_mass_regional )then - qstl = w1* sqs(j1)+w2* sqs(j2)+w3* sqs(j3)+w4* sqs(j4)+ & - w5* sqs(j5)+w6* sqs(j6)+w7* sqs(j7)+w8* sqs(j8) + qrtl = w1* sqr(j1)+w2* sqr(j2)+w3* sqr(j3)+w4* sqr(j4)+ & + w5* sqr(j5)+w6* sqr(j6)+w7* sqr(j7)+w8* sqr(j8) + if ( wrf_mass_regional )then + qstl = w1* sqs(j1)+w2* sqs(j2)+w3* sqs(j3)+w4* sqs(j4)+ & + w5* sqs(j5)+w6* sqs(j6)+w7* sqs(j7)+w8* sqs(j8) - qgtl = w1* sqg(j1)+w2* sqg(j2)+w3* sqg(j3)+w4* sqg(j4)+ & - w5* sqg(j5)+w6* sqg(j6)+w7* sqg(j7)+w8* sqg(j8) + qgtl = w1* sqg(j1)+w2* sqg(j2)+w3* sqg(j3)+w4* sqg(j4)+ & + w5* sqg(j5)+w6* sqg(j6)+w7* sqg(j7)+w8* sqg(j8) - val = (dbzptr%jqr)*qrtl + (dbzptr%jqs)*qstl + (dbzptr%jqg)*qgtl - end if + val = (dbzptr%jqr)*qrtl + (dbzptr%jqs)*qstl + (dbzptr%jqg)*qgtl + end if end if @@ -170,11 +170,11 @@ subroutine intdbz_(dbzhead,rval,sval) if (lsaveobsens) then grad = val*dbzptr%raterr2*dbzptr%err2 !-- dbzptr%diags%obssen(jiter) = grad - call obsdiagNode_set(dbzptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(dbzptr%diags,jiter=jiter,obssen=grad) else !-- if (dbzptr%luse) dbzptr%diags%tldepart(jiter)=val - if (dbzptr%luse) call obsdiagNode_set(dbzptr%diags,jiter=jiter,tldepart=val) + if (dbzptr%luse) call obsdiagnode_set(dbzptr%diags,jiter=jiter,tldepart=val) endif endif @@ -203,53 +203,53 @@ subroutine intdbz_(dbzhead,rval,sval) ! Adjoint if(dbz_exist)then - valdbz = grad - rdbz(j1)=rdbz(j1)+w1*valdbz - rdbz(j2)=rdbz(j2)+w2*valdbz - rdbz(j3)=rdbz(j3)+w3*valdbz - rdbz(j4)=rdbz(j4)+w4*valdbz - rdbz(j5)=rdbz(j5)+w5*valdbz - rdbz(j6)=rdbz(j6)+w6*valdbz - rdbz(j7)=rdbz(j7)+w7*valdbz - rdbz(j8)=rdbz(j8)+w8*valdbz + valdbz = grad + rdbz(j1)=rdbz(j1)+w1*valdbz + rdbz(j2)=rdbz(j2)+w2*valdbz + rdbz(j3)=rdbz(j3)+w3*valdbz + rdbz(j4)=rdbz(j4)+w4*valdbz + rdbz(j5)=rdbz(j5)+w5*valdbz + rdbz(j6)=rdbz(j6)+w6*valdbz + rdbz(j7)=rdbz(j7)+w7*valdbz + rdbz(j8)=rdbz(j8)+w8*valdbz else - valqr = dbzptr%jqr*grad - rqr(j1)=rqr(j1)+w1*valqr - rqr(j2)=rqr(j2)+w2*valqr - rqr(j3)=rqr(j3)+w3*valqr - rqr(j4)=rqr(j4)+w4*valqr - rqr(j5)=rqr(j5)+w5*valqr - rqr(j6)=rqr(j6)+w6*valqr - rqr(j7)=rqr(j7)+w7*valqr - rqr(j8)=rqr(j8)+w8*valqr - if ( wrf_mass_regional )then - valqs=dbzptr%jqs*grad - valqg=dbzptr%jqg*grad - - rqs(j1)=rqs(j1)+w1*valqs - rqs(j2)=rqs(j2)+w2*valqs - rqs(j3)=rqs(j3)+w3*valqs - rqs(j4)=rqs(j4)+w4*valqs - rqs(j5)=rqs(j5)+w5*valqs - rqs(j6)=rqs(j6)+w6*valqs - rqs(j7)=rqs(j7)+w7*valqs - rqs(j8)=rqs(j8)+w8*valqs + valqr = dbzptr%jqr*grad + rqr(j1)=rqr(j1)+w1*valqr + rqr(j2)=rqr(j2)+w2*valqr + rqr(j3)=rqr(j3)+w3*valqr + rqr(j4)=rqr(j4)+w4*valqr + rqr(j5)=rqr(j5)+w5*valqr + rqr(j6)=rqr(j6)+w6*valqr + rqr(j7)=rqr(j7)+w7*valqr + rqr(j8)=rqr(j8)+w8*valqr + if ( wrf_mass_regional )then + valqs=dbzptr%jqs*grad + valqg=dbzptr%jqg*grad + + rqs(j1)=rqs(j1)+w1*valqs + rqs(j2)=rqs(j2)+w2*valqs + rqs(j3)=rqs(j3)+w3*valqs + rqs(j4)=rqs(j4)+w4*valqs + rqs(j5)=rqs(j5)+w5*valqs + rqs(j6)=rqs(j6)+w6*valqs + rqs(j7)=rqs(j7)+w7*valqs + rqs(j8)=rqs(j8)+w8*valqs - rqg(j1)=rqg(j1)+w1*valqg - rqg(j2)=rqg(j2)+w2*valqg - rqg(j3)=rqg(j3)+w3*valqg - rqg(j4)=rqg(j4)+w4*valqg - rqg(j5)=rqg(j5)+w5*valqg - rqg(j6)=rqg(j6)+w6*valqg - rqg(j7)=rqg(j7)+w7*valqg - rqg(j8)=rqg(j8)+w8*valqg - end if + rqg(j1)=rqg(j1)+w1*valqg + rqg(j2)=rqg(j2)+w2*valqg + rqg(j3)=rqg(j3)+w3*valqg + rqg(j4)=rqg(j4)+w4*valqg + rqg(j5)=rqg(j5)+w5*valqg + rqg(j6)=rqg(j6)+w6*valqg + rqg(j7)=rqg(j7)+w7*valqg + rqg(j8)=rqg(j8)+w8*valqg + end if end if endif !dbzptr => dbzptr%llpoint - dbzptr => dbzNode_nextcast(dbzptr) + dbzptr => dbznode_nextcast(dbzptr) end do return end subroutine intdbz_ diff --git a/src/gsi/intdw.f90 b/src/gsi/intdw.f90 index ba2f973dc6..c3a0839906 100644 --- a/src/gsi/intdw.f90 +++ b/src/gsi/intdw.f90 @@ -12,8 +12,8 @@ module intdwmod ! 2005-11-16 Derber - remove interfaces ! 2008-11-26 Todling - remove intdw_tl ! 2009-08-13 lueken - updated documentation -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - implemented obs adjoint test +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intdw_ @@ -26,18 +26,18 @@ module intdwmod ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_dwNode , only: dwNode -use m_dwNode , only: dwNode_typecast -use m_dwNode , only: dwNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_dwnode , only: dwnode +use m_dwnode , only: dwnode_typecast +use m_dwnode , only: dwnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intdw +private +public intdw interface intdw; module procedure & - intdw_ + intdw_ end interface contains @@ -61,15 +61,15 @@ subroutine intdw_(dwhead,rval,sval) ! 2005-08-02 derber - modify for variational qc parameters for each ob ! 2005-09-28 derber - consolidate location and weight arrays ! 2006-07-28 derber - modify to use new inner loop obs data structure -! - unify NL qc +! - unify nl qc ! 2007-03-19 tremolet - binning of observations ! 2007-06-05 tremolet - use observation diagnostics structure ! 2007-07-09 tremolet - observation sensitivity ! 2008-06-02 safford - rm unused vars -! 2008-01-04 tremolet - Don't apply H^T if l_do_adjoint is false -! 2008-11-28 todling - turn FOTO optional; changed ptr%time handle +! 2008-01-04 tremolet - Don't apply h^t if l_do_adjoint is false +! 2008-11-28 todling - turn foto optional; changed ptr%time handle ! 2010-05-13 todling - update to use gsi_bundle; update interface -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - introduced ladtest_obs ! 2014-12-03 derber - modify so that use of obsdiags can be turned off ! ! usage: call intdw(ru,rv,su,sv) @@ -99,7 +99,7 @@ subroutine intdw_(dwhead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in ) :: dwhead + class(obsnode),pointer,intent(in ) :: dwhead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -110,7 +110,7 @@ subroutine intdw_(dwhead,rval,sval) real(r_kind) cg_dw,p0,grad,wnotgross,wgross real(r_kind),pointer,dimension(:) :: su,sv real(r_kind),pointer,dimension(:) :: ru,rv - type(dwNode), pointer :: dwptr + type(dwnode), pointer :: dwptr ! If no dw observations return if(.not. associated(dwhead))return @@ -124,7 +124,7 @@ subroutine intdw_(dwhead,rval,sval) if(ier/=0)return !dwptr => dwhead - dwptr => dwNode_typecast(dwhead) + dwptr => dwnode_typecast(dwhead) do while (associated(dwptr)) j1=dwptr%ij(1) j2=dwptr%ij(2) @@ -154,14 +154,14 @@ subroutine intdw_(dwhead,rval,sval) if (lsaveobsens) then grad = val * dwptr%raterr2 * dwptr%err2 !-- dwptr%diags%obssen(jiter) = grad - call obsdiagNode_set(dwptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(dwptr%diags,jiter=jiter,obssen=grad) else !-- if (dwptr%luse) dwptr%diags%tldepart(jiter)=val - if (dwptr%luse) call obsdiagNode_set(dwptr%diags,jiter=jiter,tldepart=val) + if (dwptr%luse) call obsdiagnode_set(dwptr%diags,jiter=jiter,tldepart=val) endif endif -! Do Adjoint +! Do adjoint if (l_do_adjoint) then if (.not. lsaveobsens) then if( .not. ladtest_obs) val=val-dwptr%res @@ -206,7 +206,7 @@ subroutine intdw_(dwhead,rval,sval) endif !dwptr => dwptr%llpoint - dwptr => dwNode_nextcast(dwptr) + dwptr => dwnode_nextcast(dwptr) end do diff --git a/src/gsi/intgps.f90 b/src/gsi/intgps.f90 index bc78db085e..3f7f6a4c4b 100644 --- a/src/gsi/intgps.f90 +++ b/src/gsi/intgps.f90 @@ -13,7 +13,7 @@ module intgpsmod ! 2008-11-28 Todling - add interface back ! 2009-08-13 lueken - update documentation ! 2013-10-28 todling - rename p3d to prse -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intgps_ @@ -25,18 +25,18 @@ module intgpsmod ! machine: ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_gpsNode, only: gpsNode -use m_gpsNode, only: gpsNode_typecast -use m_gpsNode, only: gpsNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_gpsnode, only: gpsnode +use m_gpsnode, only: gpsnode_typecast +use m_gpsnode, only: gpsnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intgps +private +public intgps interface intgps; module procedure & - intgps_ + intgps_ end interface contains @@ -57,8 +57,8 @@ subroutine intgps_(gpshead,rval,sval) ! 2004-10-08 parrish - add nonlinear qc option ! 2004-11-19 cucurull- add increments for surface pressure and temperature at levels ! below observation. Install non-linear forward operator. -! 2005-01-26 cucurull- Implement local GPS RO operator -! 2005-03-01 parrish - nonlinear qc change as above; correct bug in zeroing of tl_AD +! 2005-01-26 cucurull- Implement local gps ro operator +! 2005-03-01 parrish - nonlinear qc change as above; correct bug in zeroing of tl_ad ! 2005-03-23 cucurull- correct bounds for obs below the second level; place ! bounds for k1 and k2 ! 2005-04-11 treadon - merge intref and intref_qc into single routine @@ -67,15 +67,15 @@ subroutine intgps_(gpshead,rval,sval) ! 2005-12-02 cucurull - fix bug for dimensions of sp and rp ! 2006-01-03 treadon - include r_kind type in w1,w2,...,w12 declaration ! 2006-07-28 derber - modify to use new inner loop obs data structure -! - unify NL qc +! - unify nl qc ! 2006-09-06 cucurull - generalize code to hybrid vertical coordinate and modify to use ! surface pressure ! 2007-01-13 derber - clean up code and use coding standards ! 2007-03-28 derber - turn intref into generalized intgps ! 2007-07-26 cucurull - in/out 3d pressure to update code to generalized vertical coordinate ! 2008-06-02 safford - rm unused vars -! 2008-11-26 todling - add 4dvar and GSI adjoint capability (obs binning, obsens, etc) -! 2008-11-26 todling - turned FOTO optional; changed ptr%time handle +! 2008-11-26 todling - add 4dvar and gsi adjoint capability (obs binning, obsens, etc) +! 2008-11-26 todling - turned foto optional; changed ptr%time handle ! 2010-05-13 todling - update to use gsi_bundle; update interface ! 2014-12-03 derber - modify so that use of obsdiags can be turned off ! @@ -83,7 +83,7 @@ subroutine intgps_(gpshead,rval,sval) ! gpshead - obs type pointer to obs structure ! st - input temperature correction field ! sq - input q correction field -! sp - input (3D) p correction field +! sp - input (3d) p correction field ! ! output argument list: ! gpshead - obs type pointer to obs structure @@ -108,7 +108,7 @@ subroutine intgps_(gpshead,rval,sval) implicit none ! Declare passed variables - class(obsNode), pointer, intent(in ) :: gpshead + class(obsnode), pointer, intent(in ) :: gpshead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -116,14 +116,14 @@ subroutine intgps_(gpshead,rval,sval) integer(i_kind) j,ier,istatus integer(i_kind),dimension(nsig):: i1,i2,i3,i4 real(r_kind) :: w1,w2,w3,w4 - real(r_kind) :: p_TL,p_AD,t_TL,t_AD,q_TL,q_AD + real(r_kind) :: p_tl,p_ad,t_tl,t_ad,q_tl,q_ad real(r_kind) :: val,pg_gps real(r_kind) ::cg_gps,grad,p0,wnotgross,wgross real(r_kind),pointer,dimension(:) :: st,sq real(r_kind),pointer,dimension(:) :: rt,rq real(r_kind),pointer,dimension(:) :: sp real(r_kind),pointer,dimension(:) :: rp - type(gpsNode), pointer :: gpsptr + type(gpsnode), pointer :: gpsptr ! If no gps obs return if(.not. associated(gpshead))return @@ -139,7 +139,7 @@ subroutine intgps_(gpshead,rval,sval) if(ier/=0)return !gpsptr => gpshead - gpsptr => gpsNode_typecast(gpshead) + gpsptr => gpsnode_typecast(gpshead) do while (associated(gpsptr)) ! Load location information into local variables @@ -160,20 +160,20 @@ subroutine intgps_(gpshead,rval,sval) ! local refractivity (linear operator) do j=1,nsig - t_TL=w1* st(i1(j))+w2* st(i2(j))+w3* st(i3(j))+w4* st(i4(j)) - q_TL=w1* sq(i1(j))+w2* sq(i2(j))+w3* sq(i3(j))+w4* sq(i4(j)) - p_TL=w1* sp(i1(j))+w2* sp(i2(j))+w3* sp(i3(j))+w4* sp(i4(j)) - val = val + p_TL*gpsptr%jac_p(j) + t_TL*gpsptr%jac_t(j)+q_TL*gpsptr%jac_q(j) + t_tl=w1* st(i1(j))+w2* st(i2(j))+w3* st(i3(j))+w4* st(i4(j)) + q_tl=w1* sq(i1(j))+w2* sq(i2(j))+w3* sq(i3(j))+w4* sq(i4(j)) + p_tl=w1* sp(i1(j))+w2* sp(i2(j))+w3* sp(i3(j))+w4* sp(i4(j)) + val = val + p_tl*gpsptr%jac_p(j) + t_tl*gpsptr%jac_t(j)+q_tl*gpsptr%jac_q(j) end do if (luse_obsdiag)then if (lsaveobsens) then grad = val*gpsptr%raterr2*gpsptr%err2 !-- gpsptr%diags%obssen(jiter) = grad - call obsdiagNode_set(gpsptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(gpsptr%diags,jiter=jiter,obssen=grad) else !-- if (gpsptr%luse) gpsptr%diags%tldepart(jiter)=val - if (gpsptr%luse) call obsdiagNode_set(gpsptr%diags,jiter=jiter,tldepart=val) + if (gpsptr%luse) call obsdiagnode_set(gpsptr%diags,jiter=jiter,tldepart=val) endif endif @@ -205,27 +205,27 @@ subroutine intgps_(gpshead,rval,sval) ! adjoint do j=1,nsig - t_AD = grad*gpsptr%jac_t(j) - rt(i1(j))=rt(i1(j))+w1*t_AD - rt(i2(j))=rt(i2(j))+w2*t_AD - rt(i3(j))=rt(i3(j))+w3*t_AD - rt(i4(j))=rt(i4(j))+w4*t_AD - q_AD = grad*gpsptr%jac_q(j) - rq(i1(j))=rq(i1(j))+w1*q_AD - rq(i2(j))=rq(i2(j))+w2*q_AD - rq(i3(j))=rq(i3(j))+w3*q_AD - rq(i4(j))=rq(i4(j))+w4*q_AD - p_AD = grad*gpsptr%jac_p(j) - rp(i1(j))=rp(i1(j))+w1*p_AD - rp(i2(j))=rp(i2(j))+w2*p_AD - rp(i3(j))=rp(i3(j))+w3*p_AD - rp(i4(j))=rp(i4(j))+w4*p_AD + t_ad = grad*gpsptr%jac_t(j) + rt(i1(j))=rt(i1(j))+w1*t_ad + rt(i2(j))=rt(i2(j))+w2*t_ad + rt(i3(j))=rt(i3(j))+w3*t_ad + rt(i4(j))=rt(i4(j))+w4*t_ad + q_ad = grad*gpsptr%jac_q(j) + rq(i1(j))=rq(i1(j))+w1*q_ad + rq(i2(j))=rq(i2(j))+w2*q_ad + rq(i3(j))=rq(i3(j))+w3*q_ad + rq(i4(j))=rq(i4(j))+w4*q_ad + p_ad = grad*gpsptr%jac_p(j) + rp(i1(j))=rp(i1(j))+w1*p_ad + rp(i2(j))=rp(i2(j))+w2*p_ad + rp(i3(j))=rp(i3(j))+w3*p_ad + rp(i4(j))=rp(i4(j))+w4*p_ad enddo endif !gpsptr => gpsptr%llpoint - gpsptr => gpsNode_nextcast(gpsptr) + gpsptr => gpsnode_nextcast(gpsptr) end do diff --git a/src/gsi/intgust.f90 b/src/gsi/intgust.f90 index 6e34a65262..5ad02d5e1c 100644 --- a/src/gsi/intgust.f90 +++ b/src/gsi/intgust.f90 @@ -7,8 +7,8 @@ module intgustmod ! abstract: module for intgust and its tangent linear intgust_tl ! ! program history log: -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - implemented obs adjoint test +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intgust @@ -21,15 +21,15 @@ module intgustmod ! !$$$ end documentation block -use m_obsNode , only: obsNode -use m_gustNode, only: gustNode -use m_gustNode, only: gustNode_typecast -use m_gustNode, only: gustNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode , only: obsnode +use m_gustnode, only: gustnode +use m_gustnode, only: gustnode_typecast +use m_gustnode, only: gustnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intgust +private +public intgust contains @@ -44,7 +44,7 @@ subroutine intgust(gusthead,rval,sval) ! ! program history log: ! -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - introduced ladtest_obs ! 2014-12-03 derber - modify so that use of obsdiags can be turned off ! ! input argument list: @@ -71,7 +71,7 @@ subroutine intgust(gusthead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in ) :: gusthead + class(obsnode),pointer,intent(in ) :: gusthead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -84,7 +84,7 @@ subroutine intgust(gusthead,rval,sval) real(r_kind) cg_gust,p0,grad,wnotgross,wgross,pg_gust real(r_kind),pointer,dimension(:) :: sgust real(r_kind),pointer,dimension(:) :: rgust - type(gustNode), pointer :: gustptr + type(gustnode), pointer :: gustptr ! Retrieve pointers ! Simply return if any pointer not found @@ -94,7 +94,7 @@ subroutine intgust(gusthead,rval,sval) if(ier/=0)return !gustptr => gusthead - gustptr => gustNode_typecast(gusthead) + gustptr => gustnode_typecast(gusthead) do while (associated(gustptr)) j1=gustptr%ij(1) j2=gustptr%ij(2) @@ -113,10 +113,10 @@ subroutine intgust(gusthead,rval,sval) if (lsaveobsens) then grad = val*gustptr%raterr2*gustptr%err2 !-- gustptr%diags%obssen(jiter) = grad - call obsdiagNode_set(gustptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(gustptr%diags,jiter=jiter,obssen=grad) else !-- if (gustptr%luse) gustptr%diags%tldepart(jiter)=val - if (gustptr%luse) call obsdiagNode_set(gustptr%diags,jiter=jiter,tldepart=val) + if (gustptr%luse) call obsdiagnode_set(gustptr%diags,jiter=jiter,tldepart=val) endif endif @@ -149,7 +149,7 @@ subroutine intgust(gusthead,rval,sval) endif !gustptr => gustptr%llpoint - gustptr => gustNode_nextcast(gustptr) + gustptr => gustnode_nextcast(gustptr) end do diff --git a/src/gsi/inthowv.f90 b/src/gsi/inthowv.f90 index 381baf723f..35213507f5 100644 --- a/src/gsi/inthowv.f90 +++ b/src/gsi/inthowv.f90 @@ -7,8 +7,8 @@ module inthowvmod ! abstract: module for inthowv and its tangent linear inthowv_tl ! ! program history log: -! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2012-09-14 Syed RH Rizvi, ncar/nesl/mmm/das - implemented obs adjoint test +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub inthowv @@ -21,15 +21,15 @@ module inthowvmod ! !$$$ end documentation block -use m_obsNode , only: obsNode -use m_howvNode, only: howvNode -use m_howvNode, only: howvNode_typecast -use m_howvNode, only: howvNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode , only: obsnode +use m_howvnode, only: howvnode +use m_howvnode, only: howvnode_typecast +use m_howvnode, only: howvnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC inthowv +private +public inthowv contains @@ -71,7 +71,7 @@ subroutine inthowv(howvhead,rval,sval) implicit none ! Declare passed variables - class(obsNode) , pointer,intent(in ) :: howvhead + class(obsnode) , pointer,intent(in ) :: howvhead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -84,7 +84,7 @@ subroutine inthowv(howvhead,rval,sval) real(r_kind) cg_howv,p0,grad,wnotgross,wgross,pg_howv real(r_kind),pointer,dimension(:) :: showv real(r_kind),pointer,dimension(:) :: rhowv - type(howvNode), pointer :: howvptr + type(howvnode), pointer :: howvptr ! Retrieve pointers ! Simply return if any pointer not found @@ -94,7 +94,7 @@ subroutine inthowv(howvhead,rval,sval) if(ier/=0)return !howvptr => howvhead - howvptr => howvNode_typecast(howvhead) + howvptr => howvnode_typecast(howvhead) do while (associated(howvptr)) j1=howvptr%ij(1) j2=howvptr%ij(2) @@ -113,10 +113,10 @@ subroutine inthowv(howvhead,rval,sval) if (lsaveobsens) then grad = val*howvptr%raterr2*howvptr%err2 !-- howvptr%diags%obssen(jiter) = grad - call obsdiagNode_set(howvptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(howvptr%diags,jiter=jiter,obssen=grad) else !-- if (howvptr%luse) howvptr%diags%tldepart(jiter)=val - if (howvptr%luse) call obsdiagNode_set(howvptr%diags,jiter=jiter,tldepart=val) + if (howvptr%luse) call obsdiagnode_set(howvptr%diags,jiter=jiter,tldepart=val) endif endif @@ -149,7 +149,7 @@ subroutine inthowv(howvhead,rval,sval) endif !howvptr => howvptr%llpoint - howvptr => howvNode_nextcast(howvptr) + howvptr => howvnode_nextcast(howvptr) end do diff --git a/src/gsi/intjcmod.f90 b/src/gsi/intjcmod.f90 index d264cd4b5c..9b860e3578 100644 --- a/src/gsi/intjcmod.f90 +++ b/src/gsi/intjcmod.f90 @@ -4,10 +4,10 @@ module intjcmod ! module: intjcmod module for weak constraint int routines ! pgrmmr: kleist ! -! abstract: module for Jc int routines +! abstract: module for jc int routines ! ! program history log: -! 2012-01-21 kleist - consolidation of Jc int routines into single module +! 2012-01-21 kleist - consolidation of jc int routines into single module ! 2013-10-25 todling - nullify work pointers ! 2014-03-19 pondeca - add weak constraint subroutine for wspd10m ! 2014-05-07 pondeca - add weak constraint subroutine for howv @@ -30,8 +30,8 @@ module intjcmod implicit none -PRIVATE -PUBLIC intlimqc,intlimq,intlimg,intlimp,intlimv,intlimw10m,intlimhowv,intliml,intlimcldch,intjcdfi,intjcpdry,intjcpdry1,intjcpdry2 +private +public intlimqc,intlimq,intlimg,intlimp,intlimv,intlimw10m,intlimhowv,intliml,intlimcldch,intjcdfi,intjcpdry,intjcpdry1,intjcpdry2 contains @@ -80,14 +80,14 @@ subroutine intlimq(rval,sval,itbin) ! Declare passed variables type(gsi_bundle),intent(in ) :: sval type(gsi_bundle),intent(inout) :: rval - integer, intent(in) :: itbin + integer(i_kind),intent(in ) :: itbin ! Declare local variables integer(i_kind) i,j,k,ier,istatus,ii,mm1 real(r_kind) q - real(r_kind),pointer,dimension(:,:,:) :: sq=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rq=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ges_q_it=>NULL() + real(r_kind),pointer,dimension(:,:,:) :: sq=>null() + real(r_kind),pointer,dimension(:,:,:) :: rq=>null() + real(r_kind),pointer,dimension(:,:,:) :: ges_q_it=>null() if (factqmin==zero .and. factqmax==zero) return @@ -169,9 +169,9 @@ subroutine intlimqc(rval,sval,itbin,cldtype) integer(i_kind) i,j,k,ier,ier1,istatus real(r_kind) qc real(r_kind) factqc - real(r_kind),pointer,dimension(:,:,:) :: sqc=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqc=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: ges_qc_it=>NULL() + real(r_kind),pointer,dimension(:,:,:) :: sqc=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqc=>null() + real(r_kind),pointer,dimension(:,:,:) :: ges_qc_it=>null() ier=0 ier1=0 @@ -263,8 +263,8 @@ subroutine intlimg(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) gust - real(r_kind),pointer,dimension(:,:) :: sg=>NULL() - real(r_kind),pointer,dimension(:,:) :: rg=>NULL() + real(r_kind),pointer,dimension(:,:) :: sg=>null() + real(r_kind),pointer,dimension(:,:) :: rg=>null() if (factg==zero) return @@ -325,8 +325,8 @@ subroutine intlimp(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) pblh - real(r_kind),pointer,dimension(:,:) :: sp=>NULL() - real(r_kind),pointer,dimension(:,:) :: rp=>NULL() + real(r_kind),pointer,dimension(:,:) :: sp=>null() + real(r_kind),pointer,dimension(:,:) :: rp=>null() if (factp==zero) return @@ -387,8 +387,8 @@ subroutine intlimv(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) vis - real(r_kind),pointer,dimension(:,:) :: sv=>NULL() - real(r_kind),pointer,dimension(:,:) :: rv=>NULL() + real(r_kind),pointer,dimension(:,:) :: sv=>null() + real(r_kind),pointer,dimension(:,:) :: rv=>null() if (factv==zero) return @@ -453,8 +453,8 @@ subroutine intlimw10m(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) wspd10m - real(r_kind),pointer,dimension(:,:) :: sg=>NULL() - real(r_kind),pointer,dimension(:,:) :: rg=>NULL() + real(r_kind),pointer,dimension(:,:) :: sg=>null() + real(r_kind),pointer,dimension(:,:) :: rg=>null() if (factw10m==zero) return @@ -519,8 +519,8 @@ subroutine intlimhowv(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) howv - real(r_kind),pointer,dimension(:,:) :: sg=>NULL() - real(r_kind),pointer,dimension(:,:) :: rg=>NULL() + real(r_kind),pointer,dimension(:,:) :: sg=>null() + real(r_kind),pointer,dimension(:,:) :: rg=>null() if (facthowv==zero) return @@ -581,8 +581,8 @@ subroutine intliml(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) lcbas - real(r_kind),pointer,dimension(:,:) :: sv=>NULL() - real(r_kind),pointer,dimension(:,:) :: rv=>NULL() + real(r_kind),pointer,dimension(:,:) :: sv=>null() + real(r_kind),pointer,dimension(:,:) :: rv=>null() if (factl==zero) return @@ -647,8 +647,8 @@ subroutine intlimcldch(rval,sval) ! Declare local variables integer(i_kind) i,j,ier,istatus real(r_kind) cldch - real(r_kind),pointer,dimension(:,:) :: sg=>NULL() - real(r_kind),pointer,dimension(:,:) :: rg=>NULL() + real(r_kind),pointer,dimension(:,:) :: sg=>null() + real(r_kind),pointer,dimension(:,:) :: rg=>null() if (factcldch==zero) return @@ -725,22 +725,22 @@ subroutine intjcpdry(rval,sval,nbins,pjc) real(r_quad),dimension(nsig) :: mass2 real(r_quad) rcon,con,dmass integer(i_kind) i,j,k,it,ii,mm1,ier,icw,iql,iqi,istatus - real(r_kind),pointer,dimension(:,:,:) :: sq =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sc =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sql=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqi=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rq =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rc =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rql=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqi=>NULL() - real(r_kind),pointer,dimension(:,:) :: sp =>NULL() - real(r_kind),pointer,dimension(:,:) :: rp =>NULL() + real(r_kind),pointer,dimension(:,:,:) :: sq =>null() + real(r_kind),pointer,dimension(:,:,:) :: sc =>null() + real(r_kind),pointer,dimension(:,:,:) :: sql=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqi=>null() + real(r_kind),pointer,dimension(:,:,:) :: rq =>null() + real(r_kind),pointer,dimension(:,:,:) :: rc =>null() + real(r_kind),pointer,dimension(:,:,:) :: rql=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqi=>null() + real(r_kind),pointer,dimension(:,:) :: sp =>null() + real(r_kind),pointer,dimension(:,:) :: rp =>null() integer(i_kind) :: n it=ntguessig mass=zero_quad - rcon=one_quad/(two_quad*float(nlon)) + rcon=one_quad/(two_quad*real(nlon,r_quad)) mm1=mype+1 do n=1,nbins @@ -753,17 +753,17 @@ subroutine intjcpdry(rval,sval,nbins,pjc) call gsi_bundlegetpointer(sval(n),'qi',sqi,istatus);iqi=istatus+iqi call gsi_bundlegetpointer(sval(n),'ps',sp, istatus);ier=istatus+ier if(ier+icw*(iql+iqi)/=0)then - if (mype==0) write(6,*)'intjcpdry: checking ier+icw*(iql+iqi)=', ier+icw*(iql+iqi) - return + if (mype==0) write(6,*)'intjcpdry: checking ier+icw*(iql+iqi)=', ier+icw*(iql+iqi) + return end if ! Calculate mean surface pressure contribution in subdomain do j=2,lon2-1 - do i=2,lat2-1 - ii=istart(mm1)+i-2 - mass(n)=mass(n)+sp(i,j)*wgtlats(ii) - end do + do i=2,lat2-1 + ii=istart(mm1)+i-2 + mass(n)=mass(n)+sp(i,j)*wgtlats(ii) + end do end do mass2(:)=zero_quad @@ -788,7 +788,7 @@ subroutine intjcpdry(rval,sval,nbins,pjc) end do end do -! First, use MPI to get global mean increment +! First, use mpi to get global mean increment call mpl_allreduce(2*nbins,qpvals=mass) do n=1,nbins @@ -799,8 +799,8 @@ subroutine intjcpdry(rval,sval,nbins,pjc) call gsi_bundlegetpointer(rval(n),'qi',rqi,istatus);iqi=istatus+iqi call gsi_bundlegetpointer(rval(n),'ps',rp, istatus);ier=istatus+ier if(ier+icw*(iql+iqi)/=0)then - if (mype==0) write(6,*)'intjcpdry: checking ier+icw*(iql+iqi)=', ier+icw*(iql+iqi) - return + if (mype==0) write(6,*)'intjcpdry: checking ier+icw*(iql+iqi)=', ier+icw*(iql+iqi) + return end if ! Remove water-vapor contribution to get incremental dry ps ! if (mype==0) write(6,*)'intjcpdry: total mass =', mass(n) @@ -813,10 +813,10 @@ subroutine intjcpdry(rval,sval,nbins,pjc) ! Calculate mean surface pressure contribution in subdomain do j=2,lon2-1 - do i=2,lat2-1 - ii=istart(mm1)+i-2 - rp(i,j)=rp(i,j)+dmass*wgtlats(ii) - end do + do i=2,lat2-1 + ii=istart(mm1)+i-2 + rp(i,j)=rp(i,j)+dmass*wgtlats(ii) + end do end do ! Remove water to get incremental dry ps !$omp parallel do schedule(dynamic,1) private(k,j,i,ii,con) @@ -887,21 +887,21 @@ subroutine intjcpdry1(sval,nbins,mass) real(r_quad) rcon,con integer(i_kind) i,j,k,it,ii,mm1,icw,iql,iqi integer(i_kind) iq,iqr,iqs,iqg,iqh,ips - real(r_kind),pointer,dimension(:,:,:) :: sq =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sc =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sql=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqi=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqr=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqs=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqg=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: sqh=>NULL() - real(r_kind),pointer,dimension(:,:) :: sp =>NULL() + real(r_kind),pointer,dimension(:,:,:) :: sq =>null() + real(r_kind),pointer,dimension(:,:,:) :: sc =>null() + real(r_kind),pointer,dimension(:,:,:) :: sql=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqi=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqr=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqs=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqg=>null() + real(r_kind),pointer,dimension(:,:,:) :: sqh=>null() + real(r_kind),pointer,dimension(:,:) :: sp =>null() integer(i_kind) :: n it=ntguessig mass=zero_quad - rcon=one_quad/(two_quad*float(nlon)) + rcon=one_quad/(two_quad*real(nlon,r_quad)) mm1=mype+1 do n=1,nbins @@ -918,17 +918,17 @@ subroutine intjcpdry1(sval,nbins,mass) call gsi_bundlegetpointer(sval(n),'qh',sqh, iqh ) call gsi_bundlegetpointer(sval(n),'ps',sp, ips ) if ( iq*ips/=0 .or. icw*(iql+iqi)/=0 ) then - if (mype==0) write(6,*)'intjcpdry1: warning - missing some required variables' - if (mype==0) write(6,*)'intjcpdry1: constraint for dry mass constraint not performed' - return + if (mype==0) write(6,*)'intjcpdry1: warning - missing some required variables' + if (mype==0) write(6,*)'intjcpdry1: constraint for dry mass constraint not performed' + return end if ! Calculate mean surface pressure contribution in subdomain do j=2,lon2-1 - do i=2,lat2-1 - ii=istart(mm1)+i-2 - mass(n)=mass(n)+sp(i,j)*wgtlats(ii) - end do + do i=2,lat2-1 + ii=istart(mm1)+i-2 + mass(n)=mass(n)+sp(i,j)*wgtlats(ii) + end do end do mass2(:)=zero_quad @@ -1010,20 +1010,20 @@ subroutine intjcpdry2(rval,nbins,mass,pjc) real(r_quad) rcon,con,dmass integer(i_kind) i,j,k,it,ii,mm1,icw,iql,iqi integer(i_kind) iq,iqr,iqs,iqg,iqh,ips - real(r_kind),pointer,dimension(:,:,:) :: rq =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rc =>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rql=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqi=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqr=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqs=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqg=>NULL() - real(r_kind),pointer,dimension(:,:,:) :: rqh=>NULL() - real(r_kind),pointer,dimension(:,:) :: rp =>NULL() + real(r_kind),pointer,dimension(:,:,:) :: rq =>null() + real(r_kind),pointer,dimension(:,:,:) :: rc =>null() + real(r_kind),pointer,dimension(:,:,:) :: rql=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqi=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqr=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqs=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqg=>null() + real(r_kind),pointer,dimension(:,:,:) :: rqh=>null() + real(r_kind),pointer,dimension(:,:) :: rp =>null() integer(i_kind) :: n it=ntguessig - rcon=one_quad/(two_quad*float(nlon)) + rcon=one_quad/(two_quad*real(nlon,r_quad)) mm1=mype+1 do n=1,nbins @@ -1038,9 +1038,9 @@ subroutine intjcpdry2(rval,nbins,mass,pjc) call gsi_bundlegetpointer(rval(n),'qh',rqh, iqh ) call gsi_bundlegetpointer(rval(n),'ps',rp, ips ) if( iq*ips /=0 .or. icw*(iql+iqi) /=0 ) then - if (mype==0) write(6,*)'intjcpdry2: warning - missing some required variables' - if (mype==0) write(6,*)'intjcpdry2: constraint for dry mass constraint not performed' - return + if (mype==0) write(6,*)'intjcpdry2: warning - missing some required variables' + if (mype==0) write(6,*)'intjcpdry2: constraint for dry mass constraint not performed' + return end if ! Remove water-vapor contribution to get incremental dry ps ! if (mype==0) write(6,*)'intjcpdry: total mass =', mass(n) @@ -1053,10 +1053,10 @@ subroutine intjcpdry2(rval,nbins,mass,pjc) ! Calculate mean surface pressure contribution in subdomain do j=2,lon2-1 - do i=2,lat2-1 - ii=istart(mm1)+i-2 - rp(i,j)=rp(i,j)+dmass*wgtlats(ii) - end do + do i=2,lat2-1 + ii=istart(mm1)+i-2 + rp(i,j)=rp(i,j)+dmass*wgtlats(ii) + end do end do ! Remove water to get incremental dry ps !$omp parallel do schedule(dynamic,1) private(k,j,i,ii,con) @@ -1087,7 +1087,7 @@ end subroutine intjcpdry2 subroutine intjcdfi(rval,sval,pjc) !$$$ subprogram documentation block ! . . . . -! subprogram: intjcdfi calculate Jc DFI terms and contribution to gradient +! subprogram: intjcdfi calculate jc dfi terms and contribution to gradient ! prgmmr: tremolet ! ! program history log: @@ -1113,23 +1113,23 @@ subroutine intjcdfi(rval,sval,pjc) ! !$$$ end documentation block -use jcmod, only: wgtdfi,alphajc -use gsi_4dvar, only: nobs_bins -use mpimod, only: mype -use state_vectors, only : allocate_state,deallocate_state -use gsi_bundlemod, only : self_add,self_mul,assignment(=) -implicit none + use jcmod, only: wgtdfi,alphajc + use gsi_4dvar, only: nobs_bins + use mpimod, only: mype + use state_vectors, only : allocate_state,deallocate_state + use gsi_bundlemod, only : self_add,self_mul,assignment(=) + implicit none ! Declare passed variables -type(gsi_bundle),dimension(nobs_bins),intent(in ) :: sval -type(gsi_bundle),dimension(nobs_bins),intent(inout) :: rval -real(r_quad), optional, intent( out) :: pjc + type(gsi_bundle),dimension(nobs_bins),intent(in ) :: sval + type(gsi_bundle),dimension(nobs_bins),intent(inout) :: rval + real(r_quad), optional, intent( out) :: pjc ! Declare local variables -integer(i_kind) :: jj,idfi -real(r_quad),parameter :: half_quad=0.5_r_quad -type(gsi_bundle) :: sfilter,afilter -real(r_quad) :: cost + integer(i_kind) :: jj,idfi + real(r_quad),parameter :: half_quad=0.5_r_quad + type(gsi_bundle) :: sfilter,afilter + real(r_quad) :: cost !************************************************************************************ @@ -1146,11 +1146,11 @@ subroutine intjcdfi(rval,sval,pjc) ! Compute difference from filtered state call self_add(sfilter,-one,sval(idfi)) -! Apply Jc multiplicative factor +! Apply jc multiplicative factor call self_mul(sfilter,alphajc) -! Compute Jc (norm of difference) -! Jc = 1/2 * wgt * sfilter *sfilter +! Compute jc (norm of difference) +! jc = 1/2 * wgt * sfilter *sfilter ! afilter = wgt * sfilter call enorm_state(sfilter,cost,afilter) if(present(pjc))then @@ -1158,7 +1158,7 @@ subroutine intjcdfi(rval,sval,pjc) if (mype==0) write(6,*)'Jc DFI=',pjc endif -! Adjoint Jc multiplicative factor +! Adjoint jc multiplicative factor call self_mul(afilter,alphajc) ! Adjoint of difference from filtered state diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90 index e514a38a22..9c84a9e047 100644 --- a/src/gsi/intjo.f90 +++ b/src/gsi/intjo.f90 @@ -4,16 +4,16 @@ module intjomod ! module: intjo module for intjo ! prgmmr: ! -! abstract: module for H'R^{-1}H +! abstract: module for h'r^{-1}h ! ! program history log: ! 2008-12-01 Todling - wrap in module ! 2009-08-13 lueken - update documentation -! 2015-09-03 guo - obsmod::obs_handle has been replaced with obsHeadBundle, -! defined by m_obsHeadBundle. +! 2015-09-03 guo - obsmod::obs_handle has been replaced with obsheadbundle, +! defined by m_obsheadbundle. ! 2016-08-29 J Guo - Separated calls to intozlay() and intozlev() ! 2018-08-10 guo - a new generic intjo() implementation replaced type -! specific intXYZ() calls with polymorphic %intjo(). +! specific intxyz() calls with polymorphic %intjo(). ! ! subroutines included: ! sub intjo_ @@ -26,30 +26,30 @@ module intjomod ! !$$$ end documentation block -use gsi_obOperTypeManager, only: obOper_count -use gsi_obOperTypeManager, only: obOper_typeInfo -use gsi_obOperTypeManager, only: & - iobOper_t, iobOper_pw, iobOper_q, & - iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, & - iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & - iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & - iobOper_gpsbend, iobOper_gpsref, & - iobOper_rad, iobOper_pcp, iobOper_aero, iobOper_gust, & - iobOper_vis, iobOper_pblh, iobOper_wspd10m, iobOper_td2m, iobOper_mxtm, & - iobOper_mitm, iobOper_pmsl, iobOper_howv, iobOper_tcamt, iobOper_lcbas, & - iobOper_cldch, iobOper_uwnd10m, iobOper_vwnd10m, iobOper_swcp, iobOper_lwcp, & - iobOper_light +use gsi_obopertypemanager, only: oboper_count +use gsi_obopertypemanager, only: oboper_typeinfo +use gsi_obopertypemanager, only: & + ioboper_t, ioboper_pw, ioboper_q, & + ioboper_cldtot, ioboper_w, ioboper_dw, & + ioboper_rw, ioboper_dbz, & + ioboper_spd, ioboper_oz, ioboper_o3l, ioboper_colvk, & + ioboper_pm2_5, ioboper_pm10, ioboper_ps, ioboper_tcp, ioboper_sst, & + ioboper_gpsbend, ioboper_gpsref, & + ioboper_rad, ioboper_pcp, ioboper_aero, ioboper_gust, & + ioboper_vis, ioboper_pblh, ioboper_wspd10m, ioboper_td2m, ioboper_mxtm, & + ioboper_mitm, ioboper_pmsl, ioboper_howv, ioboper_tcamt, ioboper_lcbas, & + ioboper_cldch, ioboper_uwnd10m, ioboper_vwnd10m, ioboper_swcp, ioboper_lwcp, & + ioboper_light use kinds, only: i_kind use mpeu_util, only: perr,die implicit none -PRIVATE -PUBLIC intjo +private +public intjo interface intjo; module procedure & - intjo_, intjo_reduced_ + intjo_, intjo_reduced_ end interface ! This is a mapping to the exact %intjo() calling sequence in the earlier @@ -57,18 +57,18 @@ module intjomod ! ordering for rval and qpred. It is not necessary, and can be removed once ! some non-zero-diff modifications are introduced. -integer(i_kind),parameter,dimension(obOper_count):: ix_obtype = (/ & - iobOper_t, iobOper_pw, iobOper_q, & - iobOper_cldtot, iobOper_w, iobOper_dw, & - iobOper_rw, iobOper_dbz, & - iobOper_spd, iobOper_oz, iobOper_o3l, iobOper_colvk, & - iobOper_pm2_5, iobOper_pm10, iobOper_ps, iobOper_tcp, iobOper_sst, & - iobOper_gpsbend, iobOper_gpsref, & - iobOper_rad, iobOper_pcp, iobOper_aero, iobOper_gust, & - iobOper_vis, iobOper_pblh, iobOper_wspd10m, iobOper_td2m, iobOper_mxtm, & - iobOper_mitm, iobOper_pmsl, iobOper_howv, iobOper_tcamt, iobOper_lcbas, & - iobOper_cldch, iobOper_uwnd10m, iobOper_vwnd10m, iobOper_swcp, iobOper_lwcp, & - iobOper_light /) +integer(i_kind),parameter,dimension(oboper_count):: ix_obtype = (/ & + ioboper_t, ioboper_pw, ioboper_q, & + ioboper_cldtot, ioboper_w, ioboper_dw, & + ioboper_rw, ioboper_dbz, & + ioboper_spd, ioboper_oz, ioboper_o3l, ioboper_colvk, & + ioboper_pm2_5, ioboper_pm10, ioboper_ps, ioboper_tcp, ioboper_sst, & + ioboper_gpsbend, ioboper_gpsref, & + ioboper_rad, ioboper_pcp, ioboper_aero, ioboper_gust, & + ioboper_vis, ioboper_pblh, ioboper_wspd10m, ioboper_td2m, ioboper_mxtm, & + ioboper_mitm, ioboper_pmsl, ioboper_howv, ioboper_tcamt, ioboper_lcbas, & + ioboper_cldch, ioboper_uwnd10m, ioboper_vwnd10m, ioboper_swcp, ioboper_lwcp, & + ioboper_light /) !...|....1....|....2....|....3....|....4....|....5....|....6....|....7....|....8....|....9....|....0 character(len=*),parameter:: myname="intjomod" @@ -79,36 +79,36 @@ subroutine intjo_(rval,qpred,sval,sbias) !$$$ subprogram documentation block ! . . . . -! subprogram: intjo calculate RHS for analysis equation +! subprogram: intjo calculate rhs for analysis equation ! prgmmr: derber org: np23 date: 2003-12-18 ! -! abstract: calculate RHS for all variables (nonlinear qc version) +! abstract: calculate rhs for all variables (nonlinear qc version) ! ! A description of nonlinear qc follows: ! -! The observation penalty Jo is defined as +! The observation penalty jo is defined as ! -! Jo = - (sum over obs) 2*log(Po) +! jo = - (sum over obs) 2*log(po) ! ! where, ! -! Po = Wnotgross*exp(-.5*(Hn(x+xb) - yo)**2 ) + Wgross +! po = wnotgross*exp(-.5*(hn(x+xb) - yo)**2 ) + wgross ! with -! Hn = the forward model (possibly non-linear) normalized by +! hn = the forward model (possibly non-linear) normalized by ! observation error ! x = the current estimate of the analysis increment ! xb = the background state ! yo = the observation normalized by observation error ! -! Note: The factor 2 in definition of Jo is present because the -! penalty Jo as used in this code is 2*(usual definition +! Note: The factor 2 in definition of jo is present because the +! penalty jo as used in this code is 2*(usual definition ! of penalty) ! -! Wgross = Pgross*cg +! wgross = pgross*cg ! -! Wnotgross = 1 - Wgross +! wnotgross = 1 - wgross ! -! Pgross = probability of gross error for observation (assumed +! pgross = probability of gross error for observation (assumed ! here to have uniform distribution over the possible ! range of values) ! @@ -117,7 +117,7 @@ subroutine intjo_(rval,qpred,sval,sbias) ! b = possible range of variable for gross errors, normalized by ! observation error ! -! The values for the above parameters that Bill Collins used in the +! The values for the above parameters that bill collins used in the ! eta 3dvar are: ! ! cg = cg_term/b, where cg_term = sqrt(2*pi)/2 @@ -132,17 +132,17 @@ subroutine intjo_(rval,qpred,sval,sbias) ! pg_rad=.002 ! probability of gross error for radiances ! ! -! Given the above Jo, the gradient of Jo is as follows: +! Given the above jo, the gradient of jo is as follows: ! -! T -! gradx(Jo) = - (sum over observations) 2*H (Hn(x+xb)-yo)*(Po - Wgross)/Po +! t +! gradx(jo) = - (sum over observations) 2*h (hn(x+xb)-yo)*(po - wgross)/po ! ! where, ! -! H = tangent linear model of Hn about x+xb +! h = tangent linear model of hn about x+xb ! ! -! Note that if Pgross = 0.0, then Wnotgross=1.0 and Wgross=0.0. That is, +! Note that if pgross = 0.0, then wnotgross=1.0 and wgross=0.0. That is, ! the code runs as though nonlinear quality control were not present ! (which is indeed the case since the gross error probability is 0). ! @@ -159,7 +159,7 @@ subroutine intjo_(rval,qpred,sval,sbias) ! for wind components into int for st,vp ! 2004-11-30 treadon - add brightness temperatures to nonlinear ! quality control -! 2004-12-03 treadon - replace mpe_iallreduce (IBM extension) with +! 2004-12-03 treadon - replace mpe_iallreduce (ibm extension) with ! standard mpi_allreduce ! 2005-01-20 okamoto - add u,v to intrad ! 2005-02-23 wu - changes related to normalized rh option @@ -169,15 +169,15 @@ subroutine intjo_(rval,qpred,sval,sbias) ! for 2dvar only surface analysis option ! 2005-06-03 parrish - add horizontal derivatives ! 2005-07-10 kleist - add dynamic constraint term -! 2005-09-29 kleist - expand Jc term, include time derivatives vector -! 2005-11-21 kleist - separate tendencies from Jc term, add call to calctends adjoint -! 2005-12-01 cucurull - add code for GPS local bending angle, add use obsmod for ref_obs +! 2005-09-29 kleist - expand jc term, include time derivatives vector +! 2005-11-21 kleist - separate tendencies from jc term, add call to calctends adjoint +! 2005-12-01 cucurull - add code for gps local bending angle, add use obsmod for ref_obs ! 2005-12-20 parrish - add arguments to call to intt to allow for option of using boundary ! layer forward tlm. ! 2006-02-03 derber - modify to increase reproducibility ! 2006-03-17 park - correct error in call to intt--rval,sval --> rvaluv,svaluv ! in order to correctly pass wind variables. -! 2006-04-06 kleist - include both Jc formulations +! 2006-04-06 kleist - include both jc formulations ! 2006-07-26 parrish - correct inconsistency in computation of space and time derivatives of q ! currently, if derivatives computed, for q it is normalized q, but ! should be mixing ratio. @@ -185,7 +185,7 @@ subroutine intjo_(rval,qpred,sval,sbias) ! 2007-03-19 tremolet - binning of observations ! 2007-04-13 tremolet - split jo from other components of intall ! 2007-06-04 derber - use quad precision to get reproducibility over number of processors -! 2008-11-27 todling - add tendencies for FOTO support and new interface to int's +! 2008-11-27 todling - add tendencies for foto support and new interface to int's ! 2009-01-08 todling - remove reference to ozohead ! 2009-03-23 meunier - Add call to intlag (lagrangian observations) ! 2010-01-11 zhang,b - Bug fix: bias predictors need to be accumulated over nbins @@ -212,7 +212,7 @@ subroutine intjo_(rval,qpred,sval,sbias) ! qpred ! ! output argument list: -! rval - RHS on grid +! rval - rhs on grid ! qpred ! ! remarks: @@ -221,86 +221,86 @@ subroutine intjo_(rval,qpred,sval,sbias) ! ! 2) The two interfaces to the int-routines should be temporary. ! In the framework of the 4dvar-code, foto can be re-implemented as -! an approximate M and M' to the model matrices in 4dvar. Once that +! an approximate m and m' to the model matrices in 4dvar. Once that ! is done, the int-routines should no longer need the time derivatives. -! (Todling) +! (todling) ! 3) Notice that now (2010-05-13) int routines handle non-essential ! variables internally; also, when pointers non-existent, int routines -! simply return (Todling). +! simply return (todling). ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ -use kinds, only: i_kind,r_quad -use gsi_bundlemod, only: gsi_bundle -use bias_predictors, only: predictors -use m_obsdiags, only: obOper_create -use m_obsdiags, only: obOper_destroy -use gsi_obOper, only: obOper + use kinds, only: i_kind,r_quad + use gsi_bundlemod, only: gsi_bundle + use bias_predictors, only: predictors + use m_obsdiags, only: oboper_create + use m_obsdiags, only: oboper_destroy + use gsi_oboper, only: oboper -use intradmod, only: setrad + use intradmod, only: setrad -implicit none + implicit none ! Declare passed variables -type(gsi_bundle), dimension( :), intent(inout) :: rval ! (nobs_bins) -type(gsi_bundle), dimension( :), intent(in ) :: sval ! (nobs_bins) -real(r_quad ), dimension(:,:), intent(inout) :: qpred ! (:,nobs_bins) -type(predictors), intent(in ) :: sbias + type(gsi_bundle), dimension( :), intent(inout) :: rval ! (nobs_bins) + type(gsi_bundle), dimension( :), intent(in ) :: sval ! (nobs_bins) + real(r_quad ), dimension(:,:), intent(inout) :: qpred ! (:,nobs_bins) + type(predictors), intent(in ) :: sbias -character(len=*),parameter:: myname_=myname//"::intjo_" + character(len=*),parameter:: myname_=myname//"::intjo_" ! Declare local variables -integer(i_kind):: ibin,it,ix -class(obOper),pointer:: it_obOper + integer(i_kind):: ibin,it,ix + class(oboper),pointer:: it_oboper !****************************************************************************** call setrad(sval(1)) -! "RHS for jo", as it was labeled in intall(). -!$omp parallel do schedule(dynamic,1) private(ibin,it,ix,it_obOper) +! "rhs for jo", as it was labeled in intall(). +!$omp parallel do schedule(dynamic,1) private(ibin,it,ix,it_oboper) do ibin=1,size(sval) - do it=1,obOper_count - !ix=ix_obtype(it) ! Use this line to ensure the same jo summartion - ! sequence as intjo was in its early implementation, - ! for reproducibility. - - ix=it ! Using this line, jo summation sequence is not the same as - ! it used to be, nor the same if someone chooses to change - ! enumration sequence of obOpers in gsi_obOperTypeManager.F90. - ! But it would make this code more portable to new obOper - ! extensions. - - it_obOper => obOper_create(ix) - - if(.not.associated(it_obOper)) then - call perr(myname_,'unexpected obOper, associated(it_obOper) =',associated(it_obOper)) - call perr(myname_,' obOper_typeInfo(ioper) =',obOper_typeInfo(ix)) - call perr(myname_,' ioper =',ix) - call perr(myname_,' it =',it) - call perr(myname_,' obOper_count =',obOper_count) - call perr(myname_,' ibin =',ibin) - call die(myname_) + do it=1,oboper_count + !ix=ix_obtype(it) ! Use this line to ensure the same jo summartion + ! sequence as intjo was in its early implementation, + ! for reproducibility. + + ix=it ! Using this line, jo summation sequence is not the same as + ! it used to be, nor the same if someone chooses to change + ! enumration sequence of obopers in gsi_obopertypemanager.F90. + ! But it would make this code more portable to new oboper + ! extensions. + + it_oboper => oboper_create(ix) + + if(.not.associated(it_oboper)) then + call perr(myname_,'unexpected obOper, associated(it_obOper) =',associated(it_oboper)) + call perr(myname_,' obOper_typeInfo(ioper) =',oboper_typeinfo(ix)) + call perr(myname_,' ioper =',ix) + call perr(myname_,' it =',it) + call perr(myname_,' obOper_count =',oboper_count) + call perr(myname_,' ibin =',ibin) + call die(myname_) endif - if(.not.associated(it_obOper%obsLL)) then - call perr(myname_,'unexpected component, associated(%obsLL) =',associated(it_obOper%obsLL)) - call perr(myname_,' obOper_typeInfo(ioper) =',obOper_typeInfo(ix)) - call perr(myname_,' ioper =',ix) - call perr(myname_,' it =',it) - call perr(myname_,' obOper_count =',obOper_count) - call perr(myname_,' ibin =',ibin) - call die(myname_) + if(.not.associated(it_oboper%obsll)) then + call perr(myname_,'unexpected component, associated(%obsLL) =',associated(it_oboper%obsll)) + call perr(myname_,' obOper_typeInfo(ioper) =',oboper_typeinfo(ix)) + call perr(myname_,' ioper =',ix) + call perr(myname_,' it =',it) + call perr(myname_,' obOper_count =',oboper_count) + call perr(myname_,' ibin =',ibin) + call die(myname_) endif - call it_obOper%intjo(ibin,rval(ibin),sval(ibin),qpred(:,ibin),sbias) - call obOper_destroy(it_obOper) - enddo + call it_oboper%intjo(ibin,rval(ibin),sval(ibin),qpred(:,ibin),sbias) + call oboper_destroy(it_oboper) + enddo end do -return + return end subroutine intjo_ subroutine intjo_reduced_(rval,qpred,sval,sbias) @@ -326,7 +326,7 @@ subroutine intjo_reduced_(rval,qpred,sval,sbias) call intjo_(rval,qpred_bin,sval,sbias) do ibin=1,size(rval) - qpred(:)=qpred(:)+qpred_bin(:,ibin) + qpred(:)=qpred(:)+qpred_bin(:,ibin) enddo deallocate(qpred_bin) diff --git a/src/gsi/intlag.f90 b/src/gsi/intlag.f90 index de2738df5b..759de6ff1e 100644 --- a/src/gsi/intlag.f90 +++ b/src/gsi/intlag.f90 @@ -11,7 +11,7 @@ module intlagmod ! 2008-03-23 lmeunier - initial code ! 2009-08-13 lueken - update documentation ! 2011-08-01 lueken - changed F90 to f90 -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intlag @@ -23,15 +23,15 @@ module intlagmod ! machine: ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_lagNode, only: lagNode -use m_lagNode, only: lagNode_typecast -use m_lagNode, only: lagNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_lagnode, only: lagnode +use m_lagnode, only: lagnode_typecast +use m_lagnode, only: lagnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intlag +private +public intlag contains @@ -72,7 +72,7 @@ subroutine intlag(laghead,rval,sval,obsbin) use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer - use lag_fields, only: lag_gather_stateuv,lag_ADscatter_stateuv + use lag_fields, only: lag_gather_stateuv,lag_adscatter_stateuv use lag_fields, only: lag_u_full,lag_v_full,lag_uv_alloc,lag_uv_fill use lag_fields, only: lag_tl_vec,lag_ad_vec use lag_fields, only: orig_lag_num,ntotal_orig_lag,lag_kcount @@ -83,7 +83,7 @@ subroutine intlag(laghead,rval,sval,obsbin) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in ) :: laghead + class(obsnode),pointer,intent(in ) :: laghead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval integer(i_kind), intent(in ) :: obsbin @@ -97,7 +97,7 @@ subroutine intlag(laghead,rval,sval,obsbin) real(r_kind) grad_lon,grad_lat,grad_p real(r_kind),pointer,dimension(:) :: su,sv real(r_kind),pointer,dimension(:) :: ru,rv - type(lagNode), pointer :: lagptr + type(lagnode), pointer :: lagptr integer(i_kind) i,j,ier,istatus real(r_kind),dimension(:,:),allocatable:: adu_tmp,adv_tmp @@ -123,7 +123,7 @@ subroutine intlag(laghead,rval,sval,obsbin) end if end if - ! Allocate AD fields + ! Allocate ad fields if (l_do_adjoint) then allocate(adu_tmp(iglobal,lag_kcount)) allocate(adv_tmp(iglobal,lag_kcount)) @@ -131,7 +131,7 @@ subroutine intlag(laghead,rval,sval,obsbin) end if !lagptr => laghead - lagptr => lagNode_typecast(laghead) + lagptr => lagnode_typecast(laghead) do while (associated(lagptr)) ! Forward model @@ -148,8 +148,8 @@ subroutine intlag(laghead,rval,sval,obsbin) print *,'MAX INC V',maxval(abs(lag_v_full(:,:,obsbin))) end if call lag_rk2iter_tl(lagptr%speci,lagptr%specr,& - &lon_tl,lat_tl,p_tl,& - &lag_u_full(:,:,obsbin),lag_v_full(:,:,obsbin)) + lon_tl,lat_tl,p_tl,& + lag_u_full(:,:,obsbin),lag_v_full(:,:,obsbin)) if (iv_debug>=1) print *,'TL correction:',lon_tl*rad2deg,lat_tl*rad2deg if (lsaveobsens) then @@ -157,13 +157,13 @@ subroutine intlag(laghead,rval,sval,obsbin) grad_lat = lat_tl*lagptr%raterr2*lagptr%err2_lat !-- lagptr%diag_lon%obssen(jiter) = lon_tl*lagptr%raterr2*lagptr%err2_lon !-- lagptr%diag_lat%obssen(jiter) = lat_tl*lagptr%raterr2*lagptr%err2_lat - call obsdiagNode_set(lagptr%diag_lon,jiter=jiter,obssen=grad_lon) - call obsdiagNode_set(lagptr%diag_lat,jiter=jiter,obssen=grad_lat) + call obsdiagnode_set(lagptr%diag_lon,jiter=jiter,obssen=grad_lon) + call obsdiagnode_set(lagptr%diag_lat,jiter=jiter,obssen=grad_lat) else !-- if (lagptr%luse) lagptr%diag_lon%tldepart(jiter)=lon_tl !-- if (lagptr%luse) lagptr%diag_lat%tldepart(jiter)=lat_tl - if (lagptr%luse) call obsdiagNode_set(lagptr%diag_lon,jiter=jiter,tldepart=lon_tl) - if (lagptr%luse) call obsdiagNode_set(lagptr%diag_lat,jiter=jiter,tldepart=lat_tl) + if (lagptr%luse) call obsdiagnode_set(lagptr%diag_lon,jiter=jiter,tldepart=lon_tl) + if (lagptr%luse) call obsdiagnode_set(lagptr%diag_lat,jiter=jiter,tldepart=lat_tl) endif if (l_do_adjoint) then @@ -181,7 +181,7 @@ subroutine intlag(laghead,rval,sval,obsbin) wnotgross= one-lagptr%pg wgross = lagptr%pg*cg_lag/wnotgross p0 = wgross/(wgross+& - &exp(-half*(lagptr%err2_lon*lon_tl**2+lagptr%err2_lat*lat_tl**2))) + exp(-half*(lagptr%err2_lon*lon_tl**2+lagptr%err2_lat*lat_tl**2))) lon_tl = lon_tl*(one-p0) lat_tl = lat_tl*(one-p0) if (iv_debug>=1) print *,'Do nlnqc_iter' @@ -199,13 +199,13 @@ subroutine intlag(laghead,rval,sval,obsbin) ! Adjoint call lag_rk2iter_ad(lagptr%speci,lagptr%specr,& - &grad_lon,grad_lat,grad_p,adu_tmp,adv_tmp) + grad_lon,grad_lat,grad_p,adu_tmp,adv_tmp) lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)=& - &lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)+grad_lon + lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)+grad_lon lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)=& - &lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)+grad_lat + lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)+grad_lat lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)=& - &lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)+grad_p + lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)+grad_p if (iv_debug>=2) then do i=1,iglobal do j=1,lag_kcount @@ -218,15 +218,15 @@ subroutine intlag(laghead,rval,sval,obsbin) endif !lagptr => lagptr%llpoint - lagptr => lagNode_nextcast(lagptr) + lagptr => lagnode_nextcast(lagptr) end do if (l_do_adjoint) then lag_u_full(:,:,obsbin)=adu_tmp lag_v_full(:,:,obsbin)=adv_tmp - ! Scater back the UV adjoint values back to state vector - call lag_ADscatter_stateuv(ru,rv,obsbin) + ! Scatter back the uv adjoint values back to state vector + call lag_adscatter_stateuv(ru,rv,obsbin) lag_u_full(:,:,obsbin)=zero lag_v_full(:,:,obsbin)=zero diff --git a/src/gsi/intlcbas.f90 b/src/gsi/intlcbas.f90 index 7b0a695228..8bed5645e4 100644 --- a/src/gsi/intlcbas.f90 +++ b/src/gsi/intlcbas.f90 @@ -7,7 +7,7 @@ module intlcbasmod ! abstract: module for intlcbas ! ! program history log: -! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting +! 2016-05-18 guo - replaced ob_type with polymorphic obsnode through type casting ! ! subroutines included: ! sub intlcbas @@ -20,15 +20,15 @@ module intlcbasmod ! !$$$ end documentation block -use m_obsNode , only: obsNode -use m_lcbasNode, only: lcbasNode -use m_lcbasNode, only: lcbasNode_typecast -use m_lcbasNode, only: lcbasNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode , only: obsnode +use m_lcbasnode, only: lcbasnode +use m_lcbasnode, only: lcbasnode_typecast +use m_lcbasnode, only: lcbasnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intlcbas +private +public intlcbas contains @@ -68,7 +68,7 @@ subroutine intlcbas(lcbashead,rval,sval) implicit none ! Declare passed variables - class(obsNode),pointer,intent(in) :: lcbashead + class(obsnode),pointer,intent(in) :: lcbashead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -81,7 +81,7 @@ subroutine intlcbas(lcbashead,rval,sval) real(r_kind) cg_lcbas,p0,grad,wnotgross,wgross,pg_lcbas real(r_kind),pointer,dimension(:) :: slcbas real(r_kind),pointer,dimension(:) :: rlcbas - type(lcbasNode), pointer :: lcbasptr + type(lcbasnode), pointer :: lcbasptr ! If no lcbas data return if(.not. associated(lcbashead))return @@ -94,7 +94,7 @@ subroutine intlcbas(lcbashead,rval,sval) if(ier/=0)return !lcbasptr => lcbashead - lcbasptr => lcbasNode_typecast(lcbashead) + lcbasptr => lcbasnode_typecast(lcbashead) do while (associated(lcbasptr)) j1=lcbasptr%ij(1) j2=lcbasptr%ij(2) @@ -113,10 +113,10 @@ subroutine intlcbas(lcbashead,rval,sval) if (lsaveobsens) then grad = val*lcbasptr%raterr2*lcbasptr%err2 !-- lcbasptr%diags%obssen(jiter) = grad - call obsdiagNode_set(lcbasptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(lcbasptr%diags,jiter=jiter,obssen=grad) else !-- if (lcbasptr%luse) lcbasptr%diags%tldepart(jiter)=val - if (lcbasptr%luse) call obsdiagNode_set(lcbasptr%diags,jiter=jiter,tldepart=val) + if (lcbasptr%luse) call obsdiagnode_set(lcbasptr%diags,jiter=jiter,tldepart=val) endif endif @@ -146,7 +146,7 @@ subroutine intlcbas(lcbashead,rval,sval) endif !lcbasptr => lcbasptr%llpoint - lcbasptr => lcbasNode_nextcast(lcbasptr) + lcbasptr => lcbasnode_nextcast(lcbasptr) end do diff --git a/src/gsi/intlight.f90 b/src/gsi/intlight.f90 index b9e7e26244..838d39d678 100644 --- a/src/gsi/intlight.f90 +++ b/src/gsi/intlight.f90 @@ -2,19 +2,19 @@ module intlightmod !$$$ module documentation block ! . . . . -! module: intlightmod int module for the observation operator for lightning flash rate (LFR) +! module: intlightmod int module for the observation operator for lightning flash rate (lfr) ! ! prgmmr: k apodaca ! org: CSU/CIRA, Data Assimilation Group ! date: 2016-05-04 ! -! abstract: module for the tangent linear (flashrate_TL) and adjoint models (flashrate_AD) -! of LFR +! abstract: module for the tangent linear (flashrate_tl) and adjoint models (flashrate_ad) +! of lfr ! ! program history log: -! 2016-05-04 apodaca - implement TL and AD of the LFR observation operator -! 2018-02-08 apodaca - replaced ob_type with polymorphic obsNode through type casting -! 2019-03-01 j guo - encapsulated access to obsdiagNode through obsdiagNode_set() +! 2016-05-04 apodaca - implement tl and ad of the lfr observation operator +! 2018-02-08 apodaca - replaced ob_type with polymorphic obsnode through type casting +! 2019-03-01 j guo - encapsulated access to obsdiagnode through obsdiagnode_set() ! ! subroutines included: ! sub intlight_ @@ -26,18 +26,18 @@ module intlightmod ! machine: ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_lightNode, only: lightNode -use m_lightNode, only: lightNode_typecast -use m_lightNode, only: lightNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_lightnode, only: lightnode +use m_lightnode, only: lightnode_typecast +use m_lightnode, only: lightnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intlight +private +public intlight interface intlight; module procedure & - intlight_ + intlight_ end interface contains @@ -46,7 +46,7 @@ subroutine intlight_(lighthead,rval,sval) !$$$ subprogram documentation block ! . . . . -! subprogram: intlight TL and subsequent AD of the forward observation operator for LFR +! subprogram: intlight tl and subsequent ad of the forward observation operator for lfr ! prgmmr: k apodaca ! org: CSU/CIRA, Data Assimilation Group ! date: 2016-05-04 @@ -66,8 +66,8 @@ subroutine intlight_(lighthead,rval,sval) ! of the sesitivity (impact) of observations. ! ! program history log: -! 2018-01-18 k apodaca revision of AD code -! 2018-08-18 k apodaca add a the TL and AD of second oservation operator for lightning +! 2018-01-18 k apodaca revision of ad code +! 2018-08-18 k apodaca add a the tl and ad of second oservation operator for lightning ! observations suitable for non-hydrostatic, cloud-resolving models ! with additional ice-phase hydrometeor control variables ! @@ -107,7 +107,7 @@ subroutine intlight_(lighthead,rval,sval) implicit none ! Declare passed variables - class(obsNode), pointer, intent(in ) :: lighthead + class(obsnode), pointer, intent(in ) :: lighthead type(gsi_bundle), intent(in ) :: sval type(gsi_bundle), intent(inout) :: rval @@ -118,32 +118,32 @@ subroutine intlight_(lighthead,rval,sval) real(r_kind) cg_light,grad,p0,wnotgross,wgross,pg_light real(r_kind),pointer,dimension(:) :: sq,sqi,sqs,sqg,su,sv,st real(r_kind),pointer,dimension(:) :: rq,rqi,rqs,rqg,ru,rv,rt - type(lightNode), pointer :: lightptr - -! Variables for TL and AD of lightning flash rate - real(r_kind),dimension(1:nsig) :: z_TL - real(r_kind),dimension(1:nsig) :: horiz_adv_TL - real(r_kind),dimension(1:nsig) :: vert_adv_TL - real(r_kind),dimension(1:nsig) :: w_TL - real(r_kind) :: wmaxi1_TL,wmaxi2_TL,wmaxi3_TL,wmaxi4_TL - real(r_kind) :: flashrate_TL,flashratei1_TL,flashratei2_TL - real(r_kind) :: flashratei3_TL, flashratei4_TL - real(r_kind) :: h1i1_TL,h1i2_TL,h1i3_TL,h1i4_TL - real(r_kind) :: h2i1_TL,h2i2_TL,h2i3_TL,h2i4_TL - real(r_kind) :: totice_colinti1_TL,totice_colinti2_TL - real(r_kind) :: totice_colinti3_TL,totice_colinti4_TL - real(r_kind) :: htot_TL,htoti1_TL,htoti2_TL,htoti3_TL,htoti4_TL - real(r_kind) :: flashrate_AD,flashratei1_AD,flashratei2_AD - real(r_kind) :: flashratei3_AD,flashratei4_AD - real(r_kind) :: wmaxi1_AD,wmaxi2_AD,wmaxi3_AD,wmaxi4_AD - real(r_kind) :: h1i1_AD,h1i2_AD,h1i3_AD,h1i4_AD - real(r_kind) :: h2i1_AD,h2i2_AD,h2i3_AD,h2i4_AD - real(r_kind) :: totice_colinti1_AD,totice_colinti2_AD - real(r_kind) :: totice_colinti3_AD,totice_colinti4_AD - real(r_kind) :: htot_AD,htoti1_AD,htoti2_AD,htoti3_AD,htoti4_AD - real(r_kind),dimension(1:nsig) :: z_AD - real(r_kind),dimension(1:nsig) :: w_AD - real(r_kind),dimension(1:nsig) :: vert_adv_AD,horiz_adv_AD + type(lightnode), pointer :: lightptr + +! Variables for tl and ad of lightning flash rate + real(r_kind),dimension(1:nsig) :: z_tl + real(r_kind),dimension(1:nsig) :: horiz_adv_tl + real(r_kind),dimension(1:nsig) :: vert_adv_tl + real(r_kind),dimension(1:nsig) :: w_tl + real(r_kind) :: wmaxi1_tl,wmaxi2_tl,wmaxi3_tl,wmaxi4_tl + real(r_kind) :: flashrate_tl,flashratei1_tl,flashratei2_tl + real(r_kind) :: flashratei3_tl, flashratei4_tl + real(r_kind) :: h1i1_tl,h1i2_tl,h1i3_tl,h1i4_tl + real(r_kind) :: h2i1_tl,h2i2_tl,h2i3_tl,h2i4_tl + real(r_kind) :: totice_colinti1_tl,totice_colinti2_tl + real(r_kind) :: totice_colinti3_tl,totice_colinti4_tl + real(r_kind) :: htot_tl,htoti1_tl,htoti2_tl,htoti3_tl,htoti4_tl + real(r_kind) :: flashrate_ad,flashratei1_ad,flashratei2_ad + real(r_kind) :: flashratei3_ad,flashratei4_ad + real(r_kind) :: wmaxi1_ad,wmaxi2_ad,wmaxi3_ad,wmaxi4_ad + real(r_kind) :: h1i1_ad,h1i2_ad,h1i3_ad,h1i4_ad + real(r_kind) :: h2i1_ad,h2i2_ad,h2i3_ad,h2i4_ad + real(r_kind) :: totice_colinti1_ad,totice_colinti2_ad + real(r_kind) :: totice_colinti3_ad,totice_colinti4_ad + real(r_kind) :: htot_ad,htoti1_ad,htoti2_ad,htoti3_ad,htoti4_ad + real(r_kind),dimension(1:nsig) :: z_ad + real(r_kind),dimension(1:nsig) :: w_ad + real(r_kind),dimension(1:nsig) :: vert_adv_ad,horiz_adv_ad real(r_kind),dimension(1:nsig) :: diffq real(r_kind),dimension(1:nsig) :: difft real(r_kind),dimension(1:nsig) :: diffz @@ -178,7 +178,7 @@ subroutine intlight_(lighthead,rval,sval) - lightptr => lightNode_typecast(lighthead) + lightptr => lightnode_typecast(lighthead) do while (associated(lightptr)) ! Load location information into local variables @@ -205,9 +205,9 @@ subroutine intlight_(lighthead,rval,sval) ! . . . . -! In the case of lightning observations (e.g. GOES/GLM), the schematic shown below is +! In the case of lightning observations (e.g. goes/glm), the schematic shown below is ! used for bi-linear interpolation of background fields to the location of an observation -! (+) and for the finite-difference derivation method used in the calculation of the TL of +! (+) and for the finite-difference derivation method used in the calculation of the tl of ! the observation operator for lightning flash rate. Calculations are done ! at each quadrant, i.e., central, north, south, east, and west. ! @@ -230,220 +230,220 @@ subroutine intlight_(lighthead,rval,sval) ! Tangent linear of height (z) - z_TL(:)=zero - horiz_adv_TL(:)=zero + z_tl(:)=zero + horiz_adv_tl(:)=zero do k=2,nsig-1 - z_TL(i1(1))=lightptr%jac_z0i1 - z_TL(i2(1))=lightptr%jac_z0i2 - z_TL(i3(1))=lightptr%jac_z0i3 - z_TL(i4(1))=lightptr%jac_z0i4 - z_TL(i5(1))=lightptr%jac_z0i5 - z_TL(i6(1))=lightptr%jac_z0i6 - z_TL(i7(1))=lightptr%jac_z0i7 - z_TL(i8(1))=lightptr%jac_z0i8 - z_TL(i9(1))=lightptr%jac_z0i9 - z_TL(i10(1))=lightptr%jac_z0i10 - z_TL(i11(1))=lightptr%jac_z0i11 - z_TL(i12(1))=lightptr%jac_z0i12 - - - z_TL(i1(k))=z_TL(i1(k-1))+lightptr%jac_vertti1(k)*st(i1(k)) & + z_tl(i1(1))=lightptr%jac_z0i1 + z_tl(i2(1))=lightptr%jac_z0i2 + z_tl(i3(1))=lightptr%jac_z0i3 + z_tl(i4(1))=lightptr%jac_z0i4 + z_tl(i5(1))=lightptr%jac_z0i5 + z_tl(i6(1))=lightptr%jac_z0i6 + z_tl(i7(1))=lightptr%jac_z0i7 + z_tl(i8(1))=lightptr%jac_z0i8 + z_tl(i9(1))=lightptr%jac_z0i9 + z_tl(i10(1))=lightptr%jac_z0i10 + z_tl(i11(1))=lightptr%jac_z0i11 + z_tl(i12(1))=lightptr%jac_z0i12 + + + z_tl(i1(k))=z_tl(i1(k-1))+lightptr%jac_vertti1(k)*st(i1(k)) & +lightptr%jac_vertqi1(k)*sq(i1(k)) - z_TL(i2(k))=z_TL(i2(k-1))+lightptr%jac_vertti2(k)*st(i2(k)) & + z_tl(i2(k))=z_tl(i2(k-1))+lightptr%jac_vertti2(k)*st(i2(k)) & +lightptr%jac_vertqi2(k)*sq(i2(k)) - z_TL(i3(k))=z_TL(i3(k-1))+lightptr%jac_vertti3(k)*st(i3(k)) & + z_tl(i3(k))=z_tl(i3(k-1))+lightptr%jac_vertti3(k)*st(i3(k)) & +lightptr%jac_vertqi3(k)*sq(i3(k)) - z_TL(i4(k))=z_TL(i4(k-1))+lightptr%jac_vertti4(k)*st(i4(k)) & + z_tl(i4(k))=z_tl(i4(k-1))+lightptr%jac_vertti4(k)*st(i4(k)) & +lightptr%jac_vertqi4(k)*sq(i4(k)) - z_TL(i5(k))=z_TL(i5(k-1))+lightptr%jac_vertti5(k)*st(i5(k)) & + z_tl(i5(k))=z_tl(i5(k-1))+lightptr%jac_vertti5(k)*st(i5(k)) & +lightptr%jac_vertqi5(k)*sq(i5(k)) - z_TL(i6(k))=z_TL(i6(k-1))+lightptr%jac_vertti6(k)*st(i6(k)) & + z_tl(i6(k))=z_tl(i6(k-1))+lightptr%jac_vertti6(k)*st(i6(k)) & +lightptr%jac_vertqi6(k)*sq(i6(k)) - z_TL(i7(k))=z_TL(i7(k-1))+lightptr%jac_vertti7(k)*st(i7(k)) & + z_tl(i7(k))=z_tl(i7(k-1))+lightptr%jac_vertti7(k)*st(i7(k)) & +lightptr%jac_vertqi7(k)*sq(i7(k)) - z_TL(i8(k))=z_TL(i8(k-1))+lightptr%jac_vertti8(k)*st(i8(k)) & + z_tl(i8(k))=z_tl(i8(k-1))+lightptr%jac_vertti8(k)*st(i8(k)) & +lightptr%jac_vertqi8(k)*sq(i8(k)) - z_TL(i9(k))=z_TL(i9(k-1))+lightptr%jac_vertti9(k)*st(i9(k)) & + z_tl(i9(k))=z_tl(i9(k-1))+lightptr%jac_vertti9(k)*st(i9(k)) & +lightptr%jac_vertqi9(k)*sq(i9(k)) - z_TL(i10(k))=z_TL(i10(k-1))+lightptr%jac_vertti10(k)*st(i10(k)) & + z_tl(i10(k))=z_tl(i10(k-1))+lightptr%jac_vertti10(k)*st(i10(k)) & +lightptr%jac_vertqi10(k)*sq(i10(k)) - z_TL(i11(k))=z_TL(i11(k-1))+lightptr%jac_vertti11(k)*st(i11(k)) & + z_tl(i11(k))=z_tl(i11(k-1))+lightptr%jac_vertti11(k)*st(i11(k)) & +lightptr%jac_vertqi11(k)*sq(i11(k)) - z_TL(i12(k))=z_TL(i12(k-1))+lightptr%jac_vertti12(k)*st(i12(k)) & + z_tl(i12(k))=z_tl(i12(k-1))+lightptr%jac_vertti12(k)*st(i12(k)) & +lightptr%jac_vertqi12(k)*sq(i12(k)) -! Tangent Linear of the Horizontal Advection Section +! Tangent linear of the horizontal advection section - horiz_adv_TL(i1(k))=lightptr%jac_zdxi1(k)*su(i1(k)) & + horiz_adv_tl(i1(k))=lightptr%jac_zdxi1(k)*su(i1(k)) & +lightptr%jac_zdyi1(k)*sv(i1(k)) & - +lightptr%jac_udxi1(k)*(z_TL(i3(k))-z_TL(i9(k))) & - +lightptr%jac_vdyi1(k)*(z_TL(i2(k))-z_TL(i5(k))) - horiz_adv_TL(i2(k))=lightptr%jac_zdxi2(k)*su(i2(k)) & + +lightptr%jac_udxi1(k)*(z_tl(i3(k))-z_tl(i9(k))) & + +lightptr%jac_vdyi1(k)*(z_tl(i2(k))-z_tl(i5(k))) + horiz_adv_tl(i2(k))=lightptr%jac_zdxi2(k)*su(i2(k)) & +lightptr%jac_zdyi2(k)*sv(i2(k)) & - +lightptr%jac_udxi2(k)*(z_TL(i4(k))-z_TL(i10(k))) & - +lightptr%jac_vdyi2(k)*(z_TL(i6(k))-z_TL(i1 (k))) + +lightptr%jac_udxi2(k)*(z_tl(i4(k))-z_tl(i10(k))) & + +lightptr%jac_vdyi2(k)*(z_tl(i6(k))-z_tl(i1 (k))) - horiz_adv_TL(i3(k))=lightptr%jac_zdxi3(k)*su(i3(k)) & + horiz_adv_tl(i3(k))=lightptr%jac_zdxi3(k)*su(i3(k)) & +lightptr%jac_zdyi3(k)*sv(i3(k)) & - +lightptr%jac_udxi3(k)*(z_TL(i11(k))-z_TL(i1(k))) & - +lightptr%jac_vdyi3(k)*(z_TL(i4 (k))-z_TL(i7(k))) + +lightptr%jac_udxi3(k)*(z_tl(i11(k))-z_tl(i1(k))) & + +lightptr%jac_vdyi3(k)*(z_tl(i4 (k))-z_tl(i7(k))) - horiz_adv_TL(i4(k))=lightptr%jac_zdxi4(k)*su(i4(k)) & + horiz_adv_tl(i4(k))=lightptr%jac_zdxi4(k)*su(i4(k)) & +lightptr%jac_zdyi4(k)*sv(i4(k)) & - +lightptr%jac_udxi4(k)*(z_TL(i12(k))-z_TL(i2(k))) & - +lightptr%jac_vdyi4(k)*(z_TL(i8 (k))-z_TL(i3(k))) + +lightptr%jac_udxi4(k)*(z_tl(i12(k))-z_tl(i2(k))) & + +lightptr%jac_vdyi4(k)*(z_tl(i8 (k))-z_tl(i3(k))) enddo ! do k=2,nsig-1 -! Tangent Linear of the Vertical Advection Section +! Tangent linear of the vertical advection section - vert_adv_TL(:)=zero - w_TL(:)=zero + vert_adv_tl(:)=zero + w_tl(:)=zero do k=1,nsig-1 - vert_adv_TL(i1(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti1(k)* & + vert_adv_tl(i1(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti1(k)* & (((one+fv*lightptr%jac_qi1(k))*st(i1(k))) & +(lightptr%jac_ti1(k)*fv*sq(i1(k)))) - vert_adv_TL(i2(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti2(k)* & + vert_adv_tl(i2(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti2(k)* & (((one+fv*lightptr%jac_qi2(k))*st(i2(k))) & +(lightptr%jac_ti2(k)*fv*sq(i2(k)))) - vert_adv_TL(i3(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti3(k)* & + vert_adv_tl(i3(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti3(k)* & (((one+fv*lightptr%jac_qi3(k))*st(i3(k))) & +(lightptr%jac_ti3(k)*fv*sq(i3(k)))) - vert_adv_TL(i4(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti4(k)* & + vert_adv_tl(i4(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti4(k)* & (((one+fv*lightptr%jac_qi4(k))*st(i4(k))) & +(lightptr%jac_ti4(k)*fv*sq(i4(k)))) -! Tangent Linear of Vertical Velocity +! Tangent linear of vertical velocity - w_TL(i1(k))=horiz_adv_TL(i1(k))+vert_adv_TL(i1(k)) - w_TL(i2(k))=horiz_adv_TL(i2(k))+vert_adv_TL(i2(k)) - w_TL(i3(k))=horiz_adv_TL(i3(k))+vert_adv_TL(i3(k)) - w_TL(i4(k))=horiz_adv_TL(i4(k))+vert_adv_TL(i4(k)) + w_tl(i1(k))=horiz_adv_tl(i1(k))+vert_adv_tl(i1(k)) + w_tl(i2(k))=horiz_adv_tl(i2(k))+vert_adv_tl(i2(k)) + w_tl(i3(k))=horiz_adv_tl(i3(k))+vert_adv_tl(i3(k)) + w_tl(i4(k))=horiz_adv_tl(i4(k))+vert_adv_tl(i4(k)) enddo !do k=1,nsig-1 ! . . . . -! Tangent Linear of lightning flash rate +! Tangent linear of lightning flash rate ! . . . . ! Regional if (regional) then -!-- WRF-ARW +!-- wrf-arw if (wrf_mass_regional) then -! Tangent linear - Lightning flash rate as a function of +! Tangent linear - lightning flash rate as a function of ! vertical graupel flux within the mixed-phase region ! (-15 deg C) if (lightptr%kboti1 > zero) then - h1i1_TL=lightptr%jac_qgmai1(lightptr%kboti1)*sqg(i1(lightptr%kboti1))+& - lightptr%jac_qgmbi1(lightptr%kboti1)*& - (half*(w_TL(i1(lightptr%kboti1))+w_TL(i1(lightptr%kboti1+1)))) - h1i1_TL=h1i1_TL/(abs(h1i1_TL)) + h1i1_tl=lightptr%jac_qgmai1(lightptr%kboti1)*sqg(i1(lightptr%kboti1))+& + lightptr%jac_qgmbi1(lightptr%kboti1)*& + (half*(w_tl(i1(lightptr%kboti1))+w_tl(i1(lightptr%kboti1+1)))) + h1i1_tl=h1i1_tl/(abs(h1i1_tl)) else - h1i1_TL=zero + h1i1_tl=zero endif if (lightptr%kboti2 > zero) then - h1i2_TL=lightptr%jac_qgmai2(lightptr%kboti2)*sqg(i2(lightptr%kboti2))+& - lightptr%jac_qgmbi2(lightptr%kboti2)*& - (half*(w_TL(i2(lightptr%kboti2))+w_TL(i2(lightptr%kboti2+1)))) - h1i2_TL=h1i2_TL/(abs(h1i2_TL)) + h1i2_tl=lightptr%jac_qgmai2(lightptr%kboti2)*sqg(i2(lightptr%kboti2))+& + lightptr%jac_qgmbi2(lightptr%kboti2)*& + (half*(w_tl(i2(lightptr%kboti2))+w_tl(i2(lightptr%kboti2+1)))) + h1i2_tl=h1i2_tl/(abs(h1i2_tl)) else - h1i2_TL=zero + h1i2_tl=zero endif if (lightptr%kboti3 > zero) then - h1i3_TL=lightptr%jac_qgmai3(lightptr%kboti3)*sqg(i3(lightptr%kboti3))+& - lightptr%jac_qgmbi3(lightptr%kboti3)*& - (half*(w_TL(i3(lightptr%kboti3))+w_TL(i3(lightptr%kboti3+1)))) - h1i3_TL=h1i3_TL/(abs(h1i3_TL)) + h1i3_tl=lightptr%jac_qgmai3(lightptr%kboti3)*sqg(i3(lightptr%kboti3))+& + lightptr%jac_qgmbi3(lightptr%kboti3)*& + (half*(w_tl(i3(lightptr%kboti3))+w_tl(i3(lightptr%kboti3+1)))) + h1i3_tl=h1i3_tl/(abs(h1i3_tl)) else - h1i3_TL=zero + h1i3_tl=zero endif if (lightptr%kboti4 > zero) then - h1i4_TL=lightptr%jac_qgmai4(lightptr%kboti4)*sqg(i4(lightptr%kboti4))+& - lightptr%jac_qgmbi4(lightptr%kboti4)*& - (half*(w_TL(i4(lightptr%kboti4))+w_TL(i4(lightptr%kboti4+1)))) - h1i4_TL=h1i4_TL/(abs(h1i4_TL)) + h1i4_tl=lightptr%jac_qgmai4(lightptr%kboti4)*sqg(i4(lightptr%kboti4))+& + lightptr%jac_qgmbi4(lightptr%kboti4)*& + (half*(w_tl(i4(lightptr%kboti4))+w_tl(i4(lightptr%kboti4+1)))) + h1i4_tl=h1i4_tl/(abs(h1i4_tl)) else - h1i4_TL=zero + h1i4_tl=zero endif -! Tangent Linear - Lightning flash rate as a function of total column-integrated +! Tangent linear - lightning flash rate as a function of total column-integrated ! ice-phase hydrometeors - totice_colinti1_TL=zero - totice_colinti2_TL=zero - totice_colinti3_TL=zero - totice_colinti4_TL=zero + totice_colinti1_tl=zero + totice_colinti2_tl=zero + totice_colinti3_tl=zero + totice_colinti4_tl=zero do k=1,nsig-1 - totice_colinti1_TL = totice_colinti1_TL+lightptr%jac_icei1(k) * & - (sqi(i1(k))+sqs(i1(k))+sqg(i1(k)))+& - lightptr%jac_zicei1(k)*z_TL(i1(k)) + totice_colinti1_tl = totice_colinti1_tl+lightptr%jac_icei1(k) * & + (sqi(i1(k))+sqs(i1(k))+sqg(i1(k)))+& + lightptr%jac_zicei1(k)*z_tl(i1(k)) - totice_colinti2_TL = totice_colinti2_TL+lightptr%jac_icei2(k) * & - (sqi(i2(k))+sqs(i2(k))+sqg(i2(k)))+& - lightptr%jac_zicei2(k)*z_TL(i2(k)) + totice_colinti2_tl = totice_colinti2_tl+lightptr%jac_icei2(k) * & + (sqi(i2(k))+sqs(i2(k))+sqg(i2(k)))+& + lightptr%jac_zicei2(k)*z_tl(i2(k)) - totice_colinti3_TL = totice_colinti3_TL+lightptr%jac_icei3(k) * & - (sqi(i3(k))+sqs(i3(k))+sqg(i3(k)))+& - lightptr%jac_zicei3(k)*z_TL(i3(k)) + totice_colinti3_tl = totice_colinti3_tl+lightptr%jac_icei3(k) * & + (sqi(i3(k))+sqs(i3(k))+sqg(i3(k)))+& + lightptr%jac_zicei3(k)*z_tl(i3(k)) - totice_colinti4_TL = totice_colinti4_TL+lightptr%jac_icei4(k) * & - (sqi(i4(k))+sqs(i4(k))+sqg(i4(k)))+& - lightptr%jac_zicei4(k)*z_TL(i4(k)) + totice_colinti4_tl = totice_colinti4_tl+lightptr%jac_icei4(k) * & + (sqi(i4(k))+sqs(i4(k))+sqg(i4(k)))+& + lightptr%jac_zicei4(k)*z_tl(i4(k)) enddo !do k=1,nsig-1 - h2i1_TL=(1-k3)*totice_colinti1_TL - h2i2_TL=(1-k3)*totice_colinti2_TL - h2i3_TL=(1-k3)*totice_colinti3_TL - h2i4_TL=(1-k3)*totice_colinti4_TL + h2i1_tl=(1-k3)*totice_colinti1_tl + h2i2_tl=(1-k3)*totice_colinti2_tl + h2i3_tl=(1-k3)*totice_colinti3_tl + h2i4_tl=(1-k3)*totice_colinti4_tl - htoti1_TL= h1i1_TL+h2i1_TL - htoti2_TL= h1i2_TL+h2i2_TL - htoti3_TL= h1i3_TL+h2i3_TL - htoti4_TL= h1i4_TL+h2i4_TL + htoti1_tl= h1i1_tl+h2i1_tl + htoti2_tl= h1i2_tl+h2i2_tl + htoti3_tl= h1i3_tl+h2i3_tl + htoti4_tl= h1i4_tl+h2i4_tl -! Interpolation of lightning flash rate to observation location (2D field) -! Forward Model +! Interpolation of lightning flash rate to observation location (2d field) +! forward model - htot_TL = (w1*htoti1_TL + w2*htoti2_TL + & - w3*htoti3_TL + w4*htoti4_TL) - val = htot_TL + htot_tl = (w1*htoti1_tl + w2*htoti2_tl + & + w3*htoti3_tl + w4*htoti4_tl) + val = htot_tl endif ! wrf_mass_regional @@ -454,95 +454,95 @@ subroutine intlight_(lighthead,rval,sval) if (.not. regional) then ! Global -! Cloud Mask +! Cloud mask ! If clouds are present, find the maximum value of vertical velocity -! (wmax_TL) at four points sorounding an observation (+) -! and amongst all vertical levels, otherwise set wmax_TL to zero. +! (wmax_tl) at four points sorounding an observation (+) +! and amongst all vertical levels, otherwise set wmax_tl to zero. - wmaxi1_TL=zero - wmaxi2_TL=zero - wmaxi3_TL=zero - wmaxi4_TL=zero + wmaxi1_tl=zero + wmaxi2_tl=zero + wmaxi3_tl=zero + wmaxi4_tl=zero if (lightptr%jac_wmaxflagi1) then - wmax=-1.e+10_r_kind - do k=1,nsig-1 - if (w_TL(i1(k)) > wmax) then - lightptr%jac_kverti1=k - wmaxi1_TL=w_TL(i1(lightptr%jac_kverti1)) - endif - if (wmaxi1_TL < zero) then - wmaxi1_TL=zero - endif - enddo ! k loop + wmax=-1.e+10_r_kind + do k=1,nsig-1 + if (w_tl(i1(k)) > wmax) then + lightptr%jac_kverti1=k + wmaxi1_tl=w_tl(i1(lightptr%jac_kverti1)) + endif + if (wmaxi1_tl < zero) then + wmaxi1_tl=zero + endif + enddo ! k loop endif if (lightptr%jac_wmaxflagi2) then wmax=-1.e+10_r_kind do k=1,nsig-1 - if (w_TL(i2(k)) > wmax) then + if (w_tl(i2(k)) > wmax) then lightptr%jac_kverti2=k - wmaxi2_TL=w_TL(i2(lightptr%jac_kverti2)) + wmaxi2_tl=w_tl(i2(lightptr%jac_kverti2)) endif - if (wmaxi2_TL < zero) then - wmaxi2_TL=zero + if (wmaxi2_tl < zero) then + wmaxi2_tl=zero endif enddo ! k loop endif if (lightptr%jac_wmaxflagi3) then - wmax=-1.e+10_r_kind - do k=1,nsig-1 - if (w_TL(i3(k)) > wmax) then - lightptr%jac_kverti3=k - wmaxi3_TL=w_TL(i3(lightptr%jac_kverti3)) - endif - if (wmaxi3_TL < zero) then - wmaxi3_TL=zero - endif - enddo ! k loop + wmax=-1.e+10_r_kind + do k=1,nsig-1 + if (w_tl(i3(k)) > wmax) then + lightptr%jac_kverti3=k + wmaxi3_tl=w_tl(i3(lightptr%jac_kverti3)) + endif + if (wmaxi3_tl < zero) then + wmaxi3_tl=zero + endif + enddo ! k loop endif if (lightptr%jac_wmaxflagi4) then - wmax=-1.e+10_r_kind - do k=1,nsig-1 - if (w_TL(i4(k)) > wmax) then - lightptr%jac_kverti4=k - wmaxi4_TL=w_TL(i4(lightptr%jac_kverti4)) - endif - if (wmaxi4_TL < zero) then - wmaxi4_TL=zero - endif - enddo ! k loop + wmax=-1.e+10_r_kind + do k=1,nsig-1 + if (w_tl(i4(k)) > wmax) then + lightptr%jac_kverti4=k + wmaxi4_tl=w_tl(i4(lightptr%jac_kverti4)) + endif + if (wmaxi4_tl < zero) then + wmaxi4_tl=zero + endif + enddo ! k loop endif -! Tangent Linear of Lightning Flash Rate +! Tangent linear of lightning flash rate - flashratei1_TL=lightptr%jac_fratei1*wmaxi1_TL - flashratei2_TL=lightptr%jac_fratei1*wmaxi2_TL - flashratei3_TL=lightptr%jac_fratei1*wmaxi3_TL - flashratei4_TL=lightptr%jac_fratei1*wmaxi4_TL + flashratei1_tl=lightptr%jac_fratei1*wmaxi1_tl + flashratei2_tl=lightptr%jac_fratei1*wmaxi2_tl + flashratei3_tl=lightptr%jac_fratei1*wmaxi3_tl + flashratei4_tl=lightptr%jac_fratei1*wmaxi4_tl -! Interpolation of lightning flash rate to observation location (2D field) -! Forward Model +! Interpolation of lightning flash rate to observation location (2d field) +! forward model - flashrate_TL = (w1*flashratei1_TL + w2*flashratei2_TL + & - w3*flashratei3_TL + w4*flashratei4_TL) - val = flashrate_TL + flashrate_tl = (w1*flashratei1_tl + w2*flashratei2_tl + & + w3*flashratei3_tl + w4*flashratei4_tl) + val = flashrate_tl end if ! global block if (luse_obsdiag)then - if (lsaveobsens) then - grad = val*lightptr%raterr2*lightptr%err2 - !-- lightptr%diags%obssen(jiter) = grad - call obsdiagNode_set(lightptr%diags,jiter=jiter,obssen=grad) - else - !-- if (lightptr%luse) lightptr%diags%tldepart(jiter)=val - if (lightptr%luse) call obsdiagNode_set(lightptr%diags,jiter=jiter,tldepart=val) - endif - end if + if (lsaveobsens) then + grad = val*lightptr%raterr2*lightptr%err2 + !-- lightptr%diags%obssen(jiter) = grad + call obsdiagnode_set(lightptr%diags,jiter=jiter,obssen=grad) + else + !-- if (lightptr%luse) lightptr%diags%tldepart(jiter)=val + if (lightptr%luse) call obsdiagnode_set(lightptr%diags,jiter=jiter,tldepart=val) + endif + end if ! . . . . @@ -552,9 +552,9 @@ subroutine intlight_(lighthead,rval,sval) if (l_do_adjoint) then ! Difference from observation if (.not. lsaveobsens) then - if (.not. ladtest_obs) val=val-lightptr%res + if (.not. ladtest_obs) val=val-lightptr%res -! needed for gradient of nonlinear qc operator +! needed for gradient of nonlinear qc operator if (nlnqc_iter .and. lightptr%pg > tiny_r_kind .and. & lightptr%b > tiny_r_kind) then pg_light=lightptr%pg*varqc_iter @@ -570,413 +570,413 @@ subroutine intlight_(lighthead,rval,sval) else grad = val*lightptr%raterr2*lightptr%err2 end if - endif + endif ! . . . . -! Adjoint of the Lightning Flash Rate Observation Operator +! Adjoint of the lightning flash rate observation operator ! . . . . ! Variable initialization - z_AD(:)=zero - w_AD(:)=zero + z_ad(:)=zero + w_ad(:)=zero ! Regional - if (regional) then + if (regional) then -!-- WRF-ARW +!-- wrf-arw - if (wrf_mass_regional) then + if (wrf_mass_regional) then - htot_AD=grad + htot_ad=grad -! Adjoint - Total lightning flash rate +! Adjoint - total lightning flash rate - htoti1_AD=htoti1_AD+w1*htot_AD - htoti2_AD=htoti2_AD+w1*htot_AD - htoti3_AD=htoti3_AD+w1*htot_AD - htoti4_AD=htoti4_AD+w1*htot_AD + htoti1_ad=htoti1_ad+w1*htot_ad + htoti2_ad=htoti2_ad+w1*htot_ad + htoti3_ad=htoti3_ad+w1*htot_ad + htoti4_ad=htoti4_ad+w1*htot_ad - h1i1_AD=h1i1_AD+htoti1_AD - h2i1_AD=h2i1_AD+htoti1_AD + h1i1_ad=h1i1_ad+htoti1_ad + h2i1_ad=h2i1_ad+htoti1_ad - h1i2_AD=h1i2_AD+htoti2_AD - h2i2_AD=h2i2_AD+htoti2_AD + h1i2_ad=h1i2_ad+htoti2_ad + h2i2_ad=h2i2_ad+htoti2_ad - h1i3_AD=h1i3_AD+htoti3_AD - h2i3_AD=h2i3_AD+htoti3_AD + h1i3_ad=h1i3_ad+htoti3_ad + h2i3_ad=h2i3_ad+htoti3_ad - h1i4_AD=h1i4_AD+htoti4_AD - h2i4_AD=h2i4_AD+htoti4_AD + h1i4_ad=h1i4_ad+htoti4_ad + h2i4_ad=h2i4_ad+htoti4_ad - totice_colinti1_AD=totice_colinti1_AD+(1-k3)*h2i1_AD - totice_colinti2_AD=totice_colinti2_AD+(1-k3)*h2i2_AD - totice_colinti3_AD=totice_colinti3_AD+(1-k3)*h2i3_AD - totice_colinti4_AD=totice_colinti4_AD+(1-k3)*h2i4_AD + totice_colinti1_ad=totice_colinti1_ad+(1-k3)*h2i1_ad + totice_colinti2_ad=totice_colinti2_ad+(1-k3)*h2i2_ad + totice_colinti3_ad=totice_colinti3_ad+(1-k3)*h2i3_ad + totice_colinti4_ad=totice_colinti4_ad+(1-k3)*h2i4_ad -! Adjoint - Lightning flash rate as a function of total column-integrated +! Adjoint - lightning flash rate as a function of total column-integrated ! ice-phase hydrometeors - do k=nsig-1,1,-1 + do k=nsig-1,1,-1 - z_AD(i1(k))=z_AD(i1(k))+lightptr%jac_zicei1(k)*totice_colinti1_AD - rqi(i1(k))=rqi(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_AD - rqs(i1(k))=rqs(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_AD - rqg(i1(k))=rqg(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_AD - totice_colinti1_AD=two*totice_colinti1_AD + z_ad(i1(k))=z_ad(i1(k))+lightptr%jac_zicei1(k)*totice_colinti1_ad + rqi(i1(k))=rqi(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_ad + rqs(i1(k))=rqs(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_ad + rqg(i1(k))=rqg(i1(k))+lightptr%jac_icei1(k)*totice_colinti1_ad + totice_colinti1_ad=two*totice_colinti1_ad - z_AD(i2(k))=z_AD(i2(k))+lightptr%jac_zicei2(k)*totice_colinti2_AD - rqi(i2(k))=rqi(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_AD - rqs(i2(k))=rqs(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_AD - rqg(i2(k))=rqg(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_AD - totice_colinti2_AD=two*totice_colinti2_AD + z_ad(i2(k))=z_ad(i2(k))+lightptr%jac_zicei2(k)*totice_colinti2_ad + rqi(i2(k))=rqi(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_ad + rqs(i2(k))=rqs(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_ad + rqg(i2(k))=rqg(i2(k))+lightptr%jac_icei2(k)*totice_colinti2_ad + totice_colinti2_ad=two*totice_colinti2_ad - z_AD(i3(k))=z_AD(i3(k))+lightptr%jac_zicei3(k)*totice_colinti3_AD - rqi(i3(k))=rqi(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_AD - rqs(i3(k))=rqs(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_AD - rqg(i3(k))=rqg(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_AD - totice_colinti3_AD=two*totice_colinti3_AD + z_ad(i3(k))=z_ad(i3(k))+lightptr%jac_zicei3(k)*totice_colinti3_ad + rqi(i3(k))=rqi(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_ad + rqs(i3(k))=rqs(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_ad + rqg(i3(k))=rqg(i3(k))+lightptr%jac_icei3(k)*totice_colinti3_ad + totice_colinti3_ad=two*totice_colinti3_ad - z_AD(i4(k))=z_AD(i4(k))+lightptr%jac_zicei4(k)*totice_colinti4_AD - rqi(i4(k))=rqi(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_AD - rqs(i4(k))=rqs(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_AD - rqg(i4(k))=rqg(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_AD - totice_colinti4_AD=two*totice_colinti4_AD + z_ad(i4(k))=z_ad(i4(k))+lightptr%jac_zicei4(k)*totice_colinti4_ad + rqi(i4(k))=rqi(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_ad + rqs(i4(k))=rqs(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_ad + rqg(i4(k))=rqg(i4(k))+lightptr%jac_icei4(k)*totice_colinti4_ad + totice_colinti4_ad=two*totice_colinti4_ad -! Adjoint - Lightning flash rate as a function of +! Adjoint - lightning flash rate as a function of ! vertical graupel flux within the mixed-phase region -! (-15 deg C) - - if (lightptr%kboti1 > zero) then - h1i1_AD=h1i1_AD+(h1i1_AD/abs(h1i1_AD)) - rqg(i1(lightptr%kboti1))=rqg(i1(lightptr%kboti1))+& - lightptr%jac_qgmai1(lightptr%kboti1)*h1i1_AD - w_AD(i1(lightptr%kboti1))=w_AD(i1(lightptr%kboti1))+& - half*lightptr%jac_qgmbi1(lightptr%kboti1)*h1i1_AD - w_AD(i1(lightptr%kboti1+1))=w_AD(i1(lightptr%kboti1+1))+& - half*lightptr%jac_qgmbi1(lightptr%kboti1)*h1i1_AD - else - h1i1_AD=zero - rqg(i1(lightptr%kboti1))=zero - w_AD(i1(lightptr%kboti1))=zero - w_AD(i1(lightptr%kboti1+1))=zero - endif - - if (lightptr%kboti2 > zero) then - h1i2_AD=h1i2_AD+(h1i2_AD/abs(h1i2_AD)) - rqg(i2(lightptr%kboti2))=rqg(i2(lightptr%kboti2))+& - lightptr%jac_qgmai2(lightptr%kboti2)*h1i2_AD - w_AD(i2(lightptr%kboti2))=w_AD(i2(lightptr%kboti2))+& - half*lightptr%jac_qgmbi2(lightptr%kboti2)*h1i2_AD - w_AD(i2(lightptr%kboti2+1))=w_AD(i2(lightptr%kboti2+1))+& - half*lightptr%jac_qgmbi2(lightptr%kboti2)*h1i2_AD - else - h1i2_AD=zero - rqg(i2(lightptr%kboti2))=zero - w_AD(i2(lightptr%kboti2+1))=zero - endif - - if (lightptr%kboti3 > zero) then - h1i3_AD=h1i3_AD+(h1i3_AD/abs(h1i3_AD)) - rqg(i3(lightptr%kboti3))=rqg(i3(lightptr%kboti3))+& - lightptr%jac_qgmai3(lightptr%kboti3)*h1i3_AD - w_AD(i3(lightptr%kboti3))=w_AD(i3(lightptr%kboti3))+& - half*lightptr%jac_qgmbi3(lightptr%kboti3)*h1i3_AD - w_AD(i3(lightptr%kboti3+1))=w_AD(i3(lightptr%kboti3+1))+& - half*lightptr%jac_qgmbi3(lightptr%kboti3)*h1i3_AD - else - h1i3_AD=zero - rqg(i3(lightptr%kboti3))=zero - w_AD(i3(lightptr%kboti3+1))=zero - endif - - if (lightptr%kboti4 > zero) then - h1i4_AD=h1i4_AD+(h1i4_AD/abs(h1i4_AD)) - rqg(i4(lightptr%kboti4))=rqg(i4(lightptr%kboti4))+& - lightptr%jac_qgmai4(lightptr%kboti4)*h1i4_AD - w_AD(i4(lightptr%kboti4))=w_AD(i4(lightptr%kboti4))+& - half*lightptr%jac_qgmbi4(lightptr%kboti4)*h1i4_AD - w_AD(i4(lightptr%kboti4+1))=w_AD(i4(lightptr%kboti4+1))+& - half*lightptr%jac_qgmbi4(lightptr%kboti4)*h1i4_AD - else - h1i4_AD=zero - rqg(i4(lightptr%kboti4))=zero - w_AD(i4(lightptr%kboti4+1))=zero - endif - - - enddo !do k=nsig-1,1,-1 - - endif ! wrf_mass_regional - - endif !if (regional) then +! (-15 deg c) + + if (lightptr%kboti1 > zero) then + h1i1_ad=h1i1_ad+(h1i1_ad/abs(h1i1_ad)) + rqg(i1(lightptr%kboti1))=rqg(i1(lightptr%kboti1))+& + lightptr%jac_qgmai1(lightptr%kboti1)*h1i1_ad + w_ad(i1(lightptr%kboti1))=w_ad(i1(lightptr%kboti1))+& + half*lightptr%jac_qgmbi1(lightptr%kboti1)*h1i1_ad + w_ad(i1(lightptr%kboti1+1))=w_ad(i1(lightptr%kboti1+1))+& + half*lightptr%jac_qgmbi1(lightptr%kboti1)*h1i1_ad + else + h1i1_ad=zero + rqg(i1(lightptr%kboti1))=zero + w_ad(i1(lightptr%kboti1))=zero + w_ad(i1(lightptr%kboti1+1))=zero + endif + + if (lightptr%kboti2 > zero) then + h1i2_ad=h1i2_ad+(h1i2_ad/abs(h1i2_ad)) + rqg(i2(lightptr%kboti2))=rqg(i2(lightptr%kboti2))+& + lightptr%jac_qgmai2(lightptr%kboti2)*h1i2_ad + w_ad(i2(lightptr%kboti2))=w_ad(i2(lightptr%kboti2))+& + half*lightptr%jac_qgmbi2(lightptr%kboti2)*h1i2_ad + w_ad(i2(lightptr%kboti2+1))=w_ad(i2(lightptr%kboti2+1))+& + half*lightptr%jac_qgmbi2(lightptr%kboti2)*h1i2_ad + else + h1i2_ad=zero + rqg(i2(lightptr%kboti2))=zero + w_ad(i2(lightptr%kboti2+1))=zero + endif + + if (lightptr%kboti3 > zero) then + h1i3_ad=h1i3_ad+(h1i3_ad/abs(h1i3_ad)) + rqg(i3(lightptr%kboti3))=rqg(i3(lightptr%kboti3))+& + lightptr%jac_qgmai3(lightptr%kboti3)*h1i3_ad + w_ad(i3(lightptr%kboti3))=w_ad(i3(lightptr%kboti3))+& + half*lightptr%jac_qgmbi3(lightptr%kboti3)*h1i3_ad + w_ad(i3(lightptr%kboti3+1))=w_ad(i3(lightptr%kboti3+1))+& + half*lightptr%jac_qgmbi3(lightptr%kboti3)*h1i3_ad + else + h1i3_ad=zero + rqg(i3(lightptr%kboti3))=zero + w_ad(i3(lightptr%kboti3+1))=zero + endif + + if (lightptr%kboti4 > zero) then + h1i4_ad=h1i4_ad+(h1i4_ad/abs(h1i4_ad)) + rqg(i4(lightptr%kboti4))=rqg(i4(lightptr%kboti4))+& + lightptr%jac_qgmai4(lightptr%kboti4)*h1i4_ad + w_ad(i4(lightptr%kboti4))=w_ad(i4(lightptr%kboti4))+& + half*lightptr%jac_qgmbi4(lightptr%kboti4)*h1i4_ad + w_ad(i4(lightptr%kboti4+1))=w_ad(i4(lightptr%kboti4+1))+& + half*lightptr%jac_qgmbi4(lightptr%kboti4)*h1i4_ad + else + h1i4_ad=zero + rqg(i4(lightptr%kboti4))=zero + w_ad(i4(lightptr%kboti4+1))=zero + endif + + + enddo !do k=nsig-1,1,-1 + + endif ! wrf_mass_regional + + endif !if (regional) then ! . . . . ! Global - if (.not. regional) then + if (.not. regional) then - flashrate_AD=grad + flashrate_ad=grad - flashratei1_AD=flashratei1_AD+w1*flashrate_AD - flashratei2_AD=flashratei2_AD+w2*flashrate_AD - flashratei3_AD=flashratei3_AD+w3*flashrate_AD - flashratei4_AD=flashratei4_AD+w4*flashrate_AD - -! Adjoint of Maximum Vertical Velocity - - wmaxi1_AD=wmaxi1_AD+lightptr%jac_fratei1*flashratei1_AD - wmaxi2_AD=wmaxi2_AD+lightptr%jac_fratei2*flashratei2_AD - wmaxi3_AD=wmaxi3_AD+lightptr%jac_fratei3*flashratei3_AD - wmaxi4_AD=wmaxi3_AD+lightptr%jac_fratei4*flashratei4_AD - - if (lightptr%jac_wmaxflagi1) then - wmax=-1.e+10_r_kind - do k=nsig-1,1,-1 - if (wmaxi1_AD < zero) then - wmaxi1_AD=zero - endif - if (wmaxi1_AD > wmax) then - lightptr%jac_kverti1=k - w_AD(i1(lightptr%jac_kverti1))=w_AD(i1(lightptr%jac_kverti1))+wmaxi1_AD - endif - enddo - endif + flashratei1_ad=flashratei1_ad+w1*flashrate_ad + flashratei2_ad=flashratei2_ad+w2*flashrate_ad + flashratei3_ad=flashratei3_ad+w3*flashrate_ad + flashratei4_ad=flashratei4_ad+w4*flashrate_ad + +! Adjoint of maximum vertical velocity + + wmaxi1_ad=wmaxi1_ad+lightptr%jac_fratei1*flashratei1_ad + wmaxi2_ad=wmaxi2_ad+lightptr%jac_fratei2*flashratei2_ad + wmaxi3_ad=wmaxi3_ad+lightptr%jac_fratei3*flashratei3_ad + wmaxi4_ad=wmaxi3_ad+lightptr%jac_fratei4*flashratei4_ad + + if (lightptr%jac_wmaxflagi1) then + wmax=-1.e+10_r_kind + do k=nsig-1,1,-1 + if (wmaxi1_ad < zero) then + wmaxi1_ad=zero + endif + if (wmaxi1_ad > wmax) then + lightptr%jac_kverti1=k + w_ad(i1(lightptr%jac_kverti1))=w_ad(i1(lightptr%jac_kverti1))+wmaxi1_ad + endif + enddo + endif - if (lightptr%jac_wmaxflagi2) then - wmax=-1.e+10_r_kind - do k=nsig-1,1,-1 - if (wmaxi2_AD < zero) then - wmaxi2_AD=zero - endif - if (wmaxi2_AD > wmax) then - lightptr%jac_kverti2=k - w_AD(i2(lightptr%jac_kverti2))=w_AD(i2(lightptr%jac_kverti2))+wmaxi2_AD - endif - enddo - endif - - if (lightptr%jac_wmaxflagi3) then - wmax=-1.e+10_r_kind - do k=nsig-1,1,-1 - if (wmaxi3_AD < zero) then - wmaxi3_AD=zero - endif - if (wmaxi3_AD > wmax) then - lightptr%jac_kverti3=k - w_AD(i3(lightptr%jac_kverti3))=w_AD(i3(lightptr%jac_kverti3))+wmaxi3_AD - endif - enddo - endif - - if (lightptr%jac_wmaxflagi4) then - wmax=-1.e+10_r_kind - do k=nsig-1,1,-1 - if (wmaxi4_AD < zero) then - wmaxi4_AD=zero - endif - if (wmaxi4_AD > wmax) then - lightptr%jac_kverti4=k - w_AD(i4(lightptr%jac_kverti4))=w_AD(i4(lightptr%jac_kverti4))+wmaxi4_AD - endif - enddo - endif - - - endif ! end global block + if (lightptr%jac_wmaxflagi2) then + wmax=-1.e+10_r_kind + do k=nsig-1,1,-1 + if (wmaxi2_ad < zero) then + wmaxi2_ad=zero + endif + if (wmaxi2_ad > wmax) then + lightptr%jac_kverti2=k + w_ad(i2(lightptr%jac_kverti2))=w_ad(i2(lightptr%jac_kverti2))+wmaxi2_ad + endif + enddo + endif + + if (lightptr%jac_wmaxflagi3) then + wmax=-1.e+10_r_kind + do k=nsig-1,1,-1 + if (wmaxi3_ad < zero) then + wmaxi3_ad=zero + endif + if (wmaxi3_ad > wmax) then + lightptr%jac_kverti3=k + w_ad(i3(lightptr%jac_kverti3))=w_ad(i3(lightptr%jac_kverti3))+wmaxi3_ad + endif + enddo + endif + + if (lightptr%jac_wmaxflagi4) then + wmax=-1.e+10_r_kind + do k=nsig-1,1,-1 + if (wmaxi4_ad < zero) then + wmaxi4_ad=zero + endif + if (wmaxi4_ad > wmax) then + lightptr%jac_kverti4=k + w_ad(i4(lightptr%jac_kverti4))=w_ad(i4(lightptr%jac_kverti4))+wmaxi4_ad + endif + enddo + endif + + + endif ! end global block ! . . . . -! Adjoint of Vertical Velocity (from Vertical and Horizontal Advection) +! Adjoint of vertical velocity (from vertical and horizontal advection) - vert_adv_AD(:)=zero + vert_adv_ad(:)=zero - do k=nsig-1,1,-1 + do k=nsig-1,1,-1 - vert_adv_AD(i4(k))=vert_adv_AD(i4(k))+w_AD(i4(k)) - vert_adv_AD(i3(k))=vert_adv_AD(i3(k))+w_AD(i3(k)) - vert_adv_AD(i2(k))=vert_adv_AD(i2(k))+w_AD(i2(k)) - vert_adv_AD(i1(k))=vert_adv_AD(i1(k))+w_AD(i1(k)) + vert_adv_ad(i4(k))=vert_adv_ad(i4(k))+w_ad(i4(k)) + vert_adv_ad(i3(k))=vert_adv_ad(i3(k))+w_ad(i3(k)) + vert_adv_ad(i2(k))=vert_adv_ad(i2(k))+w_ad(i2(k)) + vert_adv_ad(i1(k))=vert_adv_ad(i1(k))+w_ad(i1(k)) - enddo + enddo - horiz_adv_AD(:)=zero + horiz_adv_ad(:)=zero - do k=nsig-1,2,-1 + do k=nsig-1,2,-1 - horiz_adv_AD(i4(k))=horiz_adv_AD(i4(k))+w_AD(i4(k)) - horiz_adv_AD(i4(k))=horiz_adv_AD(i3(k))+w_AD(i3(k)) - horiz_adv_AD(i2(k))=horiz_adv_AD(i2(k))+w_AD(i2(k)) - horiz_adv_AD(i1(k))=horiz_adv_AD(i1(k))+w_AD(i1(k)) + horiz_adv_ad(i4(k))=horiz_adv_ad(i4(k))+w_ad(i4(k)) + horiz_adv_ad(i4(k))=horiz_adv_ad(i3(k))+w_ad(i3(k)) + horiz_adv_ad(i2(k))=horiz_adv_ad(i2(k))+w_ad(i2(k)) + horiz_adv_ad(i1(k))=horiz_adv_ad(i1(k))+w_ad(i1(k)) - enddo + enddo -! Adjoint of q and t from the Vertical Advection Section +! Adjoint of q and t from the vertical advection section - diffq(:)=zero - difft(:)=zero + diffq(:)=zero + difft(:)=zero - do k=nsig-1,1,-1 + do k=nsig-1,1,-1 - diffq(i1(k))=-(lightptr%jac_ti1(k)*fv*lightptr%jac_vert(K) & - *lightptr%jac_sigdoti1(k))*vert_adv_AD(i1(k)) - difft(i1(k))=-((one+fv*lightptr%jac_qi1(k))*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti1(k))*vert_adv_AD(i1(k)) - diffq(i2(k))=-(lightptr%jac_ti2(k)*fv*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti2(k))*vert_adv_AD(i2(k)) - difft(i2(k))=-((one+fv*lightptr%jac_qi2(k))*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti2(k))*vert_adv_AD(i2(k)) - diffq(i3(k))=-(lightptr%jac_ti3(k)*fv*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti3(k))*vert_adv_AD(i3(k)) - difft(i3(k))=-((one+fv*lightptr%jac_qi3(k))*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti3(k))*vert_adv_AD(i3(k)) - diffq(i4(k))=-(lightptr%jac_ti4(k)*fv*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti4(k))*vert_adv_AD(i4(k)) - difft(i4(k))=-((one+fv*lightptr%jac_qi4(k))*lightptr%jac_vert(k) & - *lightptr%jac_sigdoti4(k))*vert_adv_AD(i4(k)) + diffq(i1(k))=-(lightptr%jac_ti1(k)*fv*lightptr%jac_vert(K) & + *lightptr%jac_sigdoti1(k))*vert_adv_ad(i1(k)) + difft(i1(k))=-((one+fv*lightptr%jac_qi1(k))*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti1(k))*vert_adv_ad(i1(k)) + diffq(i2(k))=-(lightptr%jac_ti2(k)*fv*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti2(k))*vert_adv_ad(i2(k)) + difft(i2(k))=-((one+fv*lightptr%jac_qi2(k))*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti2(k))*vert_adv_ad(i2(k)) + diffq(i3(k))=-(lightptr%jac_ti3(k)*fv*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti3(k))*vert_adv_ad(i3(k)) + difft(i3(k))=-((one+fv*lightptr%jac_qi3(k))*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti3(k))*vert_adv_ad(i3(k)) + diffq(i4(k))=-(lightptr%jac_ti4(k)*fv*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti4(k))*vert_adv_ad(i4(k)) + difft(i4(k))=-((one+fv*lightptr%jac_qi4(k))*lightptr%jac_vert(k) & + *lightptr%jac_sigdoti4(k))*vert_adv_ad(i4(k)) - rq(i1(k))=rq(i1(k))+diffq(i1(k)) + rq(i1(k))=rq(i1(k))+diffq(i1(k)) - rt(i1(k))=rt(i1(k))+difft(i1(k)) + rt(i1(k))=rt(i1(k))+difft(i1(k)) - rq(i2(k))=rq(i2(k))+diffq(i2(k)) + rq(i2(k))=rq(i2(k))+diffq(i2(k)) - rq(i3(k))=rq(i3(k))+diffq(i3(k)) + rq(i3(k))=rq(i3(k))+diffq(i3(k)) - rt(i3(k))=rt(i3(k))+difft(i3(k)) + rt(i3(k))=rt(i3(k))+difft(i3(k)) - rq(i4(k))=rq(i4(k))+diffq(i4(k)) + rq(i4(k))=rq(i4(k))+diffq(i4(k)) - rt(i4(k))=rt(i4(k))+difft(i4(k)) + rt(i4(k))=rt(i4(k))+difft(i4(k)) - enddo + enddo -! Adjoint of z, u, and v from the Horizontal Advection Section +! Adjoint of z, u, and v from the horizontal advection section - diffz(:)=zero - z_AD(:)=zero + diffz(:)=zero + z_ad(:)=zero - do k=nsig-1,2,-1 + do k=nsig-1,2,-1 - diffz(i5(k))=-lightptr%jac_vdyi1(k)*horiz_adv_AD(i1(k)) - diffz(i9(k))=-lightptr%jac_udxi1(k)*horiz_adv_AD(i1(k)) + diffz(i5(k))=-lightptr%jac_vdyi1(k)*horiz_adv_ad(i1(k)) + diffz(i9(k))=-lightptr%jac_udxi1(k)*horiz_adv_ad(i1(k)) - z_AD(i5(k))=z_AD(i5(k))+diffz(i5(k)) - z_AD(i2(k))=z_AD(i2(k))+(lightptr%jac_vdyi1(k)*horiz_adv_AD(i1(k))) - z_AD(i9(k))=z_AD(i9(k))+(diffz(i9(k))) - z_AD(i3(k))=z_AD(i3(k))+(lightptr%jac_udxi1(k)*horiz_adv_AD(i1(k))) + z_ad(i5(k))=z_ad(i5(k))+diffz(i5(k)) + z_ad(i2(k))=z_ad(i2(k))+(lightptr%jac_vdyi1(k)*horiz_adv_ad(i1(k))) + z_ad(i9(k))=z_ad(i9(k))+(diffz(i9(k))) + z_ad(i3(k))=z_ad(i3(k))+(lightptr%jac_udxi1(k)*horiz_adv_ad(i1(k))) - rv(i1(k))=rv(i1(k))+(lightptr%jac_zdyi1(k)*horiz_adv_AD(i1(k))) - ru(i1(k))=ru(i1(k))+(lightptr%jac_zdxi1(k)*horiz_adv_AD(i1(k))) - - diffz(i1(k)) =-lightptr%jac_vdyi2(k)*horiz_adv_AD(i2(k)) - diffz(i10(k))=-lightptr%jac_udxi2(k)*horiz_adv_AD(i2(k)) - - z_AD(i1(k))=z_AD(i1(k))+(diffz(i1(k))) - z_AD(i6(k))=z_AD(i6(k))+(lightptr%jac_vdyi2(k)*horiz_adv_AD(i2(k))) - z_AD(i10(k))=z_AD(i10(k))+(diffz(i10(k))) - z_AD(i4(k))=z_AD(i4(k))+(lightptr%jac_udxi2(k)*horiz_adv_AD(i2(k))) - rv(i2(k))=rv(i2(k))+(lightptr%jac_zdyi2(k)*horiz_adv_AD(i2(k))) - ru(i2(k))=ru(i2(k))+(lightptr%jac_zdxi2(k)*horiz_adv_AD(i2(k))) - - diffz(i7(k))= -lightptr%jac_vdyi3(k)*horiz_adv_AD(i3(k)) - diffz(i1(k))= -lightptr%jac_udxi3(k)*horiz_adv_AD(i3(k)) - - z_AD(i7(k)) = z_AD(i7(k))+diffz(i7(k)) - z_AD(i4(k)) = z_AD(i4(k))+(lightptr%jac_vdyi3(k)*horiz_adv_AD(i3(k))) - z_AD(i1(k)) = z_AD(i1(k))+diffz(i1(k)) - z_AD(i11(k))= z_AD(i11(k))+(lightptr%jac_udxi3(k)*horiz_adv_AD(i3(k))) - rv(i3(k)) = rv(i3(k))+(lightptr%jac_zdyi3(k)*horiz_adv_AD(i3(k))) - ru(i3(k)) = ru(i3(k))+(lightptr%jac_zdxi3(k)*horiz_adv_AD(i3(k))) - - diffz(i3(k))=-lightptr%jac_vdyi4(k)*horiz_adv_AD(i4(k)) - diffz(i2(k))=-z_TL(i2(k))-lightptr%jac_udxi4(k)*horiz_adv_AD(i4(k)) - - z_AD(i3(k)) = z_AD(i3(k))+diffz(i3(k)) - z_AD(i8(k)) = z_AD(i8(k))+(lightptr%jac_vdyi4(k)*horiz_adv_AD(i4(k))) - z_AD(i2(k)) = z_TL(i2(k))+diffz(i2(k)) - z_AD(i12(k))= z_AD(i12(k))+(lightptr%jac_udxi4(k)*horiz_adv_AD(i4(k))) - rv(i4(k)) = rv(i4(k))+(lightptr%jac_zdyi4(k)*horiz_adv_AD(i4(k))) - ru(i4(k)) = ru(i4(k))+(lightptr%jac_zdxi4(k)*horiz_adv_AD(i4(k))) - - enddo - -! Adjoint of q and t from the Calculation of Height (z) + rv(i1(k))=rv(i1(k))+(lightptr%jac_zdyi1(k)*horiz_adv_ad(i1(k))) + ru(i1(k))=ru(i1(k))+(lightptr%jac_zdxi1(k)*horiz_adv_ad(i1(k))) + + diffz(i1(k)) =-lightptr%jac_vdyi2(k)*horiz_adv_ad(i2(k)) + diffz(i10(k))=-lightptr%jac_udxi2(k)*horiz_adv_ad(i2(k)) + + z_ad(i1(k))=z_ad(i1(k))+(diffz(i1(k))) + z_ad(i6(k))=z_ad(i6(k))+(lightptr%jac_vdyi2(k)*horiz_adv_ad(i2(k))) + z_ad(i10(k))=z_ad(i10(k))+(diffz(i10(k))) + z_ad(i4(k))=z_ad(i4(k))+(lightptr%jac_udxi2(k)*horiz_adv_ad(i2(k))) + rv(i2(k))=rv(i2(k))+(lightptr%jac_zdyi2(k)*horiz_adv_ad(i2(k))) + ru(i2(k))=ru(i2(k))+(lightptr%jac_zdxi2(k)*horiz_adv_ad(i2(k))) + + diffz(i7(k))= -lightptr%jac_vdyi3(k)*horiz_adv_ad(i3(k)) + diffz(i1(k))= -lightptr%jac_udxi3(k)*horiz_adv_ad(i3(k)) + + z_ad(i7(k)) = z_ad(i7(k))+diffz(i7(k)) + z_ad(i4(k)) = z_ad(i4(k))+(lightptr%jac_vdyi3(k)*horiz_adv_ad(i3(k))) + z_ad(i1(k)) = z_ad(i1(k))+diffz(i1(k)) + z_ad(i11(k))= z_ad(i11(k))+(lightptr%jac_udxi3(k)*horiz_adv_ad(i3(k))) + rv(i3(k)) = rv(i3(k))+(lightptr%jac_zdyi3(k)*horiz_adv_ad(i3(k))) + ru(i3(k)) = ru(i3(k))+(lightptr%jac_zdxi3(k)*horiz_adv_ad(i3(k))) + + diffz(i3(k))=-lightptr%jac_vdyi4(k)*horiz_adv_ad(i4(k)) + diffz(i2(k))=-z_tl(i2(k))-lightptr%jac_udxi4(k)*horiz_adv_ad(i4(k)) + + z_ad(i3(k)) = z_ad(i3(k))+diffz(i3(k)) + z_ad(i8(k)) = z_ad(i8(k))+(lightptr%jac_vdyi4(k)*horiz_adv_ad(i4(k))) + z_ad(i2(k)) = z_tl(i2(k))+diffz(i2(k)) + z_ad(i12(k))= z_ad(i12(k))+(lightptr%jac_udxi4(k)*horiz_adv_ad(i4(k))) + rv(i4(k)) = rv(i4(k))+(lightptr%jac_zdyi4(k)*horiz_adv_ad(i4(k))) + ru(i4(k)) = ru(i4(k))+(lightptr%jac_zdxi4(k)*horiz_adv_ad(i4(k))) + + enddo + +! Adjoint of q and t from the calculation of height (z) - do k=nsig-1,2,-1 - - rq(i1(k))=rq(i1(k))+lightptr%jac_vertqi1(k)*z_AD(i1(k)) - rt(i1(k))=rt(i1(k))+lightptr%jac_vertti1(k)*z_AD(i1(k)) - z_AD(i1(k-1))=z_AD(i1(k-1))+z_AD(i1(k)) - z_AD(i1(k))=zero - - rq(i2(k))=rq(i2(k))+lightptr%jac_vertqi2(k)*z_AD(i2(k)) - rt(i2(k))=rt(i2(k))+lightptr%jac_vertti12(k)*z_AD(i2(k)) - z_AD(i2(k-1))=z_AD(i2(k-1))+z_AD(i2(k)) - z_AD(i2(k))=zero - - rq(i3(k))=rq(i3(k))+lightptr%jac_vertqi3(k)*z_AD(i3(k)) - rt(i3(k))=rt(i3(k))+lightptr%jac_vertti3(k)*z_AD(i3(k)) - z_AD(i3(k-1))=z_AD(i3(k-1))+z_AD(i3(k)) - z_AD(i3(k))=zero - - rq(i4(k))=rq(i4(k))+lightptr%jac_vertqi4(k)*z_AD(i4(k)) - rt(i4(k))=rt(i4(k))+lightptr%jac_vertti4(k)*z_AD(i4(k)) - z_AD(i4(k-1))=z_AD(i4(k-1))+z_AD(i4(k)) - z_AD(i4(k))=zero - - rq(i5(k))=rq(i5(k))+lightptr%jac_vertqi5(k)*z_AD(i5(k)) - rt(i5(k))=rt(i5(k))+lightptr%jac_vertti5(k)*z_AD(i5(k)) - z_AD(i5(k-1))=z_AD(i5(k-1))+z_AD(i5(k)) - z_AD(i5(k))=zero - - rq(i6(k))=rq(i6(k))+lightptr%jac_vertqi6(k)*z_AD(i6(k)) - rt(i6(k))=rt(i6(k))+lightptr%jac_vertti6(k)*z_AD(i6(k)) - z_AD(i6(k-1))=z_AD(i6(k-1))+z_AD(i6(k)) - z_AD(i6(k))=zero - - rq(i7(k))=rq(i7(k))+lightptr%jac_vertqi7(k)*z_AD(i7(k)) - rt(i7(k))=rt(i7(k))+lightptr%jac_vertti7(k)*z_AD(i7(k)) - z_AD(i7(k-1))=z_AD(i7(k-1))+z_AD(i7(k)) - z_AD(i7(k))=zero - - rq(i8(k))=rq(i8(k))+lightptr%jac_vertqi8(k)*z_AD(i8(k)) - rt(i8(k))=rt(i8(k))+lightptr%jac_vertti8(k)*z_AD(i8(k)) - z_AD(i8(k-1))=z_AD(i8(k-1))+z_AD(i8(k)) - z_AD(i8(k))=zero - - rq(i9(k))=rq(i9(k))+lightptr%jac_vertqi9(k)*z_AD(i9(k)) - rt(i9(k))=rt(i9(k))+lightptr%jac_vertti9(k)*z_AD(i9(k)) - z_AD(i9(k-1))=z_AD(i9(k-1))+z_AD(i9(k)) - z_AD(i9(k))=zero - - rq(i10(k))=rq(i10(k))+lightptr%jac_vertqi10(k)*z_AD(i10(k)) - rt(i10(k))=rt(i10(k))+lightptr%jac_vertti10(k)*z_AD(i10(k)) - z_AD(i10(k-1))=z_AD(i10(k-1))+z_AD(i10(k)) - z_AD(i10(k))=zero - - rq(i11(k))=rq(i11(k))+lightptr%jac_vertqi11(k)*z_AD(i11(k)) - rt(i11(k))=rt(i11(k))+lightptr%jac_vertti11(k)*z_AD(i11(k)) - z_AD(i11(k-1))=z_AD(i11(k-1))+z_AD(i11(k)) - z_AD(i11(k))=zero - - rq(i12(k))=rq(i12(k))+lightptr%jac_vertqi12(k)*z_AD(i12(k)) - rt(i12(k))=rt(i12(k))+lightptr%jac_vertti12(k)*z_AD(i12(k)) - z_AD(i12(k-1))=z_AD(i12(k-1))+z_AD(i12(k)) - z_AD(i12(k))=zero - - enddo + do k=nsig-1,2,-1 + + rq(i1(k))=rq(i1(k))+lightptr%jac_vertqi1(k)*z_ad(i1(k)) + rt(i1(k))=rt(i1(k))+lightptr%jac_vertti1(k)*z_ad(i1(k)) + z_ad(i1(k-1))=z_ad(i1(k-1))+z_ad(i1(k)) + z_ad(i1(k))=zero + + rq(i2(k))=rq(i2(k))+lightptr%jac_vertqi2(k)*z_ad(i2(k)) + rt(i2(k))=rt(i2(k))+lightptr%jac_vertti12(k)*z_ad(i2(k)) + z_ad(i2(k-1))=z_ad(i2(k-1))+z_ad(i2(k)) + z_ad(i2(k))=zero + + rq(i3(k))=rq(i3(k))+lightptr%jac_vertqi3(k)*z_ad(i3(k)) + rt(i3(k))=rt(i3(k))+lightptr%jac_vertti3(k)*z_ad(i3(k)) + z_ad(i3(k-1))=z_ad(i3(k-1))+z_ad(i3(k)) + z_ad(i3(k))=zero + + rq(i4(k))=rq(i4(k))+lightptr%jac_vertqi4(k)*z_ad(i4(k)) + rt(i4(k))=rt(i4(k))+lightptr%jac_vertti4(k)*z_ad(i4(k)) + z_ad(i4(k-1))=z_ad(i4(k-1))+z_ad(i4(k)) + z_ad(i4(k))=zero + + rq(i5(k))=rq(i5(k))+lightptr%jac_vertqi5(k)*z_ad(i5(k)) + rt(i5(k))=rt(i5(k))+lightptr%jac_vertti5(k)*z_ad(i5(k)) + z_ad(i5(k-1))=z_ad(i5(k-1))+z_ad(i5(k)) + z_ad(i5(k))=zero + + rq(i6(k))=rq(i6(k))+lightptr%jac_vertqi6(k)*z_ad(i6(k)) + rt(i6(k))=rt(i6(k))+lightptr%jac_vertti6(k)*z_ad(i6(k)) + z_ad(i6(k-1))=z_ad(i6(k-1))+z_ad(i6(k)) + z_ad(i6(k))=zero + + rq(i7(k))=rq(i7(k))+lightptr%jac_vertqi7(k)*z_ad(i7(k)) + rt(i7(k))=rt(i7(k))+lightptr%jac_vertti7(k)*z_ad(i7(k)) + z_ad(i7(k-1))=z_ad(i7(k-1))+z_ad(i7(k)) + z_ad(i7(k))=zero + + rq(i8(k))=rq(i8(k))+lightptr%jac_vertqi8(k)*z_ad(i8(k)) + rt(i8(k))=rt(i8(k))+lightptr%jac_vertti8(k)*z_ad(i8(k)) + z_ad(i8(k-1))=z_ad(i8(k-1))+z_ad(i8(k)) + z_ad(i8(k))=zero + + rq(i9(k))=rq(i9(k))+lightptr%jac_vertqi9(k)*z_ad(i9(k)) + rt(i9(k))=rt(i9(k))+lightptr%jac_vertti9(k)*z_ad(i9(k)) + z_ad(i9(k-1))=z_ad(i9(k-1))+z_ad(i9(k)) + z_ad(i9(k))=zero + + rq(i10(k))=rq(i10(k))+lightptr%jac_vertqi10(k)*z_ad(i10(k)) + rt(i10(k))=rt(i10(k))+lightptr%jac_vertti10(k)*z_ad(i10(k)) + z_ad(i10(k-1))=z_ad(i10(k-1))+z_ad(i10(k)) + z_ad(i10(k))=zero + + rq(i11(k))=rq(i11(k))+lightptr%jac_vertqi11(k)*z_ad(i11(k)) + rt(i11(k))=rt(i11(k))+lightptr%jac_vertti11(k)*z_ad(i11(k)) + z_ad(i11(k-1))=z_ad(i11(k-1))+z_ad(i11(k)) + z_ad(i11(k))=zero + + rq(i12(k))=rq(i12(k))+lightptr%jac_vertqi12(k)*z_ad(i12(k)) + rt(i12(k))=rt(i12(k))+lightptr%jac_vertti12(k)*z_ad(i12(k)) + z_ad(i12(k-1))=z_ad(i12(k-1))+z_ad(i12(k)) + z_ad(i12(k))=zero + + enddo endif !Adjoint - lightptr => lightNode_nextcast(lightptr) + lightptr => lightnode_nextcast(lightptr) enddo ! do while (associated(lightptr)) diff --git a/src/gsi/intlwcp.f90 b/src/gsi/intlwcp.f90 index fac9c97941..4d7d5bda63 100644 --- a/src/gsi/intlwcp.f90 +++ b/src/gsi/intlwcp.f90 @@ -20,18 +20,18 @@ module intlwcpmod ! !$$$ end documentation block -use m_obsNode, only: obsNode -use m_lwcpNode, only: lwcpNode -use m_lwcpNode, only: lwcpNode_typecast -use m_lwcpNode, only: lwcpNode_nextcast -use m_obsdiagNode, only: obsdiagNode_set +use m_obsnode, only: obsnode +use m_lwcpnode, only: lwcpnode +use m_lwcpnode, only: lwcpnode_typecast +use m_lwcpnode, only: lwcpnode_nextcast +use m_obsdiagnode, only: obsdiagnode_set implicit none -PRIVATE -PUBLIC intlwcp +private +public intlwcp interface intlwcp; module procedure & - intlwcp_ + intlwcp_ end interface contains @@ -47,10 +47,10 @@ subroutine intlwcp_(lwcphead,rval,sval) ! ! program history log: ! 2017-06-28 Ting-Chi Wu - mimic the structure in intpw.f90 and intgps.f90 -! - intlwcp.f90 includes 2 TL/ADJ options +! - intlwcp.f90 includes 2 tl/ad options ! 1) when l_wcp_cwm = .false.: -! operator = f(T,P,q) -! 2) when l_wcp_cwm = .true. and CWM partition6: +! operator = f(t,p,q) +! 2) when l_wcp_cwm = .true. and cwm partition6: ! operator = f(ql,qr) partition6 ! ! input argument list: @@ -88,7 +88,7 @@ subroutine intlwcp_(lwcphead,rval,sval) implicit none ! Declare passed variables - class(obsNode), pointer, intent(in ) :: lwcphead + class(obsnode), pointer, intent(in ) :: lwcphead type(gsi_bundle) ,intent(in ) :: sval type(gsi_bundle) ,intent(inout) :: rval @@ -96,17 +96,17 @@ subroutine intlwcp_(lwcphead,rval,sval) integer(i_kind) k,ier,istatus integer(i_kind),dimension(nsig):: i1,i2,i3,i4 ! real(r_kind) penalty - real(r_kind) :: t_TL,p_TL,q_TL - real(r_kind) :: t_AD,p_AD,q_AD - real(r_kind) :: ql_TL,qr_TL - real(r_kind) :: ql_AD,qr_AD + real(r_kind) :: t_tl,p_tl,q_tl + real(r_kind) :: t_ad,p_ad,q_ad + real(r_kind) :: ql_tl,qr_tl + real(r_kind) :: ql_ad,qr_ad real(r_kind) val,w1,w2,w3,w4 real(r_kind) cg_lwcp,grad,p0,wnotgross,wgross,pg_lwcp real(r_kind),pointer,dimension(:) :: st, sp, sq real(r_kind),pointer,dimension(:) :: sql, sqr real(r_kind),pointer,dimension(:) :: rt, rp, rq real(r_kind),pointer,dimension(:) :: rql, rqr - type(lwcpNode), pointer :: lwcpptr + type(lwcpnode), pointer :: lwcpptr ! If no lwcp data return if(.not. associated(lwcphead))return @@ -116,69 +116,69 @@ subroutine intlwcp_(lwcphead,rval,sval) if (.not.l_wcp_cwm) then - call gsi_bundlegetpointer(sval,'tsen',st,istatus);ier=istatus+ier - call gsi_bundlegetpointer(sval,'prse',sp,istatus);ier=istatus+ier - call gsi_bundlegetpointer(sval,'q' ,sq,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'tsen',rt,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'prse',rp,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'q' ,rq,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'tsen',st,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'prse',sp,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'q' ,sq,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'tsen',rt,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'prse',rp,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'q' ,rq,istatus);ier=istatus+ier else - call gsi_bundlegetpointer(sval,'ql',sql,istatus);ier=istatus+ier - call gsi_bundlegetpointer(sval,'qr',sqr,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'ql',rql,istatus);ier=istatus+ier - call gsi_bundlegetpointer(rval,'qr',rqr,istatus);ier=istatus+ier - + call gsi_bundlegetpointer(sval,'ql',sql,istatus);ier=istatus+ier + call gsi_bundlegetpointer(sval,'qr',sqr,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'ql',rql,istatus);ier=istatus+ier + call gsi_bundlegetpointer(rval,'qr',rqr,istatus);ier=istatus+ier + endif ! l_wcp_cwm if(ier/=0)return !lwcpptr => lwcphead - lwcpptr => lwcpNode_typecast(lwcphead) + lwcpptr => lwcpnode_typecast(lwcphead) do while (associated(lwcpptr)) w1=lwcpptr%wij(1) w2=lwcpptr%wij(2) w3=lwcpptr%wij(3) w4=lwcpptr%wij(4) - do k=1,nsig - i1(k)=lwcpptr%ij(1,k) - i2(k)=lwcpptr%ij(2,k) - i3(k)=lwcpptr%ij(3,k) - i4(k)=lwcpptr%ij(4,k) - enddo + do k=1,nsig + i1(k)=lwcpptr%ij(1,k) + i2(k)=lwcpptr%ij(2,k) + i3(k)=lwcpptr%ij(3,k) + i4(k)=lwcpptr%ij(4,k) + enddo val=zero ! Forward model if (.not.l_wcp_cwm) then - do k=1,nsig - t_TL=w1* st(i1(k))+w2* st(i2(k))+w3* st(i3(k))+w4* st(i4(k)) - p_TL=w1* sp(i1(k))+w2* sp(i2(k))+w3* sp(i3(k))+w4* sp(i4(k)) - q_TL=w1* sq(i1(k))+w2* sq(i2(k))+w3* sq(i3(k))+w4* sq(i4(k)) - val = val + ( t_TL*lwcpptr%jac_t(k) + & - p_TL*lwcpptr%jac_p(k) + & - q_TL*lwcpptr%jac_q(k) ) ! tpwcon*r10*(piges(k)-piges(k+1)) already did in setuplwcp.f90 - end do + do k=1,nsig + t_tl=w1* st(i1(k))+w2* st(i2(k))+w3* st(i3(k))+w4* st(i4(k)) + p_tl=w1* sp(i1(k))+w2* sp(i2(k))+w3* sp(i3(k))+w4* sp(i4(k)) + q_tl=w1* sq(i1(k))+w2* sq(i2(k))+w3* sq(i3(k))+w4* sq(i4(k)) + val = val + ( t_tl*lwcpptr%jac_t(k) + & + p_tl*lwcpptr%jac_p(k) + & + q_tl*lwcpptr%jac_q(k) ) ! tpwcon*r10*(piges(k)-piges(k+1)) already did in setuplwcp.f90 + end do else - do k=1,nsig - ql_TL=w1* sql(i1(k))+w2* sql(i2(k))+w3* sql(i3(k))+w4* sql(i4(k)) - qr_TL=w1* sqr(i1(k))+w2* sqr(i2(k))+w3* sqr(i3(k))+w4* sqr(i4(k)) - val = val + ( ql_TL*lwcpptr%jac_ql(k) + & - qr_TL*lwcpptr%jac_qr(k) ) ! tpwcon*r10*(piges(k)-piges(k+1)) already did in setuplwcp.f90 - end do + do k=1,nsig + ql_tl=w1* sql(i1(k))+w2* sql(i2(k))+w3* sql(i3(k))+w4* sql(i4(k)) + qr_tl=w1* sqr(i1(k))+w2* sqr(i2(k))+w3* sqr(i3(k))+w4* sqr(i4(k)) + val = val + ( ql_tl*lwcpptr%jac_ql(k) + & + qr_tl*lwcpptr%jac_qr(k) ) ! tpwcon*r10*(piges(k)-piges(k+1)) already did in setuplwcp.f90 + end do endif ! l_wcp_cwm if(luse_obsdiag)then if (lsaveobsens) then grad = val*lwcpptr%raterr2*lwcpptr%err2 !-- lwcpptr%diags%obssen(jiter) = grad - call obsdiagNode_set(lwcpptr%diags,jiter=jiter,obssen=grad) + call obsdiagnode_set(lwcpptr%diags,jiter=jiter,obssen=grad) else !-- if (lwcpptr%luse) lwcpptr%diags%tldepart(jiter)=val - if (lwcpptr%luse) call obsdiagNode_set(lwcpptr%diags,jiter=jiter,tldepart=val) + if (lwcpptr%luse) call obsdiagnode_set(lwcpptr%diags,jiter=jiter,tldepart=val) endif end if @@ -206,42 +206,42 @@ subroutine intlwcp_(lwcphead,rval,sval) ! Adjoint if (.not.l_wcp_cwm) then - do k=1,nsig - t_AD = grad*lwcpptr%jac_t(k) - rt(i1(k))=rt(i1(k))+w1*t_AD - rt(i2(k))=rt(i2(k))+w2*t_AD - rt(i3(k))=rt(i3(k))+w3*t_AD - rt(i4(k))=rt(i4(k))+w4*t_AD - p_AD = grad*lwcpptr%jac_p(k) - rp(i1(k))=rp(i1(k))+w1*p_AD - rp(i2(k))=rp(i2(k))+w2*p_AD - rp(i3(k))=rp(i3(k))+w3*p_AD - rp(i4(k))=rp(i4(k))+w4*p_AD - q_AD = grad*lwcpptr%jac_q(k) - rq(i1(k))=rq(i1(k))+w1*q_AD - rq(i2(k))=rq(i2(k))+w2*q_AD - rq(i3(k))=rq(i3(k))+w3*q_AD - rq(i4(k))=rq(i4(k))+w4*q_AD - enddo + do k=1,nsig + t_ad = grad*lwcpptr%jac_t(k) + rt(i1(k))=rt(i1(k))+w1*t_ad + rt(i2(k))=rt(i2(k))+w2*t_ad + rt(i3(k))=rt(i3(k))+w3*t_ad + rt(i4(k))=rt(i4(k))+w4*t_ad + p_ad = grad*lwcpptr%jac_p(k) + rp(i1(k))=rp(i1(k))+w1*p_ad + rp(i2(k))=rp(i2(k))+w2*p_ad + rp(i3(k))=rp(i3(k))+w3*p_ad + rp(i4(k))=rp(i4(k))+w4*p_ad + q_ad = grad*lwcpptr%jac_q(k) + rq(i1(k))=rq(i1(k))+w1*q_ad + rq(i2(k))=rq(i2(k))+w2*q_ad + rq(i3(k))=rq(i3(k))+w3*q_ad + rq(i4(k))=rq(i4(k))+w4*q_ad + enddo else - do k=1,nsig - ql_AD = grad*lwcpptr%jac_ql(k) - rql(i1(k))=rql(i1(k))+w1*ql_AD - rql(i2(k))=rql(i2(k))+w2*ql_AD - rql(i3(k))=rql(i3(k))+w3*ql_AD - rql(i4(k))=rql(i4(k))+w4*ql_AD - qr_AD = grad*lwcpptr%jac_qr(k) - rqr(i1(k))=rqr(i1(k))+w1*qr_AD - rqr(i2(k))=rqr(i2(k))+w2*qr_AD - rqr(i3(k))=rqr(i3(k))+w3*qr_AD - rqr(i4(k))=rqr(i4(k))+w4*qr_AD - enddo + do k=1,nsig + ql_ad = grad*lwcpptr%jac_ql(k) + rql(i1(k))=rql(i1(k))+w1*ql_ad + rql(i2(k))=rql(i2(k))+w2*ql_ad + rql(i3(k))=rql(i3(k))+w3*ql_ad + rql(i4(k))=rql(i4(k))+w4*ql_ad + qr_ad = grad*lwcpptr%jac_qr(k) + rqr(i1(k))=rqr(i1(k))+w1*qr_ad + rqr(i2(k))=rqr(i2(k))+w2*qr_ad + rqr(i3(k))=rqr(i3(k))+w3*qr_ad + rqr(i4(k))=rqr(i4(k))+w4*qr_ad + enddo endif ! l_wcp_cwm endif ! l_do_adjoint !lwcpptr => lwcpptr%llpoint - lwcpptr => lwcpNode_nextcast(lwcpptr) + lwcpptr => lwcpnode_nextcast(lwcpptr) end do