MODULE tl_convolution_mod ! !git $Id$ !svn $Id: tl_convolution.F 1180 2023-07-13 02:42:10Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group Andrew M. Moore ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! This routine performs a spatial convolution of the tangeny state ! ! solution to model the background error correlations, C, using ! ! a generalized diffusion operator. This allows the observational ! ! information to spread spatially in 4DVAR data assimilation. ! ! ! ! The background error covariance is defined as: ! ! ! ! B = S C S ! ! ! ! C = C^(1/2) C^(T/2) ! ! ! ! C^(1/2) = G L^(1/2) W^(-1/2) TLM ! ! C^(T/2) = W^(-1/2) L^(T/2) G ADM ! ! ! ! where ! ! ! ! B : background-error covariance matrix ! ! S : diagonal matrix of background-error standard deviations ! ! C : symmetric matrix of background-error correlations ! ! G : normalization coefficients matrix used to ensure that the ! ! diagonal variances of C are equal to unity. ! ! L : tangent linear and adjoint diffusion operators ! ! W : diagonal matrix of local area or volume metrics used to ! ! convert L into a symmetric matrix: LW^(-1). ! ! ! ! Here, T/2 denote the transpose if a squared-root factor. ! ! ! ! This routine is used to provide a better preconditioning of the ! ! minimization problem, which is expressed as a function of a new ! ! state vector, v, given by: ! ! ! ! v = B^(-1/2) delta_x (v-space) ! ! or ! ! delta_x = B^(1/2) v ! ! ! ! where ! ! ! ! B = B^(T/2) B^(1/2) ! ! ! ! Therefore, the cost function, J, gradient becomes: ! ! ! ! GRAD_v(J) = v + {B^(T/2)} GRAD_x(J) ! ! ! ! In incremental 4DVAR, these spatial convolutions constitutes a ! ! smoothing action on the correlation operator and they are used ! ! to transform between model space to minimization space and vice ! ! versa: ! ! ! ! ad_convolution compute GRAD_v(J) from GRAD_x(J) ! ! tl_convolution compute x from v ! ! ! ! The minimization of of J in the descent algorithm is in v-space. ! ! ! ! Reference: ! ! ! ! Weaver, A. and P. Courtier, 2001: Correlation modeling on the ! ! sphere using a generalized diffusion equation, Q.J.R. Meteo. ! ! Soc, 127, 1815-1846. ! ! ! !======================================================================! ! USE mod_kinds ! implicit none ! PRIVATE PUBLIC :: tl_convolution ! CONTAINS ! !*********************************************************************** SUBROUTINE tl_convolution (ng, tile, Linp, Lweak, ifac) !*********************************************************************** ! USE mod_param USE mod_grid USE mod_mixing USE mod_ocean USE mod_stepping ! ! Imported variable declarations. ! logical, intent(in) :: Lweak integer, intent(in) :: ng, tile, Linp, ifac ! ! Local variable declarations. ! integer :: IminS, ImaxS, JminS, JmaxS integer :: LBi, UBi, LBj, UBj, LBij, UBij ! ! Set horizontal starting and ending indices for automatic private ! storage arrays. ! IminS=BOUNDS(ng)%Istr(tile)-3 ImaxS=BOUNDS(ng)%Iend(tile)+3 JminS=BOUNDS(ng)%Jstr(tile)-3 JmaxS=BOUNDS(ng)%Jend(tile)+3 ! ! Determine array lower and upper bounds in the I- and J-directions. ! LBi=BOUNDS(ng)%LBi(tile) UBi=BOUNDS(ng)%UBi(tile) LBj=BOUNDS(ng)%LBj(tile) UBj=BOUNDS(ng)%UBj(tile) ! ! Set array lower and upper bounds for MIN(I,J) directions and ! MAX(I,J) directions. ! LBij=BOUNDS(ng)%LBij UBij=BOUNDS(ng)%UBij ! CALL tl_convolution_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), Linp, Lweak, ifac, & & GRID(ng) % pm, & & GRID(ng) % om_p, & & GRID(ng) % om_r, & & GRID(ng) % om_u, & & GRID(ng) % om_v, & & GRID(ng) % pn, & & GRID(ng) % on_p, & & GRID(ng) % on_r, & & GRID(ng) % on_u, & & GRID(ng) % on_v, & & GRID(ng) % pmon_p, & & GRID(ng) % pmon_r, & & GRID(ng) % pmon_u, & & GRID(ng) % pnom_p, & & GRID(ng) % pnom_r, & & GRID(ng) % pnom_v, & & GRID(ng) % rmask, & & GRID(ng) % pmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & & GRID(ng) % h, & & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & & MIXING(ng) % Kh, & & MIXING(ng) % Kv, & & OCEAN(ng) % b_t, & & OCEAN(ng) % b_u, & & OCEAN(ng) % b_v, & & OCEAN(ng) % b_zeta, & & OCEAN(ng) % b_ubar, & & OCEAN(ng) % b_vbar, & & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & & OCEAN(ng) % tl_zeta) RETURN END SUBROUTINE tl_convolution ! !*********************************************************************** SUBROUTINE tl_convolution_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, Linp, Lweak, ifac, & & pm, om_p, om_r, om_u, om_v, & & pn, on_p, on_r, on_u, on_v, & & pmon_p, pmon_r, pmon_u, & & pnom_p, pnom_r, pnom_v, & & rmask, pmask, umask, vmask, & & h, & & Hz, z_r, z_w, & & Kh, & & Kv, & & VnormR, VnormU, VnormV, & & HnormR, HnormU, HnormV, & & tl_t, tl_u, tl_v, & & tl_ubar, tl_vbar, & & tl_zeta) !*********************************************************************** ! USE mod_param USE mod_fourdvar USE mod_ncparam USE mod_scalars ! USE tl_conv_2d_mod USE tl_conv_3d_mod USE mp_exchange_mod USE set_depth_mod ! ! Imported variable declarations. ! logical, intent(in) :: Lweak integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew, Linp, ifac ! real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: om_p(LBi:,LBj:) real(r8), intent(in) :: om_r(LBi:,LBj:) real(r8), intent(in) :: om_u(LBi:,LBj:) real(r8), intent(in) :: om_v(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) real(r8), intent(in) :: on_p(LBi:,LBj:) real(r8), intent(in) :: on_r(LBi:,LBj:) real(r8), intent(in) :: on_u(LBi:,LBj:) real(r8), intent(in) :: on_v(LBi:,LBj:) real(r8), intent(in) :: pmon_p(LBi:,LBj:) real(r8), intent(in) :: pmon_r(LBi:,LBj:) real(r8), intent(in) :: pmon_u(LBi:,LBj:) real(r8), intent(in) :: pnom_p(LBi:,LBj:) real(r8), intent(in) :: pnom_r(LBi:,LBj:) real(r8), intent(in) :: pnom_v(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: pmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) real(r8), intent(in) :: Kh(LBi:,LBj:) real(r8), intent(in) :: Kv(LBi:,LBj:,0:) real(r8), intent(inout) :: h(LBi:,LBj:) real(r8), intent(in) :: VnormR(LBi:,LBj:,:,:,:) real(r8), intent(in) :: VnormU(LBi:,LBj:,:,:) real(r8), intent(in) :: VnormV(LBi:,LBj:,:,:) real(r8), intent(in) :: HnormR(LBi:,LBj:,:) real(r8), intent(in) :: HnormU(LBi:,LBj:,:) real(r8), intent(in) :: HnormV(LBi:,LBj:,:) real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:) real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) real(r8), intent(out) :: Hz(LBi:,LBj:,:) real(r8), intent(out) :: z_r(LBi:,LBj:,:) real(r8), intent(out) :: z_w(LBi:,LBj:,0:) ! ! Local variable declarations. ! integer :: i, ib, ir, is, it, j, k, rec real(r8) :: cff real(r8) :: fac real(r8), dimension(LBi:UBi,LBj:UBj) :: work ! ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU integer :: Iend, IendB, IendP, IendR, IendT integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV integer :: Jend, JendB, JendP, JendR, JendT integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1 integer :: Iendp1, Iendp2, Iendp2i, Iendp3 integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1 integer :: Jendp1, Jendp2, Jendp2i, Jendp3 ! Istr =BOUNDS(ng) % Istr (tile) IstrB =BOUNDS(ng) % IstrB (tile) IstrM =BOUNDS(ng) % IstrM (tile) IstrP =BOUNDS(ng) % IstrP (tile) IstrR =BOUNDS(ng) % IstrR (tile) IstrT =BOUNDS(ng) % IstrT (tile) IstrU =BOUNDS(ng) % IstrU (tile) Iend =BOUNDS(ng) % Iend (tile) IendB =BOUNDS(ng) % IendB (tile) IendP =BOUNDS(ng) % IendP (tile) IendR =BOUNDS(ng) % IendR (tile) IendT =BOUNDS(ng) % IendT (tile) Jstr =BOUNDS(ng) % Jstr (tile) JstrB =BOUNDS(ng) % JstrB (tile) JstrM =BOUNDS(ng) % JstrM (tile) JstrP =BOUNDS(ng) % JstrP (tile) JstrR =BOUNDS(ng) % JstrR (tile) JstrT =BOUNDS(ng) % JstrT (tile) JstrV =BOUNDS(ng) % JstrV (tile) Jend =BOUNDS(ng) % Jend (tile) JendB =BOUNDS(ng) % JendB (tile) JendP =BOUNDS(ng) % JendP (tile) JendR =BOUNDS(ng) % JendR (tile) JendT =BOUNDS(ng) % JendT (tile) ! Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3 Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2 Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1 IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2 IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1 Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1 Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2 Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3 Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3 Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2 Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1 JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2 JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1 Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1 Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2 Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3 ! !----------------------------------------------------------------------- ! Determine error covariance normalization factors to use. !----------------------------------------------------------------------- ! IF (Lweak) THEN rec=2 ! weak constraint ELSE rec=1 ! strong constraint END IF ! !----------------------------------------------------------------------- ! Compute time invariant depths (use zero free-surface). !----------------------------------------------------------------------- ! DO i=LBi,UBi DO j=LBj,UBj work(i,j)=0.0_r8 END DO END DO CALL set_depth_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & h, & & work, & & Hz, z_r, z_w) ! !----------------------------------------------------------------------- ! Multiply tangent linear state by the inverse squared root of its ! associated area (2D) or volume (3D). !----------------------------------------------------------------------- ! ! Tangent linear free-surface. ! DO j=JstrT,JendT DO i=IstrT,IendT tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)/ & & SQRT(om_r(i,j)*on_r(i,j)) END DO END DO ! ! Tangent linear 2D momentum. ! DO j=JstrT,JendT DO i=IstrP,IendT tl_ubar(i,j,Linp)=tl_ubar(i,j,Linp)/ & & SQRT(om_u(i,j)*on_u(i,j)) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT tl_vbar(i,j,Linp)=tl_vbar(i,j,Linp)/ & & SQRT(om_v(i,j)*on_v(i,j)) END DO END DO CALL mp_exchange2d (ng, tile, iTLM, 3, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_zeta(:,:,Linp), & & tl_ubar(:,:,Linp), & & tl_vbar(:,:,Linp)) ! ! Tangent linear 3D momentum. ! DO j=JstrT,JendT DO i=IstrP,IendT cff=om_u(i,j)*on_u(i,j)*0.5_r8 DO k=1,N(ng) tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)/ & & SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) END DO END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT cff=om_v(i,j)*on_v(i,j)*0.5_r8 DO k=1,N(ng) tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)/ & & SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) END DO END DO END DO CALL mp_exchange3d (ng, tile, iTLM, 2, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_u(:,:,:,Linp), & & tl_v(:,:,:,Linp)) ! ! Tangent linear tracers. ! DO j=JstrT,JendT DO i=IstrT,IendT cff=om_r(i,j)*on_r(i,j) DO k=1,N(ng) fac=1.0_r8/SQRT(cff*Hz(i,j,k)) DO it=1,NT(ng) tl_t(i,j,k,Linp,it)=fac*tl_t(i,j,k,Linp,it) END DO END DO END DO END DO CALL mp_exchange4d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_t(:,:,:,Linp,:)) ! !----------------------------------------------------------------------- ! Initial conditions and model error convariance: Convolve tangent ! linear state vector with a generalized diffusion equation to filter ! solution with specified horizontal scales. Convert from minimization ! space (v-space) to model space. !----------------------------------------------------------------------- ! ! Tangent linear free-surface. ! CALL tl_conv_r2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,isFsur)/ifac, & & DTsizeH(rec,isFsur), & & Kh, & & pm, pn, pmon_u, pnom_v, & & rmask, umask, vmask, & & tl_zeta(:,:,Linp)) ! ! Tangent linear 2D momentum. ! CALL tl_conv_u2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,isUbar)/ifac, & & DTsizeH(rec,isUbar), & & Kh, & & pm, pn, pmon_r, pnom_p, & & umask, pmask, & & tl_ubar(:,:,Linp)) CALL tl_conv_v2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,isVbar)/ifac, & & DTsizeH(rec,isVbar), & & Kh, & & pm, pn, pmon_p, pnom_r, & & vmask, pmask, & & tl_vbar(:,:,Linp)) ! ! Tangent linear 3D momentum. ! CALL tl_conv_u3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,isUvel)/ifac, & & NVsteps(rec,isUvel)/ifac, & & DTsizeH(rec,isUvel), & & DTsizeV(rec,isUvel), & & Kh, Kv, & & pm, pn, & & pmon_r, pnom_p, & & umask, pmask, & & Hz, z_r, & & tl_u(:,:,:,Linp)) CALL tl_conv_v3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,isVvel)/ifac, & & NVsteps(rec,isVvel)/ifac, & & DTsizeH(rec,isVvel), & & DTsizeV(rec,isVvel), & & Kh, Kv, & & pm, pn, & & pmon_p, pnom_r, & & vmask, pmask, & & Hz, z_r, & & tl_v(:,:,:,Linp)) ! ! Tangent linear tracers. ! DO it=1,NT(ng) is=isTvar(it) CALL tl_conv_r3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(rec,is)/ifac, & & NVsteps(rec,is)/ifac, & & DTsizeH(rec,is), & & DTsizeV(rec,is), & & Kh, Kv, & & pm, pn, & & pmon_u, pnom_v, & & rmask, umask, vmask, & & Hz, z_r, & & tl_t(:,:,:,Linp,it)) END DO ! !----------------------------------------------------------------------- ! Multiply convolved tangent linear state by its corresponding ! normalization factor. !----------------------------------------------------------------------- ! ! Tangent linear free-surface. ! DO j=JstrT,JendT DO i=IstrT,IendT tl_zeta(i,j,Linp)=tl_zeta(i,j,Linp)*HnormR(i,j,rec) END DO END DO ! ! Tangent linear 2D momentum. ! DO j=JstrT,JendT DO i=IstrP,IendT tl_ubar(i,j,Linp)=tl_ubar(i,j,Linp)*HnormU(i,j,rec) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT tl_vbar(i,j,Linp)=tl_vbar(i,j,Linp)*HnormV(i,j,rec) END DO END DO CALL mp_exchange2d (ng, tile, iTLM, 3, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_zeta(:,:,Linp), & & tl_ubar(:,:,Linp), & & tl_vbar(:,:,Linp)) ! ! Tangent linear 3D momentum. ! DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT tl_u(i,j,k,Linp)=tl_u(i,j,k,Linp)*VnormU(i,j,k,rec) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT tl_v(i,j,k,Linp)=tl_v(i,j,k,Linp)*VnormV(i,j,k,rec) END DO END DO END DO CALL mp_exchange3d (ng, tile, iTLM, 2, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_u(:,:,:,Linp), & & tl_v(:,:,:,Linp)) ! ! Tangent linear tracers. ! DO it=1,NT(ng) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT tl_t(i,j,k,Linp,it)=tl_t(i,j,k,Linp,it)* & & VnormR(i,j,k,rec,it) END DO END DO END DO END DO CALL mp_exchange4d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & tl_t(:,:,:,Linp,:)) RETURN END SUBROUTINE tl_convolution_tile END MODULE tl_convolution_mod