#include "cppdefs.h" #if defined OPT_PERTURBATION || defined FORCING_SV # define ENERGYNORM_SCALE #endif #if defined STOCHASTIC_OPT # ifndef HESSIAN_SO # define ENERGYNORM_SCALE # endif #endif #ifdef FULL_GRID # define IR_RANGE IstrT,IendT # define IU_RANGE IstrP,IendT # define JR_RANGE JstrT,JendT # define JV_RANGE JstrP,JendT #else # define IR_RANGE Istr,Iend # define IU_RANGE IstrU,Iend # define JR_RANGE Jstr,Jend # define JV_RANGE JstrV,Jend #endif MODULE packing_mod #ifdef PROPAGATOR ! !git $Id$ !svn $Id: packing.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 ! !======================================================================= ! ! ! These routines pack and unpack model state varaibles into/from a ! ! single vector to interface with ARPACKs Arnoldi Method for the ! ! computation Ritz eigenfunctions. ! ! ! !======================================================================= ! implicit none ! PRIVATE PUBLIC :: c_norm2 PUBLIC :: r_norm2 # ifdef ADJOINT # ifdef STOCHASTIC_OPT # ifdef STOCH_OPT_WHITE PUBLIC :: ad_so_pack # else PUBLIC :: ad_so_pack_red # endif # else PUBLIC :: ad_pack # endif PUBLIC :: ad_unpack # endif # ifdef TANGENT PUBLIC :: tl_pack PUBLIC :: tl_unpack # endif # ifdef SO_SEMI # ifdef SO_SEMI_WHITE PUBLIC :: so_semi_white # else PUBLIC :: so_semi_red # endif # endif ! CONTAINS ! SUBROUTINE c_norm2 (ng, model, Mstr, Mend, & & EvalueR, EvalueI, EvectorR, EvectorI, & & state, norm2) ! !======================================================================= ! ! ! This function computes the Euclidean norm between the propagator ! ! real/imaginary Ritz eigenvalue (EvalueR, EvalueI) and eigenvector ! ! (EvectorR, EvectorI) with state vector (state): ! ! ! ! norm2 = Euclidean NORM (state(:) + EvalueR * EvectorR(:) + ! ! EvalueI * EvectorI(:)) ! ! ! ! WARNING: This function is only intended for serial or distributed ! ! memory applications. There is not tiled partitions. All ! ! quantities are vectors. It replaces the calls to "daxpy" ! ! and "dnrm2" from the BLAS library. This "legacy" library ! ! gives different results when called inside modules and ! ! the input arguments are pointers (specially using ifort). ! ! ! !======================================================================= ! USE mod_param USE mod_parallel # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_reduce # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, model integer, intent(in) :: Mstr, Mend real(r8), intent(in) :: EvalueR real(r8), intent(in) :: EvalueI # ifdef ASSUMED_SHAPE real(r8), intent(in) :: EvectorR(Mstr:) real(r8), intent(in) :: EvectorI(Mstr:) real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: EvectorR(Mstr:Mend) real(r8), intent(in) :: EvectorI(Mstr:Mend) real(r8), intent(in) :: state(Mstr:Mend) # endif real(r8), intent(out) :: norm2 ! ! Local variable declarations. ! integer :: NSUB, is real(r8) :: cff, my_norm2 # ifdef DISTRIBUTE character (len=3) :: op_handle # endif ! !----------------------------------------------------------------------- ! Compute the Euclidean norm of: state(:) + Rvalue * Rvector(:) !----------------------------------------------------------------------- ! ! Accumulate squared sum. ! my_norm2=0.0_r8 DO is=Mstr,Mend cff=state(is)+EvalueR*EvectorR(is)+ & & EvalueI*EvectorI(is) my_norm2=my_norm2+cff*cff END DO ! ! Take sum squared-root: perform global reduction. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else NSUB=NtileX(ng)*NtileE(ng) ! tiled application # endif !$OMP CRITICAL (C_NORM) IF (tile_count.eq.0) THEN norm2=my_norm2 ELSE norm2=norm2+my_norm2 END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE op_handle='SUM' CALL mp_reduce (ng, model, 1, norm2, op_handle) # endif END IF !$OMP END CRITICAL (C_NORM) norm2=SQRT(norm2) ! RETURN END SUBROUTINE c_norm2 ! # if defined HESSIAN_FSV || defined HESSIAN_SO || defined HESSIAN_SV SUBROUTINE r_norm2 (ng, model, Mstr, Mend, & & Evalue, Evector, state, norm2) ! !======================================================================= ! ! ! This function computes the Euclidean norm between the propagator ! ! real Ritz eigenvalue (Evalue) and eigenvector (Evector) with the ! ! state vector (state): ! ! ! ! norm2 = Euclidean NORM (state(:) + Evalue * Evector(:)) ! ! ! ! WARNING: The norm is computed by the master thread and broadcasted ! ! to all the nodes in the group. It is used when the state ! ! vector is not partitioned between all nodes. ! ! ! !======================================================================= ! USE mod_param USE mod_parallel # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_bcastf # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, model integer, intent(in) :: Mstr, Mend real(r8), intent(in) :: Evalue # ifdef ASSUMED_SHAPE real(r8), intent(in) :: Evector(Mstr:) real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: Evector(Mstr:Mend) real(r8), intent(in) :: state(Mstr:Mend) # endif real(r8), intent(out) :: norm2 ! ! Local variable declarations. ! integer :: NSUB, is real(r8) :: cff, my_norm2 ! !----------------------------------------------------------------------- ! Compute the Euclidean norm of: state(:) + Rvalue * Rvector(:) !----------------------------------------------------------------------- ! ! Accumulate squared sum. ! IF (Master) THEN my_norm2=0.0_r8 DO is=Mstr,Mend cff=state(is)+Evalue*Evector(is) my_norm2=my_norm2+cff*cff END DO norm2=SQRT(my_norm2) END IF # ifdef DISTRIBUTE CALL mp_bcastf (ng, model, norm2) # endif ! RETURN END SUBROUTINE r_norm2 # else SUBROUTINE r_norm2 (ng, model, Mstr, Mend, & & Evalue, Evector, state, norm2) ! !======================================================================= ! ! ! This function computes the Euclidean norm between the propagator ! ! real Ritz eigenvalue (Evalue) and eigenvector (Evector) with the ! ! state vector (state): ! ! ! ! norm2 = Euclidean NORM (state(:) + Evalue * Evector(:)) ! ! ! ! WARNING: This function is only intended for serial or distributed ! ! memory applications. There is not tiled partitions. All ! ! quantities are vectors. It replaces the calls to "daxpy" ! ! and "dnrm2" from the BLAS library. This "legacy" library ! ! gives different results when called inside modules and ! ! the input arguments are pointers (specially using ifort). ! ! ! !======================================================================= ! USE mod_param USE mod_parallel # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_reduce # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, model integer, intent(in) :: Mstr, Mend real(r8), intent(in) :: Evalue # ifdef ASSUMED_SHAPE real(r8), intent(in) :: Evector(Mstr:) real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: Evector(Mstr:Mend) real(r8), intent(in) :: state(Mstr:Mend) # endif real(r8), intent(out) :: norm2 ! ! Local variable declarations. ! integer :: NSUB, is real(r8) :: cff, my_norm2 # ifdef DISTRIBUTE character (len=3) :: op_handle # endif ! !----------------------------------------------------------------------- ! Compute the Euclidean norm of: state(:) + Rvalue * Rvector(:) !----------------------------------------------------------------------- ! ! Accumulate squared sum. ! my_norm2=0.0_r8 DO is=Mstr,Mend cff=state(is)+Evalue*Evector(is) my_norm2=my_norm2+cff*cff END DO ! ! Take sum squared-root: perform global reduction. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else NSUB=NtileX(ng)*NtileE(ng) ! tiled application # endif !$OMP CRITICAL (R_NORM) IF (tile_count.eq.0) THEN norm2=my_norm2 ELSE norm2=norm2+my_norm2 END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE op_handle='SUM' CALL mp_reduce (ng, model, 1, norm2, op_handle) # endif END IF !$OMP END CRITICAL (R_NORM) norm2=SQRT(norm2) ! RETURN END SUBROUTINE r_norm2 # endif # if defined ADJOINT && defined FORCING_SV ! SUBROUTINE ad_pack (ng, tile, Mstr, Mend, ad_state) ! !======================================================================= ! ! ! This routine packs the adjoint variables into the state vector. ! ! The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_forces USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_pack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, ad_state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % f_t, & & OCEAN(ng) % f_u, & & OCEAN(ng) % f_v, & & FORCES(ng) % ad_stflx, & # endif & OCEAN(ng) % f_ubar, & & OCEAN(ng) % f_vbar, & & OCEAN(ng) % f_zeta, & & FORCES(ng) % ad_sustr, & & FORCES(ng) % ad_svstr) # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Scatter (global to threaded) adjoint state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iADM, Mstr, Mend, Mstate(ng), & & Swork, ad_state) # endif ! RETURN END SUBROUTINE ad_pack ! !*********************************************************************** SUBROUTINE ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, ad_state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & f_t, f_u, f_v, ad_stflx, & # endif & f_ubar, f_vbar, & & f_zeta, ad_sustr, ad_svstr) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_scalars USE mod_ocean ! # ifdef FORCING_SV USE ad_exchange_2d_mod # ifdef SOLVE3D USE ad_exchange_3d_mod # endif # ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : ad_mp_exchange2d # ifdef SOLVE3D USE mp_exchange_mod, ONLY : ad_mp_exchange3d, ad_mp_exchange4d # endif # endif # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: f_t(LBi:,LBj:,:,:) real(r8), intent(inout) :: f_u(LBi:,LBj:,:) real(r8), intent(inout) :: f_v(LBi:,LBj:,:) real(r8), intent(inout) :: ad_stflx(LBi:,LBj:,:) # endif real(r8), intent(inout) :: f_ubar(LBi:,LBj:) real(r8), intent(inout) :: f_vbar(LBi:,LBj:) real(r8), intent(inout) :: f_zeta(LBi:,LBj:) real(r8), intent(inout) :: ad_sustr(LBi:,LBj:) real(r8), intent(inout) :: ad_svstr(LBi:,LBj:) real(r8), intent(out) :: ad_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: f_t(LBi:UBi,LBj:UBj,N(ng),NT(ng)) real(r8), intent(inout) :: f_u(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: f_v(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: ad_stflx(LBi:UBi,LBj:UBj,NT(ng)) # endif real(r8), intent(inout) :: f_ubar(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: f_vbar(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: f_zeta(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: ad_sustr(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: ad_svstr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, is, itrc, j, k # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif real(r8), parameter :: Aspv = 0.0_r8 real(r8) :: cff, scale # ifdef SOLVE3D real(r8) :: cff1, cff2 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC # endif # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize adjoint state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend ad_state(is)=Aspv END DO # endif # ifdef FORCING_SV ! ! Impose adjoint periodic boundary conditions as appropriate. ! # ifdef DISTRIBUTE CALL ad_mp_exchange2d (ng, tile, iADM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & f_zeta) # ifndef SOLVE3D CALL ad_mp_exchange2d (ng, tile, iADM, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & f_ubar, f_vbar) # else CALL ad_mp_exchange3d (ng, tile, iADM, 2, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), f_u, f_v) CALL ad_mp_exchange4d (ng, tile, iADM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), f_t) # endif # endif ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL ad_exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_zeta) # ifndef SOLVE3D CALL ad_exchange_u2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_ubar) CALL ad_exchange_v2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_vbar) # else CALL ad_exchange_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), f_u) CALL ad_exchange_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), f_v) DO itrc=1,NT(ng) CALL ad_exchange_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & f_t(:,:,:,itrc)) END DO # endif END IF # endif ! !----------------------------------------------------------------------- ! Load adjoint STATE variables into full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Load adjoint of free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_state(is)=scale*f_zeta(i,j) f_zeta(i,j)=0.0_r8 ELSE f_zeta(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_state(is)=scale*f_zeta(i,j) f_zeta(i,j)=0.0_r8 # endif END DO END DO END IF # ifndef SOLVE3D ! ! Load adjoint of 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_state(is)=scale*f_ubar(i,j) f_ubar(i,j)=0.0_r8 ELSE f_ubar(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_state(is)=scale*f_ubar(i,j) f_ubar(i,j)=0.0_r8 # endif END DO END DO END IF ! ! Load adjoint of 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_state(is)=scale*f_vbar(i,j) f_vbar(i,j)=0.0_r8 ELSE f_vbar(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_state(is)=scale*f_vbar(i,j) f_vbar(i,j)=0.0_r8 # endif END DO END DO END IF # else ! ! Load adjoint of 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN ! ! Compute the adjoint forcing for tl_ubar based on f_u. ! DO j=JR_RANGE DO i=IU_RANGE DC(i,0)=0.0_r8 CF(i,0)=0.0_r8 END DO DO k=1,N(ng) DO i=IU_RANGE DC(i,k)=0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k)) DC(i,0)=DC(i,0)+DC(i,k) END DO END DO DO i=IU_RANGE cff2=f_ubar(i,j) f_ubar(i,j)=0.0_r8 # ifdef MASKING cff2=cff2*umask(i,j) # endif cff1=1.0_r8/DC(i,0) CF(i,0)=cff2*cff1 cff2=0.0_r8 END DO DO k=1,N(ng) DO i=IU_RANGE f_u(i,j,k)=f_u(i,j,k)+DC(i,k)*CF(i,0) END DO END DO DO i=IU_RANGE CF(i,0)=0.0_r8 END DO END DO # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+iadd ad_state(is)=scale*f_u(i,j,k) f_u(i,j,k)=0.0_r8 ELSE f_u(i,j,k)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*f_u(i,j,k) f_u(i,j,k)=0.0_r8 # endif END DO END DO END DO END IF ! ! Load adjoint of 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN ! ! Compute the adjoint forcing for tl_vbar based on f_v. ! DO j=JV_RANGE IF (j.ge.JstrM) THEN DO i=IR_RANGE DC(i,0)=0.0_r8 CF(i,0)=0.0_r8 END DO DO k=1,N(ng) DO i=IR_RANGE DC(i,k)=0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k)) DC(i,0)=DC(i,0)+DC(i,k) END DO END DO DO i=IR_RANGE cff2=f_vbar(i,j) f_vbar(i,j)=0.0_r8 # ifdef MASKING cff2=cff2*vmask(i,j) # endif cff1=1.0_r8/DC(i,0) CF(i,0)=cff2*cff1 cff2=0.0_r8 END DO DO k=1,N(ng) DO i=IR_RANGE f_v(i,j,k)=f_v(i,j,k)+DC(i,k)*CF(i,0) END DO END DO DO i=IR_RANGE CF(i,0)=0.0_r8 END DO END IF END DO # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+iadd ad_state(is)=scale*f_v(i,j,k) f_v(i,j,k)=0.0_r8 ELSE f_v(i,j,k)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*f_v(i,j,k) f_v(i,j,k)=0.0_r8 # endif END DO END DO END DO END IF ! ! Load adjoint of tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+iadd ad_state(is)=scale*f_t(i,j,k,itrc) f_t(i,j,k,itrc)=0.0_r8 ELSE f_t(i,j,k,itrc)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*f_t(i,j,k,itrc) f_t(i,j,k,itrc)=0.0_r8 # endif END DO END DO END DO END IF END DO # endif ! ! Load adjoint of surface U-stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) ad_state(is)=scale*ad_sustr(i,j) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) ad_state(is)=scale*ad_sustr(i,j) # endif END DO END DO END IF ! ! Load adjoint of surface V-stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=1.0_r8 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) ad_state(is)=scale*ad_svstr(i,j) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) ad_state(is)=scale*ad_svstr(i,j) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Load adjoint of surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) ad_state(is)=scale*ad_stflx(i,j,itrc) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) ad_state(is)=scale*ad_stflx(i,j,itrc) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE ad_pack_tile # elif defined SO_SEMI ! SUBROUTINE ad_pack (ng, tile, Mstr, Mend, ad_state) ! !======================================================================= ! ! ! This routine packs the adjoint variables into the state vector. ! ! The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_pack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork) # else & Mstr, Mend, ad_state) # endif # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Scatter (global to threaded) adjoint state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iADM, Mstr, Mend, Mstate(ng), & & Swork, ad_state) # endif ! RETURN END SUBROUTINE ad_pack ! !*********************************************************************** SUBROUTINE ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & Mstr, Mend, ad_state) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_scalars USE mod_ocean ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(out) :: ad_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, is, itrc, j, k integer, dimension(7+2*NT(ng)) :: offset real(r8), parameter :: Aspv = 0.0_r8 real(r8) :: cff, scale # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize adjoint state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend ad_state(is)=Aspv END DO # endif ! !----------------------------------------------------------------------- ! Load adjoint STATE and FORCING variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Load adjoint of free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_state(is)=scale*OCEAN(ng)%ad_zeta(i,j,kstp) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_state(is)=scale*OCEAN(ng)%ad_zeta(i,j,kstp) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Load adjoint of 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_state(is)=scale*OCEAN(ng)%ad_ubar(i,j,kstp) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_state(is)=scale*OCEAN(ng)%ad_ubar(i,j,kstp) # endif END DO END DO END IF ! ! Load adjoint of 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_state(is)=scale*OCEAN(ng)%ad_vbar(i,j,kstp) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_state(is)=scale*OCEAN(ng)%ad_vbar(i,j,kstp) # endif END DO END DO END IF # else ! ! Load adjoint of 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+iadd ad_state(is)=scale*OCEAN(ng)%ad_u(i,j,k,nstp) END IF # else is=(i-Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*OCEAN(ng)%ad_u(i,j,k,nstp) # endif END DO END DO END DO END IF ! ! Load adjoint of 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+iadd ad_state(is)=scale*OCEAN(ng)%ad_v(i,j,k,nstp) END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*OCEAN(ng)%ad_v(i,j,k,nstp) # endif END DO END DO END DO END IF ! ! Load adjoint of tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+iadd ad_state(is)=scale*OCEAN(ng)%ad_t(i,j,k,nstp,itrc) END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*OCEAN(ng)%ad_t(i,j,k,nstp,itrc) # endif END DO END DO END DO END IF END DO # endif ! ! Load adjoint of surface U-stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) ad_state(is)=scale*FORCES(ng)%ad_sustr(i,j) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) ad_state(is)=scale*FORCES(ng)%ad_sustr(i,j) # endif END DO END DO END IF ! ! Load adjoint of surface V-stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=1.0_r8 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) ad_state(is)=scale*FORCES(ng)%ad_svstr(i,j) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) ad_state(is)=scale*FORCES(ng)%ad_svstr(i,j) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Load adjoint of surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) ad_state(is)=scale*FORCES(ng)%ad_stflx(i,j,itrc) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) ad_state(is)=scale*FORCES(ng)%ad_stflx(i,j,itrc) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE ad_pack_tile # elif defined STOCHASTIC_OPT # ifdef STOCH_OPT_WHITE ! SUBROUTINE ad_so_pack (ng, tile, Mstr, Mend, IntTrap, ad_state) ! !======================================================================= ! ! ! This routine packs the adjoint variables into the state vector. ! ! The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_forces USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend integer, intent(in) :: IntTrap # ifdef ASSUMED_SHAPE real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_so_pack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_so_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif & IntTrap, & # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, ad_state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % ad_t, & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & & FORCES(ng) % ad_stflx, & # else & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & # endif & OCEAN(ng) % ad_zeta, & & FORCES(ng) % ad_sustr, & & FORCES(ng) % ad_svstr) # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Scatter (global to threaded) adjoint state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iADM, Mstr, Mend, Mstate(ng), & & Swork, ad_state) # endif ! RETURN END SUBROUTINE ad_so_pack ! !*********************************************************************** SUBROUTINE ad_so_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & IntTrap, & & Mstr, Mend, ad_state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & ad_t, ad_u, ad_v, ad_stflx, & # else & ad_ubar, ad_vbar, & # endif & ad_zeta, ad_sustr, ad_svstr) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_scalars USE mod_ocean USE mod_storage ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp integer, intent(in) :: IntTrap # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_stflx(LBi:,LBj:,:) # else real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:) real(r8), intent(inout) :: ad_sustr(LBi:,LBj:) real(r8), intent(inout) :: ad_svstr(LBi:,LBj:) real(r8), intent(out) :: ad_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,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) real(r8), intent(inout) :: ad_stflx(LBi:UBi,LBj:UBj,NT(ng)) # else real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_sustr(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: ad_svstr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, is, itrc, j, k # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif real(r8), parameter :: Aspv = 0.0_r8 real(r8) :: cff, cff1, scale # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize adjoint state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend ad_state(is)=Aspv END DO # endif ! !----------------------------------------------------------------------- ! Load adjoint STATE variables into full 1D state vector. !----------------------------------------------------------------------- ! ! Initialize local summations. ! IF (IntTrap.eq.1) THEN DO is=Mstr,Mend STORAGE(ng)%ad_Work(is)=0.0_r8 END DO END IF ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! cff1=dt(ng)*REAL(ntimes(ng)/Nintervals,r8) IF ((IntTrap.eq.1).or.(IntTrap.eq.Nintervals+1)) THEN cff1=0.5_r8*cff1 ENDIF ! ! Load adjoint of free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif scale=cff1 # ifdef ENERGYNORM_SCALE scale=scale/SQRT(0.5_r8*g*rho0) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_zeta(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_zeta(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Load adjoint of 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(h(i-1,j)+h(i,j))) # else scale=cff1 # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_ubar(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_ubar(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF ! ! Load adjoint of 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(h(i,j-1)+h(i,j))) # else scale=cff1 # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_vbar(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_vbar(i,j,kstp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # else ! ! Load adjoint of 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # else scale=cff1 # endif is=IJwaterU(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_u(i,j,k,nstp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # else scale=cff1 # endif is=(i-Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_u(i,j,k,nstp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF ! ! Load adjoint of 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # else scale=cff1 # endif is=IJwaterV(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_v(i,j,k,nstp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # else scale=cff1 # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_v(i,j,k,nstp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF ! ! Load adjoint of tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*Hz(i,j,k)) # else scale=cff1 # endif is=IJwaterR(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_t(i,j,k,nstp,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # ifdef ENERGYNORM_SCALE scale=cff1/SQRT(cff*Hz(i,j,k)) # else scale=cff1 # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_t(i,j,k,nstp,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF END DO # endif ! ! Load adjoint of surface U-stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=cff1 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_sustr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_sustr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF ! ! Load adjoint of surface V-stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=cff1 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_svstr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_svstr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Load adjoint of surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=cff1 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_stflx(i,j,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scale*ad_stflx(i,j,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE ad_so_pack_tile # else ! SUBROUTINE ad_so_pack_red (ng, tile, Mstr, Mend, IntTrap, & & ad_state) ! !======================================================================= ! ! ! This routine packs the adjoint variables into the state vector. ! ! The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_forces USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend integer, intent(in) :: IntTrap # ifdef ASSUMED_SHAPE real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_so_pack_red" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_so_pack_red_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif & IntTrap, & # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, ad_state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % ad_t, & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & & FORCES(ng) % ad_stflx, & # else & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & # endif & OCEAN(ng) % ad_zeta, & & FORCES(ng) % ad_sustr, & & FORCES(ng) % ad_svstr) # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Scatter (global to threaded) adjoint state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iADM, Mstr, Mend, Mstate(ng), & & Swork, ad_state) # endif ! RETURN END SUBROUTINE ad_so_pack_red ! !*********************************************************************** SUBROUTINE ad_so_pack_red_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & IntTrap, & & Mstr, Mend, ad_state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & ad_t, ad_u, ad_v, ad_stflx, & # else & ad_ubar, ad_vbar, & # endif & ad_zeta, ad_sustr, ad_svstr) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_grid USE mod_iounits USE mod_ncparam USE mod_netcdf !! USE mod_ocean # if defined PIO_LIB && defined DISTRIBUTE USE mod_pio_netcdf # endif USE mod_scalars USE mod_storage ! USE nf_fread2d_mod, ONLY : nf_fread2d # ifdef SOLVE3D USE nf_fread3d_mod, ONLY : nf_fread3d # endif USE strings_mod, ONLY : FoundError ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp integer, intent(in) :: IntTrap # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_stflx(LBi:,LBj:,:) # else real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:) real(r8), intent(inout) :: ad_sustr(LBi:,LBj:) real(r8), intent(inout) :: ad_svstr(LBi:,LBj:) real(r8), intent(out) :: ad_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,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) real(r8), intent(inout) :: ad_stflx(LBi:UBi,LBj:UBj,NT(ng)) # else real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_sustr(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: ad_svstr(LBi:UBi,LBj:UBj) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, ifield, is, itrc, j, k integer :: Fcount, Irec, Nrec integer :: gtype, status integer :: Vsize(4), ntAD, ntTL, Iinp # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif ! real(r8), parameter :: Aspv = 0.0_r8 real(dp) :: scale real(r8) :: Fmin, Fmax real(r8) :: cff, cff1, afac, scalev ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_so_pack_red_tile" # if defined PIO_LIB && defined DISTRIBUTE ! TYPE (IO_Desc_t), pointer :: ioDesc # endif # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize adjoint state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend ad_state(is)=Aspv END DO # endif ! !----------------------------------------------------------------------- ! Load adjoint STATE variables into full 1D state vector. !----------------------------------------------------------------------- ! ! Initialize local summations. ! IF (IntTrap.eq.1) THEN DO is=Mstr,Mend STORAGE(ng)%ad_Work(is)=0.0_r8 END DO END IF ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+(Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Read in adjoint state variables records. ! cff1=dt(ng)*REAL(ntimes(ng)/Nintervals,r8) cff1=cff1*cff1 ! DO i=1,4 Vsize(i)=0 END DO ! ! Loop over the number of number of ADJ netcdf file records. ! Iinp=1 scale=1.0_dp ntTL=(IntTrap-1)*nADJ(ng)+1 Fcount=ADM(ng)%Fcount Nrec=ADM(ng)%Nrec(Fcount) ! ADREC_LOOP : DO Irec=1,Nrec ! ! Autocorrelation factor. ! ntAD=(Nrec-Irec)*nADJ(ng)+1 CALL sp_bcoef (ng, ntAD, ntTL, afac) !AMM afac=afac*cff1 ! ! Process free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=r2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idFsur), & & ADM(ng)%Vid(idFsur), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_zeta(:,:,Iinp)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idFsur)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_zeta).eq.8) THEN ioDesc => ioDesc_dp_r2dvar(ng) ELSE ioDesc => ioDesc_sp_r2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idFsur), & & ADM(ng)%pioVar(idFsur), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_zeta(:,:,Iinp)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idFsur)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of free-surface. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif scalev=afac # if defined ENERGYNORM_SCALE scalev=scalev/SQRT(0.5_r8*g*rho0) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_zeta(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_zeta(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Process 2D U-momentum. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=u2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idUbar), & & ADM(ng)%Vid(idUbar), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_ubar(:,:,Iinp)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUbar)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_ubar).eq.8) THEN ioDesc => ioDesc_dp_u2dvar(ng) ELSE ioDesc => ioDesc_sp_u2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idUbar), & & ADM(ng)%pioVar(idUbar), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_ubar(:,:,Iinp)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUbar)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of 2D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO j=JR_RANGE DO i=IU_RANGE # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(h(i-1,j)+h(i,j))) # else scalev=afac # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_ubar(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_ubar(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF ! ! Process 2D V-momentum. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=v2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idVbar), & & ADM(ng)%Vid(idVbar), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_vbar(:,:,Iinp)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVbar)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_vbar).eq.8) THEN ioDesc => ioDesc_dp_v2dvar(ng) ELSE ioDesc => ioDesc_sp_v2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idVbar), & & ADM(ng)%pioVar(idVbar), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_vbar(:,:,Iinp)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVbar)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of 2D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO j=JV_RANGE DO i=IR_RANGE # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(h(i,j-1)+h(i,j))) # else scalev=afac # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_vbar(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_vbar(i,j,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # else ! ! Process 3D U-momentum. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=u3dvar status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idUvel), & & ADM(ng)%Vid(idUvel), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_u(:,:,:,Iinp)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUvel)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_u).eq.8) THEN ioDesc => ioDesc_dp_u3dvar(ng) ELSE ioDesc => ioDesc_sp_u3dvar(ng) END IF status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idUvel), & & ADM(ng)%pioVar(idUvel), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_u(:,:,:,Iinp)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUvel)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of 3D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # else scalev=afac # endif is=IJwaterU(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_u(i,j,k,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # else scalev=afac # endif is=(i-Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_u(i,j,k,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF ! ! Process 3D V-momentum. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=v3dvar status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idVvel), & & ADM(ng)%Vid(idVvel), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_v(:,:,:,Iinp)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVvel)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_v).eq.8) THEN ioDesc => ioDesc_dp_v3dvar(ng) ELSE ioDesc => ioDesc_sp_v3dvar(ng) END IF status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idVvel), & & ADM(ng)%pioVar(idVvel), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_v(:,:,:,Iinp)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVvel)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of 3D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # else scalev=afac # endif is=IJwaterV(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_v(i,j,k,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # else scalev=afac # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_v(i,j,k,Iinp) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF ! ! Process tracer type variables. ! DO itrc=1,NT(ng) ifield=isTvar(itrc) IF (SCALARS(ng)%Fstate(ifield)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=r3dvar status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,ifield), & & ADM(ng)%Tid(itrc), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_t(:,:,:,Iinp,itrc)) IF (FoundError(status, nf90_noerr, & & __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,ifield)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_t).eq.8) THEN ioDesc => ioDesc_dp_r3dvar(ng) ELSE ioDesc => ioDesc_sp_r3dvar(ng) END IF status=nf_fread3d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,ifield), & & ADM(ng)%pioTrc(itrc), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, 1, N(ng), & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_t(:,:,:,Iinp,itrc)) # endif IF (FoundError(status, PIO_noerr, & & __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,ifield)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of tracers variables. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # if defined ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*Hz(i,j,k)) # else scalev=afac # endif is=IJwaterR(i,j)+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_t(i,j,k,Iinp,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else # if defined ENERGYNORM_SCALE scalev=afac/SQRT(cff*Hz(i,j,k)) # else scalev=afac # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & scalev*ad_t(i,j,k,Iinp,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END DO END IF END DO # endif ! ! Process 2D U-momentum stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=u2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idUsms), & & ADM(ng)%Vid(idUsms), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_sustr(:,:)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUsms)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_sustr).eq.8) THEN ioDesc => ioDesc_dp_u2dvar(ng) ELSE ioDesc => ioDesc_sp_u2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idUsms), & & ADM(ng)%pioVar(idUsms), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % umask, & # endif & ad_sustr(:,:)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idUsms)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of surface U-stress. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_sustr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_sustr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF ! ! Process 2D V-momentum stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=v2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,idVsms), & & ADM(ng)%Vid(idVsms), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_svstr(:,:)) IF (FoundError(status, nf90_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVsms)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_svstr).eq.8) THEN ioDesc => ioDesc_dp_v2dvar(ng) ELSE ioDesc => ioDesc_sp_v2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,idVsms), & & ADM(ng)%pioVar(idVsms), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % vmask, & # endif & ad_svstr(:,:)) # endif IF (FoundError(status, PIO_noerr, __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,idVsms)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load adjoint of surface V-stress. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_svstr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_svstr(i,j) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Process surface tracer flux. ! DO itrc=1,NT(ng) ifield=idTsur(itrc) IF (SCALARS(ng)%Fstate(ifield)) THEN SELECT CASE (ADM(ng)%IOtype) CASE (io_nf90) gtype=r2dvar status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%ncid, Vname(1,ifield), & & ADM(ng)%Vid(ifield), Irec, & & gtype, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_stflx(:,:,itrc)) IF (FoundError(status, nf90_noerr, & & __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,ifield)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF # if defined PIO_LIB && defined DISTRIBUTE CASE (io_pio) IF (KIND(ad_stflx).eq.8) THEN ioDesc => ioDesc_dp_r2dvar(ng) ELSE ioDesc => ioDesc_sp_r2dvar(ng) END IF status=nf_fread2d(ng, iADM, ADM(ng)%name, & & ADM(ng)%pioFile, Vname(1,ifield), & & ADM(ng)%pioVar(ifield), Irec, & & ioDesc, Vsize, & & LBi, UBi, LBj, UBj, & & scale, Fmin, Fmax, & # ifdef MASKING & GRID(ng) % rmask, & # endif & ad_stflx(:,:,itrc)) # endif IF (FoundError(status, PIO_noerr, & & __LINE__, MyFile)) THEN IF (Master) WRITE (stdout,10) TRIM(Vname(1,ifield)), & & Irec, TRIM(ADM(ng)%name) exit_flag=2 ioerror=status RETURN END IF END SELECT ! ! Load surface tracer flux. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_stflx(i,j,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) ad_state(is)=STORAGE(ng)%ad_Work(is)+ & & afac*ad_stflx(i,j,itrc) STORAGE(ng)%ad_Work(is)=ad_state(is) # endif END DO END DO END IF END DO # endif END DO ADREC_LOOP ! 10 FORMAT (/,' AD_SO_PACK_RED - error while reading variable: ',a,2x,& & 'at time record = ',i3,/,17x,'in input NetCDF file: ',a) ! RETURN END SUBROUTINE ad_so_pack_red_tile # endif # elif defined ADJOINT ! SUBROUTINE ad_pack (ng, tile, Mstr, Mend, ad_state) ! !======================================================================= ! ! ! This routine packs the adjoint variables into the state vector. ! ! The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_pack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, ad_state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % ad_t, & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & # else & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & # endif & OCEAN(ng) % ad_zeta) # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Scatter (global to threaded) adjoint state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iADM, Mstr, Mend, Mstate(ng), & & Swork, ad_state) # endif ! RETURN END SUBROUTINE ad_pack ! !*********************************************************************** SUBROUTINE ad_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, ad_state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & ad_t, ad_u, ad_v, & # else & ad_ubar, ad_vbar, & # endif & ad_zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:) # else real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:) real(r8), intent(out) :: ad_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,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) # else real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, is, itrc, j, k integer, dimension(5+NT(ng)) :: offset real(r8), parameter :: Aspv = 0.0_r8 real(r8) :: cff, scale # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize adjoint state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend ad_state(is)=Aspv END DO # endif ! !----------------------------------------------------------------------- ! Load adjoint STATE variables into full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! # ifdef SOLVE3D # ifdef MASKING offset(isFsur)=0 offset(isUvel)=offset(isFsur)+NwaterR(ng) offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END DO # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END DO # else offset(isFsur)=0 offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END DO # endif # endif # else # ifdef MASKING offset(isFsur)=0 offset(isUbar)=offset(isFsur)+NwaterR(ng) offset(isVbar)=offset(isUbar)+NwaterU(ng) # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) # else offset(isFsur)=0 offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) # endif # endif # endif ! ! Load adjoint free-surface. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_state(is)=scale*ad_zeta(i,j,kstp) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_state(is)=scale*ad_zeta(i,j,kstp) # endif END DO END DO # ifndef SOLVE3D ! ! Load adjoint 2D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_state(is)=scale*ad_ubar(i,j,kstp) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_state(is)=scale*ad_ubar(i,j,kstp) # endif END DO END DO ! ! Load adjoint 2D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_state(is)=scale*ad_vbar(i,j,kstp) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_state(is)=scale*ad_vbar(i,j,kstp) # endif END DO END DO # else ! ! Load adjoint 3D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd ad_state(is)=scale*ad_u(i,j,k,nstp) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*ad_u(i,j,k,nstp) # endif END DO END DO END DO ! ! Load adjoint 3D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd ad_state(is)=scale*ad_v(i,j,k,nstp) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*ad_v(i,j,k,nstp) # endif END DO END DO END DO ! ! Load adjoint tracers variables. For now, use salinity scale for ! passive tracers. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif DO itrc=1,NT(ng) # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd ad_state(is)=scale*ad_t(i,j,k,nstp,itrc) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_state(is)=scale*ad_t(i,j,k,nstp,itrc) # endif END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE ad_pack_tile # endif # if defined ADJOINT && (defined SO_SEMI || defined STOCHASTIC_OPT) ! SUBROUTINE ad_unpack (ng, tile, Mstr, Mend, state) ! !======================================================================= ! ! ! This routine unpacks the requested adjoint state and/or surface ! ! forcing variables used in stochastic optimals. The state vector ! ! contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_forces USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_gather_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_unpak" ! # include "tile.h" ! # ifdef DISTRIBUTE ! ! Gather (threaded to global) adjoint state solution from all ! distributed nodes. ! CALL mp_gather_state (ng, iNLM, Mstr, Mend, Mstate(ng), & & state, Swork) ! # endif # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif CALL ad_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & # ifdef STOCHASTIC_OPT & knew(ng), & # else & kstp(ng), & # endif # ifdef SOLVE3D & nstp(ng), & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif # ifdef ENERGYNORM_SCALE & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & # endif # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork) # else & Mstr, Mend, state) # endif # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE ad_unpack ! !*********************************************************************** SUBROUTINE ad_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kout, & # ifdef SOLVE3D & nout, & # endif # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif # ifdef ENERGYNORM_SCALE & h, & # ifdef SOLVE3D & Hz, & # endif # endif & Mstr, Mend, state) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_ocean USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: kout # ifdef SOLVE3D integer, intent(in) :: nout # endif integer, intent(in) :: Mstr, Mend ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif # ifdef ENERGYNORM_SCALE real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) # endif # endif real(r8), intent(in) :: state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif # ifdef ENERGYNORM_SCALE real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) # endif # endif real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, is, itrc, j, k # ifdef SOLVE3D # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif # else integer, dimension(5) :: offset # endif real(r8) :: cff, scale # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Extract adjoint FORCING variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Extract adjoint free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) OCEAN(ng)%ad_zeta(i,j,kout)=scale*state(is) ELSE OCEAN(ng)%ad_zeta(i,j,kout)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isFsur) OCEAN(ng)%ad_zeta(i,j,kout)=scale*state(is) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Extract adjoint 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) OCEAN(ng)%ad_ubar(i,j,kout)=scale*state(is) ELSE OCEAN(ng)%ubar(i,j,kout)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) OCEAN(ng)%ad_ubar(i,j,kout)=scale*state(is) # endif END DO END DO END IF ! ! Extract adjoint 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) OCEAN(ng)%ad_vbar(i,j,kout)=scale*state(is) ELSE OCEAN(ng)%ad_vbar(i,j,kout)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) OCEAN(ng)%ad_vbar(i,j,kout)=scale*state(is) # endif END DO END DO END IF # else ! ! Extract adjoint 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd OCEAN(ng)%ad_u(i,j,k,nout)=scale*state(is) ELSE OCEAN(ng)%ad_u(i,j,k,nout)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd OCEAN(ng)%ad_u(i,j,k,nout)=scale*state(is) # endif END DO END DO END DO END IF ! ! Extract adjoint 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd OCEAN(ng)%ad_v(i,j,k,nout)=scale*state(is) ELSE OCEAN(ng)%ad_v(i,j,k,nout)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd OCEAN(ng)%ad_v(i,j,k,nout)=scale*state(is) # endif END DO END DO END DO END IF ! ! Extract adjoint tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd OCEAN(ng)%ad_t(i,j,k,nout,itrc)=scale*state(is) ELSE OCEAN(ng)%ad_t(i,j,k,nout,itrc)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd OCEAN(ng)%ad_t(i,j,k,nout,itrc)=scale*state(is) # endif END DO END DO END DO END IF END DO # endif ! ! Extract adjoint surface U-momentum stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) FORCES(ng)%ad_sustr(i,j)=scale*state(is) ELSE FORCES(ng)%ad_sustr(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) FORCES(ng)%ad_sustr(i,j)=scale*state(is) # endif END DO END DO END IF ! ! Extract adjoint surface V-momentum stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=1.0_r8 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) FORCES(ng)%ad_svstr(i,j)=scale*state(is) ELSE FORCES(ng)%ad_svstr(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) FORCES(ng)%ad_svstr(i,j)=scale*state(is) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Extract adjoint surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTvar(itrc)) FORCES(ng)%ad_stflx(i,j,itrc)=scale*state(is) ELSE FORCES(ng)%ad_stflx(i,j,itrc)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTvar(itrc)) FORCES(ng)%ad_stflx(i,j,itrc)=scale*state(is) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE ad_unpack_tile # elif defined ADJOINT ! SUBROUTINE ad_unpack (ng, tile, Mstr, Mend, state) ! !======================================================================= ! ! ! This routine unpacks the adjoint model variables from the state ! ! vector. If applicable, the state vector includes only unmasked ! ! water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_gather_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", ad_unpack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iADM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Gather (threaded to global) adjoint state solution from all ! distributed nodes. ! CALL mp_gather_state (ng, iTLM, Mstr, Mend, Mstate(ng), & & state, Swork) ! # endif CALL ad_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & knew(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % ad_t, & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & # else & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & # endif & OCEAN(ng) % ad_zeta) # ifdef PROFILE CALL wclock_off (ng, iADM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE ad_unpack ! !*********************************************************************** SUBROUTINE ad_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & knew, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & ad_t, ad_u, ad_v, & # else & ad_ubar, ad_vbar, & # endif & ad_zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: knew # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: state(Mstr:) real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:) # else real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,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) # else real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, is, itrc, j, k integer, dimension(5+NT(ng)) :: offset real(r8) :: cff, scale # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Extract adjoint state variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! # ifdef SOLVE3D # ifdef MASKING offset(isFsur)=0 offset(isUvel)=offset(isFsur)+NwaterR(ng) offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END DO # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END DO # else offset(isFsur)=0 offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END DO # endif # endif # else # ifdef MASKING offset(isFsur)=0 offset(isUbar)=offset(isFsur)+NwaterR(ng) offset(isVbar)=offset(isUbar)+NwaterU(ng) # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) # else offset(isFsur)=0 offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) # endif # endif # endif ! ! Extract adjoint free-surface. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) ad_zeta(i,j,knew)=scale*state(is) ELSE ad_zeta(i,j,knew)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) ad_zeta(i,j,knew)=scale*state(is) # endif END DO END DO # ifndef SOLVE3D ! ! Extract adjoint 2D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # if define ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) ad_ubar(i,j,knew)=scale*state(is) ELSE ad_ubar(i,j,knew)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) ad_ubar(i,j,knew)=scale*state(is) # endif END DO END DO ! ! Extract adjoint 2D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) ad_vbar(i,j,knew)=scale*state(is) ELSE ad_vbar(i,j,knew)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) ad_vbar(i,j,knew)=scale*state(is) # endif END DO END DO # else ! ! Extract adjoint 3D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd ad_u(i,j,k,nstp)=scale*state(is) ELSE ad_u(i,j,k,nstp)=0.0_r8 END IF # else # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd ad_u(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO ! ! Extract adjoint 3D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # if defined ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd ad_v(i,j,k,nstp)=scale*state(is) ELSE ad_v(i,j,k,nstp)=0.0_r8 END IF # else # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_v(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO ! ! Extract tangent linear tracers variables. For now, use salinity scale ! for passive tracers. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif DO itrc=1,NT(ng) # if defined ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd ad_t(i,j,k,nstp,itrc)=scale*state(is) ELSE ad_t(i,j,k,nstp,itrc)=0.0_r8 END IF # else # if defined ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd ad_t(i,j,k,nstp,itrc)=scale*state(is) # endif END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE ad_unpack_tile # endif # ifdef TANGENT ! SUBROUTINE tl_pack (ng, tile, Mstr, Mend, tl_state) ! !======================================================================= ! ! ! This routine packs the tangent linear variables into the state ! ! vetor. The state vector contains only interior water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_scatter_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(out) :: tl_state(Mstr:) # else real(r8), intent(out) :: tl_state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", tl_pack" # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iTLM, 2, __LINE__, MyFile) # endif CALL tl_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & krhs(ng), kstp(ng), knew(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, tl_state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & # else & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & # endif & OCEAN(ng) % tl_zeta) # ifdef DISTRIBUTE ! ! Scatter (global to threaded) tangent linear state solution to all ! distributed nodes. ! CALL mp_scatter_state (ng, iTLM, Mstr, Mend, Mstate(ng), & & Swork, tl_state) # endif # ifdef PROFILE CALL wclock_off (ng, iTLM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE tl_pack ! !*********************************************************************** SUBROUTINE tl_pack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & krhs, kstp, knew, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, tl_state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & tl_t, tl_u, tl_v, & # else & tl_ubar, tl_vbar, & # endif & tl_zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: krhs, kstp, knew # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(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:,:,:) # else real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) real(r8), intent(out) :: tl_state(Mstr:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) 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) # else real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: tl_state(Mstr:Mend) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, is, itrc, j, k integer, dimension(5+NT(ng)) :: offset real(r8), parameter :: Aspv = 0.0_r8 real(r8) :: cff, scale # include "set_bounds.h" # ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Initialize tangent linear state vector with special value (zero) to ! facilitate gathering/scattering communications between all nodes. ! This is achieved by summing all the buffers. !----------------------------------------------------------------------- ! DO is=Mstr,Mend tl_state(is)=Aspv END DO # endif ! !----------------------------------------------------------------------- ! Load tangent linear state variables into full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! # ifdef SOLVE3D # ifdef MASKING offset(isFsur)=0 offset(isUvel)=offset(isFsur)+NwaterR(ng) offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END DO # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END DO # else offset(isFsur)=0 offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END DO # endif # endif # else # ifdef MASKING offset(isFsur)=0 offset(isUbar)=offset(isFsur)+NwaterR(ng) offset(isVbar)=offset(isUbar)+NwaterU(ng) # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) # else offset(isFsur)=0 offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) # endif # endif # endif ! ! Load tangent linear free-surface. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) tl_state(is)=scale*tl_zeta(i,j,knew) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) tl_state(is)=scale*tl_zeta(i,j,knew) # endif END DO END DO # ifndef SOLVE3D ! ! Load tangent linear 2D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) tl_state(is)=scale*tl_ubar(i,j,knew) END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) tl_state(is)=scale*tl_ubar(i,j,knew) # endif END DO END DO ! ! Load tangent linear 2D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) tl_state(is)=scale*tl_vbar(i,j,knew) END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) tl_state(is)=scale*tl_vbar(i,j,knew) # endif END DO END DO # else ! ! Load tangent linear 3D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd tl_state(is)=scale*tl_u(i,j,k,nstp) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd tl_state(is)=scale*tl_u(i,j,k,nstp) # endif END DO END DO END DO ! ! Load tangent linear 3D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd tl_state(is)=scale*tl_v(i,j,k,nstp) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_state(is)=scale*tl_v(i,j,k,nstp) # endif END DO END DO END DO ! ! Load tangent linear tracers variables. For now, use salinity scale ! for passive tracers. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif DO itrc=1,NT(ng) # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd tl_state(is)=scale*tl_t(i,j,k,nstp,itrc) END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_state(is)=scale*tl_t(i,j,k,nstp,itrc) # endif END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE tl_pack_tile # endif # if defined TANGENT && (defined STOCHASTIC_OPT || \ defined HESSIAN_SV ) ! SUBROUTINE tl_unpack (ng, tile, Mstr, Mend, state) ! !======================================================================= ! ! ! This routine unpacks the tangent linear variables from the state ! ! vector. If applicable, the state vector includes only unmasked ! ! water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_forces USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_gather_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", tl_unpack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iTLM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Gather (threaded to global) tangent linear state solution from all ! distributed nodes. ! CALL mp_gather_state (ng, iTLM, Mstr, Mend, Mstate(ng), & & state, Swork) ! # endif CALL tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & # else & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & # endif & OCEAN(ng) % tl_zeta, & # ifdef SOLVE3D & FORCES(ng) % tl_stflx, & # endif & FORCES(ng) % tl_sustr, & & FORCES(ng) % tl_svstr) # ifdef PROFILE CALL wclock_off (ng, iTLM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE tl_unpack ! !*********************************************************************** SUBROUTINE tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & tl_t, tl_u, tl_v, & # else & tl_ubar, tl_vbar, & # endif & tl_zeta, & # ifdef SOLVE3D & tl_stflx, & # endif & tl_sustr, tl_svstr) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_ocean USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: state(Mstr:) real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(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_stflx(LBi:,LBj:,:) # else real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) real(r8), intent(inout) :: tl_sustr(LBi:,LBj:) real(r8), intent(inout) :: tl_svstr(LBi:,LBj:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) 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_stflx(LBi:UBi,LBj:UBj,NT(ng)) # else real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_sustr(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: tl_svstr(LBi:UBi,LBj:UBj) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, is, itrc, j, k # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif real(r8) :: cff, scale # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Extract tangent linear FORCING variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Extract tangent linear free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) tl_zeta(i,j,kstp)=scale*state(is) ELSE tl_zeta(i,j,kstp)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isFsur) tl_zeta(i,j,kstp)=scale*state(is) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Extract tangent linear 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) tl_ubar(i,j,kstp)=scale*state(is) ELSE tl_ubar(i,j,kstp)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) tl_ubar(i,j,kstp)=scale*state(is) # endif END DO END DO END IF ! ! Extract tangent linear 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) tl_vbar(i,j,kstp)=scale*state(is) ELSE tl_vbar(i,j,kstp)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) tl_vbar(i,j,kstp)=scale*state(is) # endif END DO END DO END IF # else ! ! Extract tangent linear 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd tl_u(i,j,k,nstp)=scale*state(is) ELSE tl_u(i,j,k,nstp)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd tl_u(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO END IF ! ! Extract tangent linear 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd tl_v(i,j,k,nstp)=scale*state(is) ELSE tl_v(i,j,k,nstp)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_v(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO END IF ! ! Extract tangent linear tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd tl_t(i,j,k,nstp,itrc)=scale*state(is) ELSE tl_t(i,j,k,nstp,itrc)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_t(i,j,k,nstp,itrc)=scale*state(is) # endif END DO END DO END DO END IF END DO # endif ! ! Extract tangent linear surface U-momentum stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) tl_sustr(i,j)=scale*state(is) ELSE tl_sustr(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) tl_sustr(i,j)=scale*state(is) # endif END DO END DO END IF ! ! Extract tangent linear surface V-momentum stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=1.0_r8 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) tl_svstr(i,j)=scale*state(is) ELSE tl_svstr(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) tl_svstr(i,j)=scale*state(is) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Extract tangent linear surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) tl_stflx(i,j,itrc)=scale*state(is) ELSE tl_stflx(i,j,itrc)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) tl_stflx(i,j,itrc)=scale*state(is) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE tl_unpack_tile # elif defined TANGENT && defined FORCING_SV ! SUBROUTINE tl_unpack (ng, tile, Mstr, Mend, state) ! !======================================================================= ! ! ! This routine unpacks the tangent linear variables from the state ! ! vector. If applicable, the state vector includes only unmasked ! ! water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_forces USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_gather_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", tl_unpack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iTLM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Gather (threaded to global) tangent linear state solution from all ! distributed nodes. ! CALL mp_gather_state (ng, iTLM, Mstr, Mend, Mstate(ng), & & state, Swork) ! # endif CALL tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % f_t, & & OCEAN(ng) % f_u, & & OCEAN(ng) % f_v, & # endif & OCEAN(ng) % f_ubar, & & OCEAN(ng) % f_vbar, & & OCEAN(ng) % f_zeta, & # ifdef SOLVE3D & FORCES(ng) % tl_stflx, & # endif & FORCES(ng) % tl_sustr, & & FORCES(ng) % tl_svstr) # ifdef PROFILE CALL wclock_off (ng, iTLM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE tl_unpack ! !*********************************************************************** SUBROUTINE tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & f_t, f_u, f_v, & # endif & f_ubar, f_vbar, & & f_zeta, & # ifdef SOLVE3D & tl_stflx, & # endif & tl_sustr, tl_svstr) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_forces USE mod_ncparam USE mod_ocean USE mod_scalars ! # ifdef FORCING_SV USE exchange_2d_mod # ifdef SOLVE3D USE exchange_3d_mod # endif # ifdef DISTRIBUTE USE mp_exchange_mod, ONLY : mp_exchange2d # ifdef SOLVE3D USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d # endif # endif # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: state(Mstr:) real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(inout) :: f_t(LBi:,LBj:,:,:) real(r8), intent(inout) :: f_u(LBi:,LBj:,:) real(r8), intent(inout) :: f_v(LBi:,LBj:,:) real(r8), intent(inout) :: tl_stflx(LBi:,LBj:,:) # endif real(r8), intent(inout) :: f_ubar(LBi:,LBj:) real(r8), intent(inout) :: f_vbar(LBi:,LBj:) real(r8), intent(inout) :: f_zeta(LBi:,LBj:) real(r8), intent(inout) :: tl_sustr(LBi:,LBj:) real(r8), intent(inout) :: tl_svstr(LBi:,LBj:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: f_t(LBi:UBi,LBj:UBj,N(ng),NT(ng)) real(r8), intent(inout) :: f_u(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: f_v(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: tl_stflx(LBi:UBi,LBj:UBj,NT(ng)) # endif real(r8), intent(inout) :: f_ubar(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: f_vbar(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: f_zeta(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: tl_sustr(LBi:UBi,LBj:UBj) real(r8), intent(inout) :: tl_svstr(LBi:UBi,LBj:UBj) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, icount, is, itrc, j, k # ifdef SALINITY integer, dimension(7+2*NT(ng)) :: offset # else integer, dimension(7+2*(NT(ng)+1)) :: offset # endif real(r8) :: cff, scale # ifdef SOLVE3D real(r8) :: cff1, cff2 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF real(r8), dimension(IminS:ImaxS,0:N(ng)) :: DC # endif # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Extract tangent linear FORCING variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! ! First clear the "offset" array. ! offset=0 ! # ifdef SOLVE3D # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) END IF iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+NwaterV(ng) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+NwaterR(ng) END IF END IF END DO # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) END IF iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+ & & (Lm(ng)+2)*(Mm(ng)+1) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+ & & (Lm(ng)+2)*(Mm(ng)+2) END IF END IF END DO # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUvel)) THEN offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVvel)) THEN offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) END IF iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END IF END DO IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN IF (itrc.eq.1) THEN offset(isTsur(itrc))=offset(isVstr)+Lm(ng)*(Mm(ng)-Voff) ELSE offset(isTsur(itrc))=offset(isTsur(itrc-1))+Lm(ng)*Mm(ng) END IF END IF END DO # endif # endif # else # ifdef MASKING IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+NwaterR(ng) END IF IF (SCALARS(ng)%Fstate(isVbar)) THEN offset(isVbar)=offset(isUbar)+NwaterU(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+NwaterU(ng) END IF # else # ifdef FULL_GRID IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isVstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)+1)*(Mm(ng)+2) END IF # else IF (SCALARS(ng)%Fstate(isFsur)) THEN offset(isFsur)=0 END IF IF (SCALARS(ng)%Fstate(isUbar)) THEN offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isVbar) THEN offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isUstr)=0 END IF IF (SCALARS(ng)%Fstate(isUstr)) THEN offset(isVstr)=offset(isUstr)+(Lm(ng)-Uoff)*Mm(ng) END IF # endif # endif # endif ! ! Extract tangent linear free-surface. ! IF (SCALARS(ng)%Fstate(isFsur)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) f_zeta(i,j)=scale*state(is) ELSE f_zeta(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isFsur) f_zeta(i,j)=scale*state(is) # endif END DO END DO END IF # ifndef SOLVE3D ! ! Extract tangent linear 2D U-velocity. ! IF (SCALARS(ng)%Fstate(isUbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) f_ubar(i,j)=scale*state(is) ELSE f_ubar(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) f_ubar(i,j)=scale*state(is) # endif END DO END DO END IF ! ! Extract tangent linear 2D V-velocity. ! IF (SCALARS(ng)%Fstate(isVbar)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) f_vbar(i,j)=scale*state(is) ELSE f_vbar(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) f_vbar(i,j)=scale*state(is) # endif END DO END DO END IF # else ! ! Extract tangent linear 3D U-velocity. ! IF (SCALARS(ng)%Fstate(isUvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+iadd f_u(i,j,k)=scale*state(is) ELSE f_u(i,j,k)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+iadd f_u(i,j,k)=scale*state(is) # endif END DO END DO END DO ! ! Compute the forcing forcing for tl_ubar based on f_u. ! DO j=JR_RANGE DO i=IU_RANGE DC(i,0)=0.0_r8 CF(i,0)=0.0_r8 END DO DO k=1,N(ng) DO i=IU_RANGE DC(i,k)=0.5_r8*(Hz(i,j,k)+Hz(i-1,j,k)) DC(i,0)=DC(i,0)+DC(i,k) CF(i,0)=CF(i,0)+DC(i,k)*f_u(i,j,k) END DO END DO DO i=IU_RANGE cff1=1.0_r8/DC(i,0) cff2=CF(i,0)*cff1 # ifdef MASKING cff2=cff2*umask(i,j) # endif f_ubar(i,j)=cff2 END DO END DO END IF ! ! Extract tangent linear 3D V-velocity. ! IF (SCALARS(ng)%Fstate(isVvel)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+iadd f_v(i,j,k)=scale*state(is) ELSE f_v(i,j,k)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd f_v(i,j,k)=scale*state(is) # endif END DO END DO END DO ! ! Compute the forcing forcing for tl_vbar based on f_v. ! DO j=JV_RANGE IF (j.ge.JstrM) THEN DO i=IR_RANGE DC(i,0)=0.0_r8 CF(i,0)=0.0_r8 END DO DO k=1,N(ng) DO i=IR_RANGE DC(i,k)=0.5_r8*(Hz(i,j,k)+Hz(i,j-1,k)) DC(i,0)=DC(i,0)+DC(i,k) CF(i,0)=CF(i,0)+DC(i,k)*f_v(i,j,k) END DO END DO DO i=IR_RANGE cff1=1.0_r8/DC(i,0) cff2=CF(i,0)*cff1 # ifdef MASKING cff2=cff2*vmask(i,j) # endif f_vbar(i,j)=cff2 END DO END IF END DO END IF ! ! Extract tangent linear tracers variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTvar(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+iadd f_t(i,j,k,itrc)=scale*state(is) ELSE f_t(i,j,k,itrc)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+iadd f_t(i,j,k,itrc)=scale*state(is) # endif END DO END DO END DO END IF END DO # endif # ifdef FORCING_SV ! ! Impose periodic boundary conditions as appropriate. ! IF (EWperiodic(ng).or.NSperiodic(ng)) THEN CALL exchange_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_zeta) # ifndef SOLVE3D CALL exchange_u2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_ubar) CALL exchange_v2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, f_vbar) # else CALL exchange_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), f_u) CALL exchange_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), f_v) DO itrc=1,NT(ng) CALL exchange_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & f_t(:,:,:,itrc)) END DO # endif END IF # ifdef DISTRIBUTE CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & f_zeta) # ifndef SOLVE3D CALL mp_exchange2d (ng, tile, iTLM, 2, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & f_ubar, f_vbar) # else CALL mp_exchange3d (ng, tile, iTLM, 2, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), f_u, f_v) CALL mp_exchange4d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), 1, NT(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), f_t) # endif # endif # endif ! ! Extract tangent linear surface U-momentum stress. ! IF (SCALARS(ng)%Fstate(isUstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUstr) tl_sustr(i,j)=scale*state(is) ELSE tl_sustr(i,j)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUstr) tl_sustr(i,j)=scale*state(is) # endif END DO END DO END IF ! ! Extract tangent linear surface V-momentum stress. ! IF (SCALARS(ng)%Fstate(isVstr)) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif scale=1.0_r8 DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVstr) tl_svstr(i,j)=scale*state(is) ELSE tl_svstr(i,j)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVstr) tl_svstr(i,j)=scale*state(is) # endif END DO END DO END IF # ifdef SOLVE3D ! ! Extract tangent linear surface tracer flux variables. ! DO itrc=1,NT(ng) IF (SCALARS(ng)%Fstate(isTsur(itrc))) THEN # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif scale=1.0_r8 DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isTsur(itrc)) tl_stflx(i,j,itrc)=scale*state(is) ELSE tl_stflx(i,j,itrc)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isTsur(itrc)) tl_stflx(i,j,itrc)=scale*state(is) # endif END DO END DO END IF END DO # endif ! RETURN END SUBROUTINE tl_unpack_tile # elif defined TANGENT ! SUBROUTINE tl_unpack (ng, tile, Mstr, Mend, state) ! !======================================================================= ! ! ! This routine unpacks the tangent linear variables from the state ! ! vector. If applicable, the state vector includes only unmasked ! ! water points. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ocean USE mod_stepping # ifdef DISTRIBUTE USE mod_storage # endif # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_gather_state # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) # endif ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__//", tl_unpack" ! # include "tile.h" ! # ifdef PROFILE CALL wclock_on (ng, iTLM, 2, __LINE__, MyFile) # endif # ifdef DISTRIBUTE ! ! Gather (threaded to global) tangent linear state solution from all ! distributed nodes. ! CALL mp_gather_state (ng, iTLM, Mstr, Mend, Mstate(ng), & & state, Swork) ! # endif CALL tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & # endif # ifdef DISTRIBUTE & 1, Mstate(ng), Swork, & # else & Mstr, Mend, state, & # endif # ifdef MASKING & GRID(ng) % IJwaterR, & & GRID(ng) % IJwaterU, & & GRID(ng) % IJwaterV, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & # endif & GRID(ng) % h, & # ifdef SOLVE3D & GRID(ng) % Hz, & & OCEAN(ng) % tl_t, & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & # else & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & # endif & OCEAN(ng) % tl_zeta) # ifdef PROFILE CALL wclock_off (ng, iTLM, 2, __LINE__, MyFile) # endif ! RETURN END SUBROUTINE tl_unpack ! !*********************************************************************** SUBROUTINE tl_unpack_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & # endif & Mstr, Mend, state, & # ifdef MASKING & IJwaterR, IJwaterU, IJwaterV, & & rmask, umask, vmask, & # endif & h, & # ifdef SOLVE3D & Hz, & & tl_t, tl_u, tl_v, & # else & tl_ubar, tl_vbar, & # endif & tl_zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_ncparam USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: Mstr, Mend integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:,LBj:) integer, intent(in) :: IJwaterU(LBi:,LBj:) integer, intent(in) :: IJwaterV(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) # endif real(r8), intent(in) :: state(Mstr:) real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SOLVE3D real(r8), intent(in) :: Hz(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:,:,:) # else real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:) # else # ifdef MASKING integer, intent(in) :: IJwaterR(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterU(LBi:UBi,LBj:UBj) integer, intent(in) :: IJwaterV(LBi:UBi,LBj:UBj) real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj) real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) 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) # else real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:) # endif real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:) # endif ! ! Local variable declarations. ! # ifndef MASKING integer :: Imax, Ioff, Jmax, Joff # endif integer :: Uoff, Voff integer :: i, iadd, is, itrc, j, k integer, dimension(5+NT(ng)) :: offset real(r8) :: cff, scale # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Extract tangent linear STATE variables from full 1D state vector. !----------------------------------------------------------------------- ! ! Set offsets for momentum variables due to periodic boundary ! conditions. Recall that in East-West periodic boundary conditions ! IstrU=1. Otherwise, IstrU=2. Similarly, in North-South periodic ! applications IstrV=1 or else IstrV=2. ! IF (EWperiodic(ng)) THEN Uoff=0 ELSE Uoff=1 END IF ! IF (NSperiodic(ng)) THEN Voff=0 ELSE Voff=1 END IF ! ! Determine the index offset for each variable in the state vector. # ifdef MASKING ! Notice that in Land/Sea masking application the state vector only ! contains water points to avoid large null space. # endif ! # ifdef SOLVE3D # ifdef MASKING offset(isFsur)=0 offset(isUvel)=offset(isFsur)+NwaterR(ng) offset(isVvel)=offset(isUvel)+NwaterU(ng)*N(ng) iadd=NwaterV(ng)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=NwaterR(ng)*N(ng) END DO # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUvel)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVvel)=offset(isUvel)+(Lm(ng)+1)*(Mm(ng)+2)*N(ng) iadd=(Lm(ng)+2)*(Mm(ng)+1)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=(Lm(ng)+2)*(Mm(ng)+2)*N(ng) END DO # else offset(isFsur)=0 offset(isUvel)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVvel)=offset(isUvel)+(Lm(ng)-Uoff)*Mm(ng)*N(ng) iadd=Lm(ng)*(Mm(ng)-Voff)*N(ng) DO itrc=1,NT(ng) offset(isTvar(itrc))=offset(isTvar(itrc)-1)+iadd iadd=Lm(ng)*Mm(ng)*N(ng) END DO # endif # endif # else # ifdef MASKING offset(isFsur)=0 offset(isUbar)=offset(isFsur)+NwaterR(ng) offset(isVbar)=offset(isUbar)+NwaterU(ng) # else # ifdef FULL_GRID offset(isFsur)=0 offset(isUbar)=offset(isFsur)+(Lm(ng)+2)*(Mm(ng)+2) offset(isVbar)=offset(isUbar)+(Lm(ng)+1)*(Mm(ng)+2) # else offset(isFsur)=0 offset(isUbar)=offset(isFsur)+Lm(ng)*Mm(ng) offset(isVbar)=offset(isUbar)+(Lm(ng)-Uoff)*Mm(ng) # endif # endif # endif ! ! Extract tangent linear free-surface. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Ioff=0 Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(0.5_r8*g*rho0) # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN is=IJwaterR(i,j)+offset(isFsur) tl_zeta(i,j,kstp)=scale*state(is) ELSE tl_zeta(i,j,kstp)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isFsur) tl_zeta(i,j,kstp)=scale*state(is) # endif END DO END DO # ifndef SOLVE3D ! ! Extract tangent linear 2D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i-1,j)+h(i,j))) # endif # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN is=IJwaterU(i,j)+offset(isUbar) tl_ubar(i,j,kstp)=scale*state(is) ELSE tl_ubar(i,j,kstp)=0.0_r8 END IF # else is=(i-Ioff)+(j-Joff)*Imax+offset(isUbar) tl_ubar(i,j,kstp)=scale*state(is) # endif END DO END DO ! ! Extract tangent linear 2D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Ioff=1 Joff=1 # else Imax=Lm(ng) Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(h(i,j-1)+h(i,j))) # endif # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN is=IJwaterV(i,j)+offset(isVbar) tl_vbar(i,j,kstp)=scale*state(is) ELSE tl_vbar(i,j,kstp)=0.0_r8 END IF # else is=(i+Ioff)+(j-Joff)*Imax+offset(isVbar) tl_vbar(i,j,kstp)=scale*state(is) # endif END DO END DO # else ! ! Extract tangent linear 3D U-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+1 Jmax=Mm(ng)+2 Ioff=0 Joff=0 # else Imax=Lm(ng)-Uoff Jmax=Mm(ng) Ioff=Uoff Joff=1 # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterU(ng)+offset(isUvel) # else iadd=(k-1)*Imax*Jmax+offset(isUvel) # endif DO j=JR_RANGE DO i=IU_RANGE # ifdef MASKING IF (umask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=IJwaterU(i,j)+iadd tl_u(i,j,k,nstp)=scale*state(is) ELSE tl_u(i,j,k,nstp)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) # endif is=(i-Ioff)+(j-Joff)*Imax+iadd tl_u(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO ! ! Extract tangent linear 3D V-velocity. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+1 Ioff=1 Joff=1 # else Imax=Lm(ng) Jmax=Mm(ng)-Voff Ioff=0 Joff=1+Voff # endif # endif # ifdef ENERGYNORM_SCALE cff=0.25_r8*rho0 # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterV(ng)+offset(isVvel) # else iadd=(k-1)*Imax*Jmax+offset(isVvel) # endif DO j=JV_RANGE DO i=IR_RANGE # ifdef MASKING IF (vmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=IJwaterV(i,j)+iadd tl_v(i,j,k,nstp)=scale*state(is) ELSE tl_v(i,j,k,nstp)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_v(i,j,k,nstp)=scale*state(is) # endif END DO END DO END DO ! ! Extract tangent linear tracers variables. For now, use salinity scale ! for passive tracers. ! # ifndef MASKING # ifdef FULL_GRID Imax=Lm(ng)+2 Jmax=Mm(ng)+2 Ioff=1 Joff=0 # else Imax=Lm(ng) Jmax=Mm(ng) Ioff=0 Joff=1 # endif # endif DO itrc=1,NT(ng) # ifdef ENERGYNORM_SCALE IF (itrc.eq.itemp) THEN cff=0.5_r8*rho0*Tcoef(ng)*Tcoef(ng)*g*g/bvf_bak ELSE IF (itrc.eq.isalt) THEN cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak ELSE cff=0.5_r8*rho0*Scoef(ng)*Scoef(ng)*g*g/bvf_bak END IF # else scale=1.0_r8 # endif DO k=1,N(ng) # ifdef MASKING iadd=(k-1)*NwaterR(ng)+offset(isTvar(itrc)) # else iadd=(k-1)*Imax*Jmax+offset(isTvar(itrc)) # endif DO j=JR_RANGE DO i=IR_RANGE # ifdef MASKING IF (rmask(i,j).gt.0.0_r8) THEN # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=IJwaterR(i,j)+iadd tl_t(i,j,k,nstp,itrc)=scale*state(is) ELSE tl_t(i,j,k,nstp,itrc)=0.0_r8 END IF # else # ifdef ENERGYNORM_SCALE scale=1.0_r8/SQRT(cff*Hz(i,j,k)) # endif is=(i+Ioff)+(j-Joff)*Imax+iadd tl_t(i,j,k,nstp,itrc)=scale*state(is) # endif END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE tl_unpack_tile # endif # ifdef SO_SEMI # ifdef SO_SEMI_WHITE ! SUBROUTINE so_semi_white (ng, tile, Mstr, Mend, state, ad_state) ! !======================================================================= ! ! ! This routine computes a new stochastic optimals perturbation vector ! ! (seminorm estimation) assuming white noise forcing using ARPACK. ! ! ! !======================================================================= ! USE mod_param USE mod_parallel USE mod_scalars USE mod_storage # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_reduce # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! integer :: NSUB, is, rec real(r8) :: SOnorm, my_SOnorm, my_TRnorm real(r8) :: SOnorm1, my_SOnorm1 real(r8) :: cff, cff1, cff2 # ifdef DISTRIBUTE real(r8), dimension(3) :: rbuffer character (len=3), dimension(3) :: op_handle # endif # include "tile.h" ! !----------------------------------------------------------------------- ! Compute seminorm, stochastic optimals adjoint perturbation vector. !----------------------------------------------------------------------- ! ! Initialize adjoint state vector. ! DO is=Mstr,Mend ad_state(is)=0.0_r8 END DO ! ! Sum over all adjoint surface forcing records available in "so_state'. ! IF (Master) THEN WRITE (stdout,'(/)') END IF my_TRnorm=0.0_r8 ! DO rec=1,Nsemi(ng) ! ! Compute normalization factor. ! cff=REAL((nADJ(ng)-1)*(2*nADJ(ng)-1),r8)/REAL(6*nADJ(ng),r8) cff1=1.0_r8+cff cff2=0.5_r8*REAL((nADJ(ng)-1))-cff ! my_SOnorm=0.0_r8 my_SOnorm1=0.0_r8 DO is=Mstr,Mend my_SOnorm=my_SOnorm+ & & STORAGE(ng)%so_state(is,rec)*state(is) END DO ! IF (rec.ne.Nsemi(ng)) THEN DO is=Mstr,Mend my_SOnorm1=my_SOnorm1+ & & STORAGE(ng)%so_state(is,rec+1)*state(is) my_TRnorm=my_TRnorm+ & & cff1*STORAGE(ng)%so_state(is,rec)* & & STORAGE(ng)%so_state(is,rec)+ & & 2.0_r8*cff2*STORAGE(ng)%so_state(is,rec )* & & STORAGE(ng)%so_state(is,rec+1)+ & & cff*STORAGE(ng)%so_state(is,rec+1)* & & STORAGE(ng)%so_state(is,rec+1) END DO ELSE DO is=Mstr,Mend my_TRnorm=my_TRnorm+ & & STORAGE(ng)%so_state(is,rec)* & & STORAGE(ng)%so_state(is,rec) END DO END IF ! ! Global reduction of normalization factor. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else IF (DOMAIN(ng)%SouthWest_Corner(tile).and. & & DOMAIN(ng)%NorthEast_Corner(tile)) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF # endif !$OMP CRITICAL (SO_NORM) IF (tile_count.eq.0) THEN SOnorm=0.0_r8 SOnorm1=0.0_r8 IF (rec.eq.1) THEN TRnorm(ng)=0.0_r8 END IF END IF SOnorm=SOnorm+my_SOnorm SOnorm1=SOnorm1+my_SOnorm1 tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE rbuffer(1)=SOnorm rbuffer(2)=SOnorm1 op_handle(1)='SUM' op_handle(2)='SUM' CALL mp_reduce (ng, iADM, 2, rbuffer, op_handle) SOnorm=rbuffer(1) SOnorm1=rbuffer(2) # endif END IF !$OMP END CRITICAL (SO_NORM) ! ! Report normalization factors. ! IF (Master) THEN WRITE (stdout,10) rec, SOnorm, SOnorm1 10 FORMAT (3x,'Rec = ',i2.2,2x,'SOnorm = ',1p,e15.8,0p, & & 2x,'SOnorm1 = ',1p,e15.8) END IF ! ! Compute new perturbation vector. ! IF (rec.ne.Nsemi(ng)) THEN DO is=Mstr,Mend ad_state(is)=ad_state(is)+ & & cff1*SOnorm *STORAGE(ng)%so_state(is,rec )+ & & cff2*SOnorm1*STORAGE(ng)%so_state(is,rec )+ & & cff2*SOnorm *STORAGE(ng)%so_state(is,rec+1)+ & & cff *SOnorm1*STORAGE(ng)%so_state(is,rec+1) END DO ELSE DO is=Mstr,Mend ad_state(is)=ad_state(is)+ & & SOnorm*STORAGE(ng)%so_state(is,rec) END DO END IF END DO ! ! Global reduction of normalization factor, TRnorm. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else IF (DOMAIN(ng)%SouthWest_Corner(tile).and. & & DOMAIN(ng)%NorthEast_Corner(tile)) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF # endif !$OMP CRITICAL (TR_NORM) IF (tile_count.eq.0) THEN TRnorm(ng)=0.0_r8 END IF TRnorm(ng)=TRnorm(ng)+my_TRnorm tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE op_handle(1)='SUM' CALL mp_reduce (ng, iADM, 1, TRnorm(ng), op_handle(1)) # endif END IF !$OMP END CRITICAL (TR_NORM) ! RETURN END SUBROUTINE so_semi_white # else ! SUBROUTINE so_semi_red (ng, tile, Mstr, Mend, state, ad_state) ! !======================================================================= ! ! ! This routine computes a new stochastic optimals perturbation vector ! ! (seminorm estimation) assuming red noise forcing using ARPACK. ! ! ! !======================================================================= ! USE mod_param USE mod_parallel USE mod_iounits USE mod_storage USE mod_scalars # ifdef DISTRIBUTE ! USE distribute_mod, ONLY : mp_reduce # endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: Mstr, Mend # ifdef ASSUMED_SHAPE real(r8), intent(in) :: state(Mstr:) real(r8), intent(out) :: ad_state(Mstr:) # else real(r8), intent(in) :: state(Mstr:Mend) real(r8), intent(out) :: ad_state(Mstr:Mend) # endif ! ! Local variable declarations. ! integer :: NSUB, is, ntAD, ntTL, rec, rec1 real(r8) :: SOnorm, my_TRnorm real(r8), dimension(Nsemi(ng)) :: Bcoef real(r8), dimension(Nsemi(ng)) :: SOdotprod real(r8), dimension(Nsemi(ng)) :: my_dotprod # ifdef DISTRIBUTE character (len=3), dimension(Nsemi(ng)) :: op_handle # endif # include "tile.h" ! !----------------------------------------------------------------------- ! Compute seminorm, stochastic optimals adjoint perturbation vector. !----------------------------------------------------------------------- ! ! Initialize adjoint state vector. ! DO is=Mstr,Mend ad_state(is)=0.0_r8 END DO ! ! First compute the dot-products. ! DO rec=1,Nsemi(ng) my_dotprod(rec)=0.0_r8 DO is=Mstr,Mend my_dotprod(rec)=my_dotprod(rec)+ & & STORAGE(ng)%so_state(is,rec)*state(is) END DO END DO ! ! Global reduction of dot products. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else IF (DOMAIN(ng)%SouthWest_Corner(tile).and. & & DOMAIN(ng)%NorthEast_Corner(tile)) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF # endif !$OMP CRITICAL (SO_DOT) IF (tile_count.eq.0) THEN DO rec=1,Nsemi(ng) SOdotprod(rec)=0.0_r8 END DO END IF DO rec=1,Nsemi(ng) SOdotprod(rec)=SOdotprod(rec)+my_dotprod(rec) END DO tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE DO rec=1,Nsemi(ng) op_handle(rec)='SUM' END DO CALL mp_reduce (ng, iADM, Nsemi(ng), SOdotprod, op_handle) # endif END IF !$OMP END CRITICAL (SO_DOT) ! ! Second, loop over time twice allowing for the decorrelation due to the ! red noise AR(1) process with assumed decorrelation time SOdecay. ! IF (Master) THEN WRITE (stdout,'(/)') END IF my_TRnorm=0.0_r8 ! DO rec=1,Nsemi(ng) ntAD=(rec-1)*nADJ(ng)+1 SOnorm=0.0_r8 DO rec1=1,Nsemi(ng) ntTL=(rec1-1)*nADJ(ng)+1 CALL sp_bcoef (ng, ntAD, ntTL, Bcoef(rec1)) SOnorm=SOnorm+Bcoef(rec1)*SOdotprod(rec1) DO is=Mstr,Mend my_TRnorm=my_TRnorm+ & & STORAGE(ng)%so_state(is,rec )*Bcoef(rec1)* & & STORAGE(ng)%so_state(is,rec1) END DO END DO ! ! Report normalization factors. ! IF (Master) THEN WRITE (stdout,10) rec, SOdotprod(rec), Bcoef(rec), SOnorm 10 FORMAT (1x,'Rec = ',i2.2,1x,'SOdotprod = ',1p,e13.6,0p, & & 1x,'Bcoef = ',1p,e13.6,0p,1x,'SOnorm = ',1p,e13.6) END IF ! ! Compute new perturbation vector. ! DO is=Mstr,Mend ad_state(is)=ad_state(is)+ & & SOnorm*STORAGE(ng)%so_state(is,rec) END DO END DO ! ! Global reduction of normalization factor, TRnorm. ! # ifdef DISTRIBUTE NSUB=1 ! distributed-memory # else IF (DOMAIN(ng)%SouthWest_Corner(tile).and. & & DOMAIN(ng)%NorthEast_Corner(tile)) THEN NSUB=1 ! non-tiled application ELSE NSUB=NtileX(ng)*NtileE(ng) ! tiled application END IF # endif !$OMP CRITICAL (TR_NORM) IF (tile_count.eq.0) THEN TRnorm(ng)=0.0_r8 END IF TRnorm(ng)=TRnorm(ng)+my_TRnorm tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 # ifdef DISTRIBUTE op_handle(1)='SUM' CALL mp_reduce (ng, iADM, 1, TRnorm(ng), op_handle(1)) # endif END IF !$OMP END CRITICAL (TR_NORM) ! RETURN END SUBROUTINE so_semi_red # endif # endif # if defined SO_SEMI || !defined STOCH_OPT_WHITE ! SUBROUTINE sp_bcoef (ng, ntAD, ntTL, Bcoef) ! !======================================================================= ! ! ! This routine is used to compute red noise stochastic processes ! ! time-lagged coefficient, Bcoef, used to evaluate discrete ! ! double-time integrals. Currently, a discrete-time Markov chain ! ! model is assumed with autoregressive order-one processes, AR(1). ! ! Notice that the routine SP_ACOEF is called to compute the inner ! ! integral. ! ! ! !======================================================================= ! USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, ntAD, ntTL real(r8), intent(out):: Bcoef ! ! Local variable declarations. ! integer :: i, it1, it2 real(r8) :: Acoef, Acoef1, Acoef2, df1, rov ! !----------------------------------------------------------------------- ! Compute red noise stochastic process time-lagged coefficient to ! evaluate discrete double time-integrals. Currently, only Markov ! processes, AR(1), are considered. !----------------------------------------------------------------------- ! ! Here, ntAD is the current model timestep and ntTL is the timestep ! associated with forcing. ! rov=1.0_r8/REAL(nADJ(ng),r8) ! IF ((ntAD.gt.1).and.(ntAD.lt.ntimes(ng)+1)) THEN it1=ntAD it2=ntAD-nADJ(ng) CALL sp_acoef (ng, it1, ntTL, Acoef) Bcoef=Acoef DO i=1,nADJ(ng)-1 CALL sp_acoef (ng, it1+i, ntTL, Acoef1) CALL sp_acoef (ng, it2+i, ntTL, Acoef2) df1=REAL(i,r8)*rov Bcoef=Bcoef+(1.0_r8-df1)*Acoef1+df1*Acoef2 END DO ELSE IF (ntAD.eq.1) THEN CALL sp_acoef (ng, 1, ntTL, Acoef) Bcoef=Acoef DO i=1,nADJ(ng)-1 CALL sp_acoef (ng, 1+i, ntTL, Acoef1) df1=REAL(i,r8)*rov Bcoef=Bcoef+(1.0_r8-df1)*Acoef1 END DO ELSE IF (ntAD.eq.ntimes(ng)+1) THEN CALL sp_acoef (ng, ntimes(ng)+1, ntTL, Acoef) Bcoef=Acoef DO i=1,nADJ(ng)-1 CALL sp_acoef (ng, ntimes(ng)+1-nADJ(ng)+i, ntTL, Acoef1) df1=REAL(i,r8)*rov Bcoef=Bcoef+df1*Acoef1 END DO END IF ! RETURN END SUBROUTINE sp_bcoef ! SUBROUTINE sp_acoef (ng, ntAD, ntTL, Acoef) ! !======================================================================= ! ! ! This routine is used to compute red noise stochastic processes ! ! time-lagged coefficient, Acoef, used to evaluate inner time ! ! integral. Currently, a discrete-time Markov chain model is ! ! assumed with autoregressive order-one processes, AR(1). Notice ! ! that function SP_AUTOC is used to set autocorrelation model. ! ! ! !======================================================================= ! USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, ntAD, ntTL real(r8), intent(out):: Acoef ! ! Local variable declarations. ! integer :: i, idf1, idf2, idf4 real(r8) :: df3, rov ! !----------------------------------------------------------------------- ! Compute red noise stochastic process time-lagged coefficients to ! evaluate discrete inner time-integral. Currently, only Markov ! processes, AR(1), are considered. !----------------------------------------------------------------------- ! ! Here, ntAD is the current timestep corresponding to time when ! solution is saved. ! rov=1.0_r8/REAL(nADJ(ng),r8) IF ((ntTL.gt.1).and.(ntTL.lt.ntimes(ng)+1)) THEN Acoef=0.0_r8 DO i=1,nADJ(ng)-1 idf1=IABS(ntAD-ntTL-i)+1 idf2=IABS(ntAD-(ntTL-nADJ(ng))-i)+1 df3=REAL(i,r8)*rov Acoef=Acoef+sp_autoc(ng,idf1)*(1.0_r8-df3)+ & & sp_autoc(ng,idf2)*df3 END DO idf4=IABS(ntAD-ntTL)+1 Acoef=Acoef+sp_autoc(ng,idf4) ELSE IF (ntTL.eq.1) THEN Acoef=0.0_r8 DO i=1,nADJ(ng)-1 idf1=IABS(ntAD-1-i)+1 df3=REAL(i,r8)*rov Acoef=Acoef+sp_autoc(ng,idf1)*(1.0_r8-df3) END DO idf4=IABS(ntAD-1)+1 Acoef=Acoef+sp_autoc(ng,idf4) ELSE IF (ntTL.eq.ntimes(ng)+1) THEN Acoef=0.0_r8 DO i=1,nADJ(ng)-1 idf2=IABS(ntAD-ntimes(ng)-1+nADJ(ng)-i)+1 df3=REAL(i,r8)*rov Acoef=Acoef+sp_autoc(ng,idf2)*df3 END DO idf4=IABS(ntAD-ntimes(ng)-1)+1 Acoef=Acoef+sp_autoc(ng,idf4) END IF ! RETURN END SUBROUTINE sp_acoef ! FUNCTION sp_autoc (ng, idf) ! !======================================================================= ! ! ! This routine is used to compute red noise stochastic processes ! ! autocorrelation model. Notice that only AR(1) processes are ! ! considered. However, other models can be easily implemented in ! ! terms of look tables. ! ! !======================================================================= ! USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, idf ! ! Function result. ! real(r8) :: sp_autoc ! !----------------------------------------------------------------------- ! Set autocorrelation model. !----------------------------------------------------------------------- # ifdef SO_NON_AR1 ! ! Use a user-defined temporal decorrelation function such as in the ! form of a look-up table computed from actual data. ! sp_autoc=0.0_r8 # else ! ! Assume an AR(1) process with decorrelation time SO_decay. ! sp_autoc=EXP(-ABS(REAL(idf-1,r8))*dt(ng)/SO_decay(ng)) # endif ! RETURN END FUNCTION sp_autoc # endif #endif #undef IR_RANGE #undef IU_RANGE #undef JR_RANGE #undef JV_RANGE END MODULE packing_mod