module stplightmod !$$$ module documentation block ! . . . . ! module: stplightmod module for stplight and its tangent linear stplight_tl ! prgmmr: k apodaca ! org: CSU/CIRA, Data Assimilation Group ! date: 2016-05-19 ! ! abstract: module for calculating stplight and its tangent linear stplight_tl ! ! program history log: ! 2016-05-19 apodaca - original version ! ! subroutines included: ! sub stplight ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block implicit none PRIVATE PUBLIC stplight contains subroutine stplight(lighthead,rval,sval,out,sges,nstep) !$$$ subprogram documentation block ! . . . . ! subprogram: stplight compute contribution to penalty and stepsize ! from lightning, using nonlinear qc ! prgmmr: apodaca org: CSU/CIRA date: 2016-05-19 ! ! abstract: This routine applies the (linear) operator of the ! lightning flash rate model and finds an optimal estimate ! of step size as done in steepest descent or conjugate ! gradient algorithms. Note that in the case for lightning ! flash rate dk=-delJ(h(xk)) ! This version includes nonlinear qc. ! ! program history log: ! 2016-05-19 k apodaca ! 2016-06-21 k apodaca - update documentation ! 2018-03-07 k apodaca - replaced ob_type with polymorphic obsNode through type casting ! 2018-08-27 k apodaca - add TL-related components of a second oservation operator for ! lightning observations suitable for non-hydrostatic, ! cloud-resolving models with additional ice-phase hydrometeor ! control variables ! ! input argument list: ! lighthead ! rt - search direction (gradxJ) for virtual temperature ! rq - search direction (gradxJ) for specific humidity ! rqi - search direction (gradxJ) for cloud ice ! rqs - search direction (gradxJ) for snow ! rqg - search direction (gradxJ) for graupel ! ru - search direction (gradxJ) for the u-component of wind ! rv - search direction (gradxJ) for the v-component of wind ! st - analysis increment (correction) for virtual temperature ! sq - analysis increment (correction) for specific humidity ! sqi - analysis increment (correction) for cloud ice ! sqs - analysis increment (correction) for snow ! sqg - analysis increment (correction) for graupel ! su - analysis increment (correction) for u-component of wind ! sv - analysis increment (correction) for the v-component of wind ! sges - stepsize estimates (nstep) ! nstep - number of stepsize estimates (==0 means use outer iteration values) ! ! output argument list: ! out(1:nstep)- contribution to penalty from lightning - sges(1:nstep) ! ! attributes: ! language: Fortran 90 and/or above ! machine: ! !$$$ use kinds, only: r_kind,i_kind,r_quad use qcmod, only: nlnqc_iter,varqc_iter use constants, only: zero,one,two,half,tiny_r_kind,cg_term,zero_quad,r3600,fv use gridmod, only: nsig use gridmod, only: wrf_mass_regional,regional use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use m_obsNode , only: obsNode use m_lightNode , only: lightNode use m_lightNode , only: lightNode_typecast use m_lightNode , only: lightNode_nextcast implicit none ! Declare passed variables class(obsNode ),pointer ,intent(in ) :: lighthead integer(i_kind) ,intent(in ) :: nstep real(r_quad),dimension(max(1,nstep)),intent(inout) :: out type(gsi_bundle) ,intent(in ) :: rval,sval real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges ! Declare local variables integer(i_kind) k,kk,ier,istatus integer(i_kind),dimension(nsig) :: i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12 real(r_kind) :: val,val2 real(r_kind) :: w1,w2,w3,w4 real(r_kind),pointer,dimension(:) :: st,sq,sqi,sqs,sqg,su,sv real(r_kind),pointer,dimension(:) :: rt,rq,rqi,rqs,rqg,ru,rv type(lightNode), pointer :: lightptr ! Variables for TL and r_TL 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) :: htot_TL real(r_kind) :: htoti1_TL,htoti2_TL,htoti3_TL,htoti4_TL real(r_kind) :: totice_colinti1_TL,totice_colinti2_TL,totice_colinti3_TL real(r_kind) :: totice_colinti4_TL real(r_kind) :: wmax real(r_kind),parameter :: k3=0.95_r_kind real(r_kind),dimension(1:nsig) :: rz_TL real(r_kind),dimension(1:nsig) :: rhoriz_adv_TL real(r_kind),dimension(1:nsig) :: rvert_adv_TL real(r_kind),dimension(1:nsig) :: rw_TL real(r_kind) :: rwmaxi1_tl,rwmaxi2_tl,rwmaxi3_tl,rwmaxi4_tl real(r_kind) :: rflashrate_tl,rflashratei1_tl,rflashratei2_tl real(r_kind) :: rflashratei3_tl,rflashratei4_tl real(r_kind) :: rh1i1_TL,rh1i2_TL,rh1i3_TL,rh1i4_TL real(r_kind) :: rh2i1_TL,rh2i2_TL,rh2i3_TL,rh2i4_TL real(r_kind) :: rhtot_TL real(r_kind) :: rhtoti1_TL,rhtoti2_TL,rhtoti3_TL,rhtoti4_TL real(r_kind) :: rtotice_colinti1_TL,rtotice_colinti2_TL,rtotice_colinti3_TL real(r_kind) :: rtotice_colinti4_TL real(r_kind) cg_light,wgross,wnotgross real(r_kind) pg_light,nref real(r_kind),dimension(max(1,nstep))::pen ! Initialize penalty, b1, and b3 to zero out=zero_quad ! If no light data return if(.not. associated(lighthead))return ! Retrieve pointers ! Simply return if any pointer not found ier=0 call gsi_bundlegetpointer(sval,'tv', st,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'q', sq,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'u', su,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'v', sv,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'tv', rt,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'q', rq,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'u', ru,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'v', rv,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'qi',sqi,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'qi',rqi,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'qg',sqg,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'qg',rqg,istatus);ier=istatus+ier call gsi_bundlegetpointer(sval,'qs',sqs,istatus);ier=istatus+ier call gsi_bundlegetpointer(rval,'qs',rqs,istatus);ier=istatus+ier if(ier/=0)return ! Loop over observations lightptr => lightNode_typecast(lighthead) do while (associated(lightptr)) if(lightptr%luse)then val2=-lightptr%res if(nstep > 0)then do k=1,nsig i1(k)=lightptr%ij(1,k) i2(k)=lightptr%ij(2,k) i3(k)=lightptr%ij(3,k) i4(k)=lightptr%ij(4,k) i5(k)=lightptr%ij(5,k) i6(k)=lightptr%ij(6,k) i7(k)=lightptr%ij(7,k) i8(k)=lightptr%ij(8,k) i9(k)=lightptr%ij(9,k) i10(k)=lightptr%ij(10,k) i11(k)=lightptr%ij(11,k) i12(k)=lightptr%ij(12,k) end do w1=lightptr%wij(1) w2=lightptr%wij(2) w3=lightptr%wij(3) w4=lightptr%wij(4) val=zero ! . . . . ! 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 ! the observation operator for lightning flash rate. Calculations are done ! at each quadrant, i.e., central, north, south, east, and west. ! ! i6-------i8 ! | | ! | | ! i10-----i2-------i4------i12 ! | | | | ! | | + | | ! i9------i1-------i3------i11 ! | | ! | | ! i5-------i7 ! ! . . . . ! In the following section, the tangent linear of the lightning flash rate observation ! operator is calculated by being broken into parts. ! Tangent linear of height (z) 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)) & +lightptr%jac_vertqi1(k)*sq(i1(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)) & +lightptr%jac_vertqi3(k)*sq(i3(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)) & +lightptr%jac_vertqi5(k)*sq(i5(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)) & +lightptr%jac_vertqi7(k)*sq(i7(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)) & +lightptr%jac_vertqi9(k)*sq(i9(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)) & +lightptr%jac_vertqi11(k)*sq(i11(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 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_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))) 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))) 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))) enddo ! do k=2,nsig-1 ! Tangent Linear of the Vertical Advection Section ! Variable Initialization 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)* & (((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)* & (((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)* & (((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)* & (((one+fv*lightptr%jac_qi4(k))*st(i4(k))) & +(lightptr%jac_ti4(k)*fv*sq(i4(k)))) ! 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)) enddo !do k=1,nsig-1 ! . . . . ! Tangent Linear of lightning flash rate ! . . . . ! Regional if (regional) then ! WRF-ARW if (wrf_mass_regional) then ! Tangent linear - Lightning flash rate as a function of ! vertical graupel flux within the mixed-phase region ! (-15 lightptr%jac_qgmbi1(lightptr%kboti) 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)) else 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)) else 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)) else 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)) else h1i4_TL=zero endif ! 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 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_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_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 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 ! 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) val2 = val2 + htot_TL endif ! wrf_mass_regional endif !if (regional) then ! . . . . ! Global if (.not. regional) then ! 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. 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 endif if (lightptr%jac_wmaxflagi2) then wmax=-1.e+10_r_kind do k=1,nsig-1 if (w_TL(i2(k)) > wmax) then lightptr%jac_kverti2=k wmaxi2_TL=w_TL(i2(lightptr%jac_kverti2)) endif 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 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 endif ! 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 ! 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) val2 = val2 + flashrate_TL endif ! end regional/global block ! . . . . ! Search direction (gradxJ) for lightning flash rate observation ! operator ! gradxJ: Tangent linear of height (z) ! Variable Initialization rz_TL(:)=zero rhoriz_adv_TL(:)=zero rvert_adv_TL(:)=zero rw_TL(:)=zero do k=2,nsig-1 rz_TL(i1(1))=lightptr%jac_z0i1 rz_TL(i2(1))=lightptr%jac_z0i2 rz_TL(i3(1))=lightptr%jac_z0i3 rz_TL(i4(1))=lightptr%jac_z0i4 rz_TL(i5(1))=lightptr%jac_z0i5 rz_TL(i6(1))=lightptr%jac_z0i6 rz_TL(i7(1))=lightptr%jac_z0i7 rz_TL(i8(1))=lightptr%jac_z0i8 rz_TL(i9(1))=lightptr%jac_z0i9 rz_TL(i10(1))=lightptr%jac_z0i10 rz_TL(i11(1))=lightptr%jac_z0i11 rz_TL(i12(1))=lightptr%jac_z0i12 rz_TL(i1(k))=rz_TL(i1(k-1))+lightptr%jac_vertti1(k)*rt(i1(k)) & +lightptr%jac_vertqi1(k)*rq(i1(k)) rz_TL(i2(k))=rz_TL(i2(k-1))+lightptr%jac_vertti2(k)*rt(i2(k)) & +lightptr%jac_vertqi2(k)*rq(i2(k)) rz_TL(i3(k))=rz_TL(i3(k-1))+lightptr%jac_vertti3(k)*rt(i3(k)) & +lightptr%jac_vertqi3(k)*rq(i3(k)) rz_TL(i4(k))=rz_TL(i4(k-1))+lightptr%jac_vertti4(k)*rt(i4(k)) & +lightptr%jac_vertqi4(k)*rq(i4(k)) rz_TL(i5(k))=rz_TL(i5(k-1))+lightptr%jac_vertti5(k)*rt(i5(k)) & +lightptr%jac_vertqi5(k)*rq(i5(k)) rz_TL(i6(k))=rz_TL(i6(k-1))+lightptr%jac_vertti6(k)*rt(i6(k)) & +lightptr%jac_vertqi6(k)*rq(i6(k)) rz_TL(i7(k))=rz_TL(i7(k-1))+lightptr%jac_vertti7(k)*rt(i7(k)) & +lightptr%jac_vertqi7(k)*rq(i7(k)) rz_TL(i8(k))=rz_TL(i8(k-1))+lightptr%jac_vertti8(k)*rt(i8(k)) & +lightptr%jac_vertqi8(k)*rq(i8(k)) rz_TL(i9(k))=rz_TL(i9(k-1))+lightptr%jac_vertti9(k)*rt(i9(k)) & +lightptr%jac_vertqi9(k)*rq(i9(k)) rz_TL(i10(k))=rz_TL(i10(k-1))+lightptr%jac_vertti10(k)*rt(i10(k)) & +lightptr%jac_vertqi10(k)*rq(i10(k)) rz_TL(i11(k))=rz_TL(i11(k-1))+lightptr%jac_vertti11(k)*rt(i11(k)) & +lightptr%jac_vertqi11(k)*rq(i11(k)) rz_TL(i12(k))=rz_TL(i12(k-1))+lightptr%jac_vertti12(k)*rt(i12(k)) & +lightptr%jac_vertqi12(k)*rq(i12(k)) ! gradxJ: Tangent Linear of the Horizontal Advection Section rhoriz_adv_TL(i1(k))=lightptr%jac_zdxi1(k)*ru(i1(k)) & +lightptr%jac_zdyi1(k)*rv(i1(k)) & +lightptr%jac_udxi1(k)*(rz_TL(i3(k))-rz_TL(i9(k))) & +lightptr%jac_vdyi1(k)*(rz_TL(i2(k))-rz_TL(i5(k))) rhoriz_adv_TL(i2(k))=lightptr%jac_zdxi2(k)*ru(i2(k)) & +lightptr%jac_zdyi2(k)*rv(i2(k)) & +lightptr%jac_udxi2(k)*(rz_TL(i4(k))-rz_TL(i10(k)))& +lightptr%jac_vdyi2(k)*(rz_TL(i6(k))-rz_TL(i1 (k))) rhoriz_adv_TL(i3(k))=lightptr%jac_zdxi3(k)*ru(i3(k)) & +lightptr%jac_zdyi3(k)*rv(i3(k)) & +lightptr%jac_udxi3(k)*(rz_TL(i11(k))-rz_TL(i1(k))) & +lightptr%jac_vdyi3(k)*(rz_TL(i4 (k))-rz_TL(i7(k))) rhoriz_adv_TL(i4(k))=lightptr%jac_zdxi4(k)*ru(i4(k)) & +lightptr%jac_zdyi4(k)*rv(i4(k)) & +lightptr%jac_udxi4(k)*(rz_TL(i12(k))-rz_TL(i2(k))) & +lightptr%jac_vdyi4(k)*(rz_TL(i8 (k))-rz_TL(i3(k))) enddo !do k=2,nsig-1 ! gradxJ: Tangent Linear of the Vertical Advection Section do k=1,nsig-1 rvert_adv_TL(i1(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti1(k)* & (((one+fv*lightptr%jac_qi1(k))*rt(i1(k))) & +(lightptr%jac_ti1(k)*fv*rq(i1(k)))) rvert_adv_TL(i2(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti2(k)* & (((one+fv*lightptr%jac_qi2(k))*rt(i2(k))) & +(lightptr%jac_ti2(k)*fv*rq(i2(k)))) rvert_adv_TL(i3(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti3(k)* & (((one+fv*lightptr%jac_qi3(k))*rt(i3(k))) & +(lightptr%jac_ti3(k)*fv*rq(i3(k)))) rvert_adv_TL(i4(k))=-lightptr%jac_vert(k)*lightptr%jac_sigdoti4(k)* & (((one+fv*lightptr%jac_qi4(k))*rt(i4(k))) & +(lightptr%jac_ti4(k)*fv*rq(i4(k)))) ! gradxJ: Tangent Linear of Vertical Velocity rw_TL(i1(k))=rhoriz_adv_TL(i1(k))+rvert_adv_TL(i1(k)) rw_TL(i2(k))=rhoriz_adv_TL(i2(k))+rvert_adv_TL(i2(k)) rw_TL(i3(k))=rhoriz_adv_TL(i3(k))+rvert_adv_TL(i3(k)) rw_TL(i4(k))=rhoriz_adv_TL(i4(k))+rvert_adv_TL(i4(k)) enddo ! . . . . ! gradxJ: Tangent Linear of lightning flash rate ! . . . . ! Regional if (regional) then ! WRF-ARW if (wrf_mass_regional) then ! gradxJ: Tangent linear - Lightning flash rate as a function of ! vertical graupel flux within the mixed-phase region ! (-15 lightptr%jac_qgmbi1(lightptr%kboti1(k)deg C) if (lightptr%kboti1 > zero) then rh1i1_TL=lightptr%jac_qgmai1(lightptr%kboti1)*rqg(i1(lightptr%kboti1))+& lightptr%jac_qgmbi1(lightptr%kboti1)*& (half*(rw_TL(i1(lightptr%kboti1))+rw_TL(i1(lightptr%kboti1+1)))) rh1i1_TL=rh1i1_TL/(abs(rh1i1_TL)) else rh1i1_TL=zero endif if (lightptr%kboti2 > 0) then rh1i2_TL=lightptr%jac_qgmai2(lightptr%kboti2)*rqg(i2(lightptr%kboti2))+& lightptr%jac_qgmbi2(lightptr%kboti2)*& (half*(rw_TL(i2(lightptr%kboti2))+rw_TL(i2(lightptr%kboti2+1)))) rh1i2_TL=rh1i2_TL/(abs(rh1i2_TL)) else rh1i2_TL=zero endif if (lightptr%kboti3 > zero) then rh1i3_TL=lightptr%jac_qgmai3(lightptr%kboti3)*rqg(i3(lightptr%kboti3))+& lightptr%jac_qgmbi3(lightptr%kboti3)*& (half*(rw_TL(i3(lightptr%kboti3))+rw_TL(i3(lightptr%kboti3+1)))) rh1i3_TL=rh1i3_TL/(abs(rh1i3_TL)) else rh1i3_TL=zero endif if (lightptr%kboti4 > zero) then rh1i4_TL=lightptr%jac_qgmai4(lightptr%kboti4)*rqg(i4(lightptr%kboti4))+& lightptr%jac_qgmbi4(lightptr%kboti4)*& (half*(rw_TL(i4(lightptr%kboti4))+rw_TL(i4(lightptr%kboti4+1)))) rh1i4_TL=rh1i4_TL/(abs(rh1i4_TL)) else rh1i4_TL=zero endif ! gradxJ: Tangent Linear - Lightning flash rate as a function of total column-integrated ! ice-phase hydrometeors do k=1,nsig-1 rtotice_colinti1_TL = rtotice_colinti1_TL+lightptr%jac_icei1(k) * & (rqi(i1(k))+rqs(i1(k))+rqg(i1(k)))+& lightptr%jac_zicei1(k)*rz_TL(i1(k)) rtotice_colinti2_TL = rtotice_colinti2_TL+lightptr%jac_icei2(k) * & (rqi(i2(k))+rqs(i2(k))+rqg(i2(k)))+& lightptr%jac_zicei2(k)*rz_TL(i2(k)) rtotice_colinti3_TL = rtotice_colinti3_TL+lightptr%jac_icei3(k) * & (rqi(i3(k))+rqs(i3(k))+rqg(i3(k)))+& lightptr%jac_zicei3(k)*rz_TL(i3(k)) rtotice_colinti4_TL = rtotice_colinti4_TL+lightptr%jac_icei4(k) * & (rqi(i4(k))+rqs(i4(k))+rqg(i4(k)))+& lightptr%jac_zicei4(k)*rz_TL(i4(k)) enddo rh2i1_TL=(1-k3)*rtotice_colinti1_TL rh2i2_TL=(1-k3)*rtotice_colinti2_TL rh2i3_TL=(1-k3)*rtotice_colinti3_TL rh2i4_TL=(1-k3)*rtotice_colinti4_TL rhtoti1_TL= rh1i1_TL+rh2i1_TL rhtoti2_TL= rh1i2_TL+rh2i2_TL rhtoti3_TL= rh1i3_TL+rh2i3_TL rhtoti4_TL= rh1i4_TL+rh2i4_TL ! Interpolation of lightning flash rate to observation location (2D field) ! Forward Model rhtot_TL = (w1*rhtoti1_TL + w2*rhtoti2_TL + & w3*rhtoti3_TL + w4*rhtoti4_TL) val = val + rhtot_TL endif ! wrf_mass_regional endif !if (regional) then ! . . . . ! Global if (.not. regional) then ! Cloud Mask (gradxJ) ! 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. rwmaxi1_TL=zero rwmaxi2_TL=zero rwmaxi3_TL=zero rwmaxi4_TL=zero if (lightptr%jac_wmaxflagi1) then wmax=-1.e+10_r_kind do k=1,nsig-1 if (rw_TL(i1(k)) > wmax) then lightptr%jac_kverti1=k rwmaxi1_TL=rw_TL(i1(lightptr%jac_kverti1)) endif if (rwmaxi1_TL < zero) then rwmaxi1_TL=zero endif enddo ! k loop endif if (lightptr%jac_wmaxflagi2) then wmax=-1.e+10_r_kind do k=1,nsig-1 if (rw_TL(i2(k)) > wmax) then lightptr%jac_kverti2=k rwmaxi2_TL=rw_TL(i2(lightptr%jac_kverti2)) endif if (rwmaxi2_TL < zero) then rwmaxi2_TL=zero endif enddo ! k loop endif if (lightptr%jac_wmaxflagi3) then wmax=-1.e+10_r_kind do k=1,nsig-1 if (rw_TL(i3(k)) > wmax) then lightptr%jac_kverti3=k rwmaxi3_TL=rw_TL(i3(lightptr%jac_kverti3)) endif if (rwmaxi3_TL < zero) then rwmaxi3_TL=zero endif enddo ! k loop endif if (lightptr%jac_wmaxflagi4) then wmax=-1.e+10_r_kind do k=1,nsig-1 if (rw_TL(i4(k)) > wmax) then lightptr%jac_kverti4=k rwmaxi4_TL=rw_TL(i4(lightptr%jac_kverti4)) endif if (rwmaxi4_TL < zero) then rwmaxi4_TL=zero endif enddo ! k loop endif ! gradxJ: Tangent Linear of Lightning Flash Rate rflashratei1_TL=lightptr%jac_fratei1*wmaxi1_TL rflashratei2_TL=lightptr%jac_fratei1*wmaxi2_TL rflashratei3_TL=lightptr%jac_fratei1*wmaxi3_TL rflashratei4_TL=lightptr%jac_fratei1*wmaxi4_TL ! Interpolation of lightning flash rate TL to observation location (2D field) ! Forward Model rflashrate_TL = (w1*rflashratei1_TL + w2*rflashratei2_TL & +w3*rflashratei3_TL + w4*rflashratei4_TL) val = val + rflashrate_TL endif ! end regional/global block ! . . . . ! penalty and gradient do kk=1,nstep nref=val2+sges(kk)*val pen(kk)=nref*nref*lightptr%err2 end do else pen(1)=val2*val2*lightptr%err2 end if !if(nstep > 0)then ! Modify penalty term if nonlinear QC if (nlnqc_iter .and. lightptr%pg > tiny_r_kind .and. lightptr%b > tiny_r_kind) then pg_light=lightptr%pg*varqc_iter cg_light=cg_term/lightptr%b wnotgross= one-pg_light wgross = pg_light*cg_light/wnotgross do kk=1,max(1,nstep) pen(kk) = -two*log((exp(-half*pen(kk)) + wgross)/(one+wgross)) end do endif ! Cost function out(1) = out(1)+pen(1)*lightptr%raterr2 do kk=2,nstep out(kk) = out(kk)+(pen(kk)-pen(1))*lightptr%raterr2 end do endif ! if(lightptr%luse) then lightptr => lightNode_nextcast(lightptr) end do !do while (associated(lightptr)) return end subroutine stplight end module stplightmod