SUBROUTINE ana_perturb (ng, tile, model) ! !! git $Id$ !! svn $Id: ana_perturb.h 1180 2023-07-13 02:42:10Z arango $ !!====================================================================== !! Copyright (c) 2002-2023 The ROMS/TOMS Group ! !! Licensed under a MIT/X style license ! !! See License_ROMS.md ! !======================================================================= ! ! ! This routine perturbs initial conditions for momentum and tracers ! ! type variables using analytical expressions. ! ! ! ! It is also used to perturb the tangent linear and adjoint models ! ! at specified state variable and spatial (i,j,k) point to verify ! ! the correctness of these algorithms. This is activated with the ! ! SANITY_CHECK CPP switch. ! ! ! ! If each interior point is perturbed at one time, the resulting ! ! tangent linear (T) and adjoint (A) M-by-N matrices yield: ! ! ! ! T - tranpose(A) = 0 within round off ! ! ! ! That is, their inner product give a symmetric matrix. Here, M is ! ! the number of state points and N is the number of perturbations. ! ! In realistic applications, it is awkward to perturb all interior ! ! points for each state variable. Alternatively, random check at a ! ! specified points is inexpensive. The standard input "User" array ! ! is used to specify such point: ! ! ! ! INT(user(1)) => state variable to perturb ! ! INT(user(2)) => I-index to perturb ! ! INT(user(3)) => J-index to perturb ! ! INT(user(4)) => K-index to perturb (3D state fields) ! ! ! !======================================================================= ! USE mod_param USE mod_ncparam #ifdef ADJUST_BOUNDARY USE mod_boundary #endif USE mod_ocean #if defined ADJUST_STFLUX || defined ADJUST_WSTRESS USE mod_forces #endif USE mod_stepping ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__ ! #include "tile.h" ! CALL ana_perturb_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), krhs(ng), knew(ng), & #ifdef SOLVE3D & nstp(ng), nrhs(ng), nnew(ng), & #endif #ifdef ADJUST_BOUNDARY # ifdef SOLVE3D & BOUNDARY(ng) % ad_t_obc, & & BOUNDARY(ng) % ad_u_obc, & & BOUNDARY(ng) % ad_v_obc, & # endif & BOUNDARY(ng) % ad_ubar_obc, & & BOUNDARY(ng) % ad_vbar_obc, & & BOUNDARY(ng) % ad_zeta_obc, & #endif #ifdef ADJUST_WSTRESS & FORCES(ng) % ad_ustr, & & FORCES(ng) % ad_vstr, & #endif #if defined ADJUST_STFLUX && defined SOLVE3D & FORCES(ng) % ad_tflux, & #endif #ifdef SOLVE3D & OCEAN(ng) % ad_t, & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & #endif & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & & OCEAN(ng) % ad_zeta, & #ifdef ADJUST_BOUNDARY # ifdef SOLVE3D & BOUNDARY(ng) % tl_t_obc, & & BOUNDARY(ng) % tl_u_obc, & & BOUNDARY(ng) % tl_v_obc, & # endif & BOUNDARY(ng) % tl_ubar_obc, & & BOUNDARY(ng) % tl_vbar_obc, & & BOUNDARY(ng) % tl_zeta_obc, & #endif #ifdef ADJUST_WSTRESS & FORCES(ng) % tl_ustr, & & FORCES(ng) % tl_vstr, & #endif #if defined ADJUST_STFLUX && defined SOLVE3D & FORCES(ng) % tl_tflux, & #endif #ifdef SOLVE3D & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & #endif & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & & OCEAN(ng) % tl_zeta) ! ! Set analytical header file name used. ! #ifdef DISTRIBUTE IF (Lanafile) THEN #else IF (Lanafile.and.(tile.eq.0)) THEN #endif ANANAME(19)=MyFile END IF ! RETURN END SUBROUTINE ana_perturb ! !*********************************************************************** SUBROUTINE ana_perturb_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & kstp, krhs, knew, & #ifdef SOLVE3D & nstp, nrhs, nnew, & #endif #ifdef ADJUST_BOUNDARY # ifdef SOLVE3D & ad_t_obc, ad_u_obc, ad_v_obc, & # endif & ad_ubar_obc, ad_vbar_obc, & & ad_zeta_obc, & #endif #ifdef ADJUST_WSTRESS & ad_ustr, ad_vstr, & #endif #if defined ADJUST_STFLUX && defined SOLVE3D & ad_tflux, & #endif #ifdef SOLVE3D & ad_t, ad_u, ad_v, & #endif & ad_ubar, ad_vbar, ad_zeta, & #ifdef ADJUST_BOUNDARY # ifdef SOLVE3D & tl_t_obc, tl_u_obc, tl_v_obc, & # endif & tl_ubar_obc, tl_vbar_obc, & & tl_zeta_obc, & #endif #ifdef ADJUST_WSTRESS & tl_ustr, tl_vstr, & #endif #if defined ADJUST_STFLUX && defined SOLVE3D & tl_tflux, & #endif #ifdef SOLVE3D & tl_t, tl_u, tl_v, & #endif & tl_ubar, tl_vbar, tl_zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_iounits USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: kstp, krhs, knew #ifdef SOLVE3D integer, intent(in) :: nstp, nrhs, nnew #endif ! #ifdef ASSUMED_SHAPE # ifdef ADJUST_BOUNDARY # ifdef SOLVE3D real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:) real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:) real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:) # endif real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:) real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:) real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:) # endif # ifdef ADJUST_WSTRESS real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:) # endif # if defined ADJUST_STFLUX && defined SOLVE3D real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:) # endif # ifdef SOLVE3D real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:) # endif real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:) # ifdef ADJUST_BOUNDARY # ifdef SOLVE3D real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:) real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:) real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:) # endif real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:) real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:) real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:) # endif # ifdef ADJUST_WSTRESS real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:) real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:) # endif # if defined ADJUST_STFLUX && defined SOLVE3D real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:) # endif # ifdef SOLVE3D real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:) # endif real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) #else # ifdef ADJUST_BOUNDARY # ifdef SOLVE3D real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, & & Nbrec(ng),2,NT(ng)) real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2) real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2) # endif real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2) real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2) real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2) # endif # ifdef ADJUST_WSTRESS real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2) real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2) # endif # if defined ADJUST_STFLUX && defined SOLVE3D real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, & & Nfrec(ng),2,NT(ng)) # endif # ifdef SOLVE3D real(r8), intent(inout) :: ad_t(LBi:UBI,LBj:UBj,N(ng),2,NT(ng)) real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2) # endif real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:) # ifdef ADJUST_BOUNDARY # ifdef SOLVE3D real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, & & Nbrec(ng),2,NT(ng)) real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2) real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2) # endif real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2) real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2) real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2) # endif # ifdef ADJUST_WSTRESS real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2) real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2) # endif # if defined ADJUST_STFLUX && defined SOLVE3D real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, & & Nfrec(ng),2,NT(ng)) # endif # ifdef SOLVE3D real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) # endif real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:) #endif ! ! Local variable declarations. ! #ifdef ADJUST_BOUNDARY logical :: Lperturb(4) ! #endif integer :: IperAD, JperAD, KperAD, ivarAD integer :: IperTL, JperTL, KperTL, ivarTL integer :: i, ib, ir, itrc, j, k #include "set_bounds.h" ! !----------------------------------------------------------------------- ! Set tangent and adjoint variable and random point to perturb. !----------------------------------------------------------------------- ! ivarTL=INT(user(1)) ivarAD=INT(user(2)) IperTL=INT(user(3)) IperAD=INT(user(4)) JperTL=INT(user(5)) JperAD=INT(user(6)) #ifdef SOLVE3D KperTL=INT(user(7)) KperAD=INT(user(8)) #endif IF (Master) THEN IF (TLmodel) THEN IF (ivarTL.eq.isUbar) THEN WRITE (stdout,10) 'tl_ubar perturbed at (i,j) = ', & & IperTL, JperTL ELSE IF (ivarTL.eq.isVbar) THEN WRITE (stdout,10) 'tl_vbar perturbed at (i,j) = ', & & IperTL, JperTL ELSE IF (ivarTL.eq.isFsur) THEN WRITE (stdout,10) 'tl_zeta perturbed at (i,j) = ', & & IperTL, JperTL #ifdef ADJUST_WSTRESS ELSE IF (ivarTL.eq.isUstr) THEN WRITE (stdout,10) 'tl_ustr perturbed at (i,j) = ', & & IperTL, JperTL ELSE IF (ivarTL.eq.isVstr) THEN WRITE (stdout,10) 'tl_vstr perturbed at (i,j) = ', & & IperTL, JperTL #endif #ifdef SOLVE3D ELSE IF (ivarTL.eq.isUvel) THEN WRITE (stdout,20) 'tl_u perturbed at (i,j,k) = ', & & IperTL, JperTL, KperTL ELSE IF (ivarTL.eq.isVvel) THEN WRITE (stdout,20) 'tl_v perturbed at (i,j,k) = ', & & IperTL, JperTL, KperTL #endif END IF #ifdef SOLVE3D DO itrc=1,NT(ng) IF (ivarTL.eq.isTvar(itrc)) THEN WRITE (stdout,30) 'tl_t perturbed at (i,j,k,itrc) = ', & & IperTL, JperTL, KperTL, itrc # ifdef ADJUST_STFLUX ELSE IF (ivarTL.eq.isTsur(itrc)) THEN WRITE (stdout,20) 'tl_tflux perturbed at (i,j,k,itrc) = ',& & IperTL, JperTL, KperTL, itrc # endif END IF END DO #endif #ifdef ADJUST_BOUNDARY IF (ivarTL.eq.isUbar) THEN WRITE (stdout,10) 'tl_ubar_obc (S/N) perturbed at (i) = ', & & IperTL WRITE (stdout,10) 'tl_ubar_obc (E/W) perturbed at (j) = ', & & JperTL ELSE IF (ivarTL.eq.isVbar) THEN WRITE (stdout,10) 'tl_vbar_obc (S/N) perturbed at (i) = ', & & IperTL WRITE (stdout,10) 'tl_vbar_obc (E/W) perturbed at (j) = ', & & JperTL ELSE IF (ivarTL.eq.isFsur) THEN WRITE (stdout,10) 'tl_zeta_obc (S/N) perturbed at (i) = ', & & IperTL WRITE (stdout,10) 'tl_zeta_obc (E/W) perturbed at (j) = ', & & JperTL # ifdef SOLVE3D ELSE IF (ivarTL.eq.isUvel) THEN WRITE (stdout,10) 'tl_u_obc (S/N) perturbed at (i,k) = ', & & IperTL, KperTL WRITE (stdout,10) 'tl_u_obc (E/W) perturbed at (j,k) = ', & & JperTL, KperTL ELSE IF (ivarTL.eq.isUvel) THEN WRITE (stdout,10) 'tl_u_obc (S/N) perturbed at (i,k) = ', & & IperTL, KperTL WRITE (stdout,10) 'tl_u_obc (E/W) perturbed at (j,k) = ', & & JperTL, KperTL # endif END IF # ifdef SOLVE3D DO itrc=1,NT(ng) IF (ivarTL.eq.isTvar(itrc)) THEN WRITE (stdout,20) 'tl_t_obc perturbed at (i,k,itrc) = ', & & IperTL, KperTL, itrc WRITE (stdout,20) 'tl_t_obc perturbed at (j,k,itrc) = ', & & JperTL, KperTL, itrc END IF END DO # endif #endif END IF IF (ADmodel) THEN IF (ivarAD.eq.isUbar) THEN WRITE (stdout,40) 'ad_ubar perturbed at (i,j) = ', & & IperAD, JperAD ELSE IF (ivarAD.eq.isVbar) THEN WRITE (stdout,40) 'ad_vbar perturbed at (i,j) = ', & & IperAD, JperAD ELSE IF (ivarAD.eq.isFsur) THEN WRITE (stdout,40) 'ad_zeta perturbed at (i,j) = ', & & IperAD, JperAD #ifdef ADJUST_WSTRESS ELSE IF (ivarAD.eq.isUstr) THEN WRITE (stdout,40) 'ad_ustr perturbed at (i,j,k) = ', & & IperAD, JperAD ELSE IF (ivarAD.eq.isVstr) THEN WRITE (stdout,40) 'ad_vstr perturbed at (i,j,k) = ', & & IperAD, JperAD #endif #ifdef SOLVE3D ELSE IF (ivarAD.eq.isUvel) THEN WRITE (stdout,50) 'ad_u perturbed at (i,j,k) = ', & & IperAD, JperAD, KperAD ELSE IF (ivarAD.eq.isVvel) THEN WRITE (stdout,50) 'ad_v perturbed at (i,j,k) = ', & & IperAD, JperAD, KperAD #endif END IF #ifdef SOLVE3D DO itrc=1,NT(ng) IF (ivarAD.eq.isTvar(itrc)) THEN WRITE (stdout,60) 'ad_t perturbed at (i,j,k,itrc) = ', & & IperAD, JperAD, KperAD, itrc # ifdef ADJUST_STFLUX ELSE IF (ivarAD.eq.isTsur(itrc)) THEN WRITE (stdout,50) 'ad_tflux perturbed at (i,j,k,itrc) = ',& & IperAD, JperAD, KperAD, itrc # endif END IF END DO #endif #ifdef ADJUST_BOUNDARY IF (ivarAD.eq.isUbar) THEN WRITE (stdout,40) 'ad_ubar_obc (S/N) perturbed at (i) = ', & & IperAD WRITE (stdout,40) 'ad_ubar_obc (E/W) perturbed at (j) = ', & & JperAD ELSE IF (ivarAD.eq.isVbar) THEN WRITE (stdout,40) 'ad_vbar_obc (S/N) perturbed at (i) = ', & & IperAD WRITE (stdout,40) 'ad_vbar_obc (E/W) perturbed at (j) = ', & & JperAD ELSE IF (ivarAD.eq.isFsur) THEN WRITE (stdout,40) 'ad_zeta_obc (S/N) perturbed at (i) = ', & & IperAD WRITE (stdout,40) 'ad_zeta_obc (E/W) perturbed at (j) = ', & & JperAD # ifdef SOLVE3D ELSE IF (ivarAD.eq.isUvel) THEN WRITE (stdout,40) 'ad_u_obc (S/N) perturbed at (i,k) = ', & & IperAD, KperAD WRITE (stdout,40) 'ad_u_obc (E/W) perturbed at (j,k) = ', & & JperAD, KperAD ELSE IF (ivarAD.eq.isUvel) THEN WRITE (stdout,40) 'ad_u_obc (S/N) perturbed at (i,k) = ', & & IperAD, KperAD WRITE (stdout,40) 'ad_u_obc (E/W) perturbed at (j,k) = ', & & JperAD, KperAD # endif END IF # ifdef SOLVE3D DO itrc=1,NT(ng) IF (ivarAD.eq.isTvar(itrc)) THEN WRITE (stdout,50) 'ad_t_obc perturbed at (i,k,itrc) = ', & & IperAD, KperAD, itrc WRITE (stdout,50) 'ad_t_obc perturbed at (j,k,itrc) = ', & & JperAD, KperAD, itrc END IF END DO # endif #endif END IF END IF ! !----------------------------------------------------------------------- ! Peturb initial conditions for 2D momentum (m/s) components. !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarTL.eq.isUbar).and. & & (i.eq.IperTL).and.(j.eq.JperTL)) THEN tl_ubar(i,j,kstp)=1.0_r8 ELSE tl_ubar(i,j,kstp)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isVbar).and. & & (i.eq.IperTL).and.(j.eq.JperTL)) THEN tl_vbar(i,j,kstp)=1.0_r8 ELSE tl_vbar(i,j,kstp)=0.0_r8 END IF END DO END DO END IF ! IF (ADmodel) THEN DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarAD.eq.isUbar).and. & & (i.eq.IperAD).and.(j.eq.JperAD)) THEN ad_ubar(i,j,knew)=1.0_r8 ELSE ad_ubar(i,j,knew)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isVbar).and. & & (i.eq.IperAD).and.(j.eq.JperAD)) THEN ad_vbar(i,j,knew)=1.0_r8 ELSE ad_vbar(i,j,knew)=0.0_r8 END IF END DO END DO END IF #ifdef ADJUST_WSTRESS ! !----------------------------------------------------------------------- ! Peturb initial conditions for surface momentum stress (UNIT ????). !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO ir=1,Nfrec(ng) DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarTL.eq.isUstr).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (ir.eq.KperTL)) THEN tl_ustr(i,j,ir,kstp)=1.0_r8 ELSE tl_ustr(i,j,ir,kstp)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isVstr).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (ir.eq.KperTL)) THEN tl_vstr(i,j,ir,kstp)=1.0_r8 ELSE tl_vstr(i,j,ir,kstp)=0.0_r8 END IF END DO END DO END DO END IF ! IF (ADmodel) THEN DO ir=1,Nfrec(ng) DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarAD.eq.isUstr).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (ir.eq.KperAD)) THEN ad_ustr(i,j,ir,knew)=1.0_r8 ELSE ad_ustr(i,j,ir,knew)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isVstr).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (ir.eq.KperAD)) THEN ad_vstr(i,j,ir,knew)=1.0_r8 ELSE ad_vstr(i,j,ir,knew)=0.0_r8 END IF END DO END DO END DO END IF #endif ! !----------------------------------------------------------------------- ! Perturb initial conditions for free-surface (m). !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isFsur).and. & & (i.eq.IperTL).and.(j.eq.JperTL)) THEN tl_zeta(i,j,kstp)=1.0_r8 ELSE tl_zeta(i,j,kstp)=0.0_r8 END IF END DO END DO END IF ! IF (ADmodel) THEN DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isFsur).and. & & (i.eq.IperAD).and.(j.eq.JperAD)) THEN ad_zeta(i,j,knew)=1.0_r8 ELSE ad_zeta(i,j,knew)=0.0_r8 END IF END DO END DO END IF #ifdef SOLVE3D ! !----------------------------------------------------------------------- ! Initial conditions for 3D momentum components (m/s). !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarTL.eq.isUvel).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (k.eq.KperTL)) THEN tl_u(i,j,k,nstp)=1.0_r8 ELSE tl_u(i,j,k,nstp)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isVvel).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (k.eq.KperTL)) THEN tl_v(i,j,k,nstp)=1.0_r8 ELSE tl_v(i,j,k,nstp)=0.0_r8 END IF END DO END DO END DO END IF ! IF (ADmodel) THEN DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT IF ((ivarAD.eq.isUvel).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (k.eq.KperAD)) THEN ad_u(i,j,k,nstp)=1.0_r8 ELSE ad_u(i,j,k,nstp)=0.0_r8 END IF END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isVvel).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (k.eq.KperAD)) THEN ad_v(i,j,k,nstp)=1.0_r8 ELSE ad_v(i,j,k,nstp)=0.0_r8 END IF END DO END DO END DO END IF ! !----------------------------------------------------------------------- ! Perturb initial conditions for tracer type variables. !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO itrc=1,NT(ng) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isTvar(itrc)).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (k.eq.KperTL)) THEN tl_t(i,j,k,nstp,itrc)=1.0_r8 ELSE tl_t(i,j,k,nstp,itrc)=0.0_r8 END IF END DO END DO END DO END DO END IF ! IF (ADmodel) THEN DO itrc=1,NT(ng) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isTvar(itrc)).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (k.eq.KperAD)) THEN ad_t(i,j,k,nstp,itrc)=1.0_r8 ELSE ad_t(i,j,k,nstp,itrc)=0.0_r8 END IF END DO END DO END DO END DO END IF # ifdef ADJUST_STFLUX ! !----------------------------------------------------------------------- ! Perturb initial conditions for surface tracer flux. !----------------------------------------------------------------------- ! IF (TLmodel) THEN DO itrc=1,NT(ng) DO ir=1,Nfrec(ng) DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarTL.eq.isTsur(itrc)).and. & & (i.eq.IperTL).and.(j.eq.JperTL).and. & & (ir.eq.KperTL)) THEN tl_tflux(i,j,ir,nstp,itrc)=1.0_r8 ELSE tl_tflux(i,j,ir,nstp,itrc)=0.0_r8 END IF END DO END DO END DO END DO END IF ! IF (ADmodel) THEN DO itrc=1,NT(ng) DO ir=1,Nfrec(ng) DO j=JstrT,JendT DO i=IstrT,IendT IF ((ivarAD.eq.isTsur(itrc)).and. & & (i.eq.IperAD).and.(j.eq.JperAD).and. & & (ir.eq.KperAD)) THEN ad_tflux(i,j,ir,nstp,itrc)=1.0_r8 ELSE ad_tflux(i,j,ir,nstp,itrc)=0.0_r8 END IF END DO END DO END DO END DO END IF # endif # ifdef ADJUST_BOUNDARY ! !----------------------------------------------------------------------- ! Perturb open boundary conditions. !----------------------------------------------------------------------- ! Lperturb(iwest )=DOMAIN(ng)%Western_Edge (tile) Lperturb(ieast )=DOMAIN(ng)%Eastern_Edge (tile) Lperturb(isouth)=DOMAIN(ng)%Southern_Edge(tile) Lperturb(inorth)=DOMAIN(ng)%Northern_Edge(tile) DO ir=1,Nbrec(ng) DO ib=1,4 ! ! Perturb free-surface open boundaries. ! IF (Lperturb(ib).and.Lobc(ib,isFsur,ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isFsur)) THEN DO j=Jstr,Jend IF (j.eq.JperTL) THEN tl_zeta_obc(j,ib,ir,1)=1.0_r8 ELSE tl_zeta_obc(j,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isFsur)) THEN DO j=Jstr,Jend IF (j.eq.JperAD) THEN ad_zeta_obc(j,ib,ir,1)=1.0_r8 ELSE ad_zeta_obc(j,ib,ir,1)=0.0_r8 END IF END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isFsur)) THEN DO i=Istr,Iend IF (i.eq.IperTL) THEN tl_zeta_obc(i,ib,ir,1)=1.0_r8 ELSE tl_zeta_obc(i,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isFsur)) THEN DO i=Istr,Iend IF (i.eq.IperAD) THEN ad_zeta_obc(i,ib,ir,1)=1.0_r8 ELSE ad_zeta_obc(i,ib,ir,1)=0.0_r8 END IF END DO END IF END IF END IF ! ! Perturb 2D U-momentum open boundaries. ! IF (Lperturb(ib).and.Lobc(ib,isUbar,ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isUbar)) THEN DO j=Jstr,Jend IF (j.eq.JperTL) THEN tl_ubar_obc(j,ib,ir,1)=1.0_r8 ELSE tl_ubar_obc(j,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isUbar)) THEN DO j=Jstr,Jend IF (j.eq.JperAD) THEN ad_ubar_obc(j,ib,ir,1)=1.0_r8 ELSE ad_ubar_obc(j,ib,ir,1)=0.0_r8 END IF END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isUbar)) THEN DO i=IstrU,Iend IF (i.eq.IperTL) THEN tl_ubar_obc(i,ib,ir,1)=1.0_r8 ELSE tl_ubar_obc(i,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isUbar)) THEN DO i=IstrU,Iend IF (i.eq.IperAD) THEN ad_ubar_obc(i,ib,ir,1)=1.0_r8 ELSE ad_ubar_obc(i,ib,ir,1)=0.0_r8 END IF END DO END IF END IF END IF ! ! Perturb 2D V-momentum open boundaries. ! IF (Lperturb(ib).and.Lobc(ib,isVbar,ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isVbar)) THEN DO j=JstrV,Jend IF (j.eq.JperTL) THEN tl_vbar_obc(j,ib,ir,1)=1.0_r8 ELSE tl_vbar_obc(j,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isVbar)) THEN DO j=JstrV,Jend IF (j.eq.JperAD) THEN ad_vbar_obc(j,ib,ir,1)=1.0_r8 ELSE ad_vbar_obc(j,ib,ir,1)=0.0_r8 END IF END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isVbar)) THEN DO i=Istr,Iend IF (i.eq.IperTL) THEN tl_vbar_obc(i,ib,ir,1)=1.0_r8 ELSE tl_vbar_obc(i,ib,ir,1)=0.0_r8 END IF END DO ELSE IF (ADmodel.and.(ivarAD.eq.isVbar)) THEN DO i=Istr,Iend IF (i.eq.IperAD) THEN ad_vbar_obc(i,ib,ir,1)=1.0_r8 ELSE ad_vbar_obc(i,ib,ir,1)=0.0_r8 END IF END DO END IF END IF END IF # ifdef SOLVE3D ! ! Perturb 3D U-momentum open boundaries. ! IF (Lperturb(ib).and.Lobc(ib,isUvel,ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isUvel)) THEN DO k=1,N(ng) DO j=Jstr,Jend IF ((j.eq.JperTL).and.(k.eq.KperTL)) THEN tl_u_obc(j,k,ib,ir,1)=1.0_r8 ELSE tl_u_obc(j,k,ib,ir,1)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isUvel)) THEN DO k=1,N(ng) DO j=Jstr,Jend IF ((j.eq.JperAD).and.(k.eq.KperAD)) THEN ad_u_obc(j,k,ib,ir,1)=1.0_r8 ELSE ad_u_obc(j,k,ib,ir,1)=0.0_r8 END IF END DO END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isUvel)) THEN DO k=1,N(ng) DO i=IstrU,Iend IF ((i.eq.IperTL).and.(k.eq.KperTL)) THEN tl_u_obc(i,k,ib,ir,1)=1.0_r8 ELSE tl_u_obc(i,k,ib,ir,1)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isUvel)) THEN DO k=1,N(ng) DO i=IstrU,Iend IF ((i.eq.IperAD).and.(k.eq.KperAD)) THEN ad_u_obc(i,k,ib,ir,1)=1.0_r8 ELSE ad_u_obc(i,k,ib,ir,1)=0.0_r8 END IF END DO END DO END IF END IF END IF ! ! Perturb 3D V-momentum open boundaries. ! IF (Lperturb(ib).and.Lobc(ib,isVvel,ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isVvel)) THEN DO k=1,N(ng) DO j=JstrV,Jend IF ((j.eq.JperTL).and.(k.eq.KperTL)) THEN tl_v_obc(j,k,ib,ir,1)=1.0_r8 ELSE tl_v_obc(j,k,ib,ir,1)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isVvel)) THEN DO k=1,N(ng) DO j=JstrV,Jend IF ((j.eq.JperAD).and.(k.eq.KperAD)) THEN ad_v_obc(j,k,ib,ir,1)=1.0_r8 ELSE ad_v_obc(j,k,ib,ir,1)=0.0_r8 END IF END DO END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isVvel)) THEN DO k=1,N(ng) DO i=Istr,Iend IF ((i.eq.IperTL).and.(k.eq.KperTL)) THEN tl_v_obc(i,k,ib,ir,1)=1.0_r8 ELSE tl_v_obc(i,k,ib,ir,1)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isVvel)) THEN DO k=1,N(ng) DO i=Istr,Iend IF ((i.eq.IperAD).and.(k.eq.KperAD)) THEN ad_v_obc(i,k,ib,ir,1)=1.0_r8 ELSE ad_v_obc(i,k,ib,ir,1)=0.0_r8 END IF END DO END DO END IF END IF END IF ! ! Perturb tracers open boundaries. ! DO itrc=1,NT(ng) IF (Lperturb(ib).and.Lobc(ib,isTvar(itrc),ng)) THEN IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN IF (TLmodel.and.(ivarTL.eq.isTvar(itrc))) THEN DO k=1,N(ng) DO j=Jstr,Jend IF ((j.eq.JperTL).and.(k.eq.KperTL)) THEN tl_t_obc(j,k,ib,ir,1,itrc)=1.0_r8 ELSE tl_t_obc(j,k,ib,ir,1,itrc)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isTvar(itrc))) THEN DO k=1,N(ng) DO j=Jstr,Jend IF ((j.eq.JperAD).and.(k.eq.KperAD)) THEN ad_t_obc(j,k,ib,ir,1,itrc)=1.0_r8 ELSE ad_t_obc(j,k,ib,ir,1,itrc)=0.0_r8 END IF END DO END DO END IF ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN IF (TLmodel.and.(ivarTL.eq.isTvar(itrc))) THEN DO k=1,N(ng) DO i=Istr,Iend IF ((i.eq.IperTL).and.(k.eq.KperTL)) THEN tl_t_obc(i,k,ib,ir,1,itrc)=1.0_r8 ELSE tl_t_obc(i,k,ib,ir,1,itrc)=0.0_r8 END IF END DO END DO ELSE IF (ADmodel.and.(ivarAD.eq.isTvar(itrc))) THEN DO k=1,N(ng) DO i=Istr,Iend IF ((i.eq.IperAD).and.(k.eq.KperAD)) THEN ad_t_obc(i,k,ib,ir,1,itrc)=1.0_r8 ELSE ad_t_obc(i,k,ib,ir,1,itrc)=0.0_r8 END IF END DO END DO END IF END IF END IF END DO # endif END DO END DO # endif #endif ! 10 FORMAT (/,' ANA_PERTURB - Tangent ', a, 2i4,/) #ifdef SOLVE3D 20 FORMAT (/,' ANA_PERTURB - Tangent ', a, 3i4,/) 30 FORMAT (/,' ANA_PERTURB - Tangent ', a, 4i4,/) #endif 40 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 2i4,/) #ifdef SOLVE3D 50 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 3i4,/) 60 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 4i4,/) #endif ! RETURN END SUBROUTINE ana_perturb_tile