SUBROUTINE ad_main3d (RunInterval) ! !git $Id$ !svn $Id: ad_main3d.F 1166 2023-05-17 20:11:58Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group Andrew M. Moore ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! This subroutine is the main driver for adjoint ROMS/TOMS when ! ! configurated as a full 3D baroclinic ocean model. It advances ! ! backwards the adjoint model equations for all nested grids, if ! ! any, by the specified time interval (seconds), RunInterval. ! ! ! !======================================================================= ! USE mod_param USE mod_parallel USE mod_fourdvar USE mod_iounits USE mod_scalars USE mod_stepping ! USE ad_diag_mod, ONLY : ad_diag USE ad_ini_fields_mod, ONLY : ad_ini_fields, ad_ini_zeta USE ad_force_dual_mod, ONLY : ad_force_dual USE ad_htobs_mod, ONLY : ad_htobs USE ad_omega_mod, ONLY : ad_omega USE ad_rho_eos_mod, ONLY : ad_rho_eos USE ad_rhs3d_mod, ONLY : ad_rhs3d USE ad_set_depth_mod, ONLY : ad_set_depth USE ad_set_massflux_mod, ONLY : ad_set_massflux USE ad_set_vbc_mod, ONLY : ad_set_vbc USE ad_set_zeta_mod, ONLY : ad_set_zeta USE ad_step2d_mod, ONLY : ad_step2d USE ad_step3d_t_mod, ONLY : ad_step3d_t USE ad_step3d_uv_mod, ONLY : ad_step3d_uv USE dateclock_mod, ONLY : time_string USE frc_weak_mod, ONLY : frc_ADgather, frc_clear USE obs_read_mod, ONLY : obs_read USE omega_mod, ONLY : omega USE rho_eos_mod, ONLY : rho_eos USE set_depth_mod, ONLY : set_depth USE set_massflux_mod USE strings_mod, ONLY : FoundError ! implicit none ! ! Imported variable declarations. ! real(dp), intent(in) :: RunInterval ! ! Local variable declarations. ! logical :: backward = .TRUE. logical :: ad_advance logical :: DoNestLayer, Time_Step integer :: Nsteps, Rsteps integer :: ig, il, is, istep, ng, nl, tile integer :: my_iif, next_indx1 integer :: ks, kt ! real(r8) :: HalfDT, MaxDT, my_StepTime ! character (len=*), parameter :: MyFile = & & "ROMS/Adjoint/ad_main3d.F" ! !======================================================================= ! Time-step adjoint 3D primitive equations backwards. !======================================================================= ! ! Initialize. ! Time_Step=.TRUE. DoNestLayer=.TRUE. ! KERNEL_LOOP : DO WHILE (Time_Step) ! ! In nesting applications, the number of nesting layers (NestLayers) is ! used to facilitate refinement grids and composite/refinament grids ! combinations. Otherwise, the solution it is looped once for a single ! grid application (NestLayers = 1). ! nl=0 ! NEST_LAYER : DO WHILE (DoNestLayer) ! ! Determine number of time steps to compute in each nested grid layer ! based on the specified time interval (seconds), RunInterval. Non ! nesting applications have NestLayers=1. Notice that RunInterval is ! set in the calling driver. Its value may span the full period of the ! simulation, a multi-model coupling interval (RunInterval > ifac*dt), ! or just a single step (RunInterval=0). ! CALL ntimesteps (iADM, RunInterval, nl, Nsteps, Rsteps) IF (FoundError(exit_flag, NoError, 261, MyFile)) RETURN ! IF ((nl.le.0).or.(nl.gt.NestLayers)) EXIT ! ! Time-step governing equations for Nsteps. ! STEP_LOOP : DO istep=Nsteps,1,-1 ! ! Set time indices and time clock. ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) tdays(ng)=time(ng)*sec2day IF (step_counter(ng).eq.1) Time_Step=.FALSE. IF ((iic(ng).eq.ntstart(ng)).and.(ng.eq.Ngrids)) THEN ad_advance=.FALSE. ELSE ad_advance=.TRUE. END IF END DO ! !----------------------------------------------------------------------- ! Read in required data, if any, from input NetCDF files. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) CALL ad_get_data (ng) IF (FoundError(exit_flag, NoError, & & 300, MyFile)) RETURN END DO ! !----------------------------------------------------------------------- ! Process input data, if any: time interpolate between snapshots. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_set_data (ng, tile) CALL set_depth (ng, tile, iADM) END DO END DO IF (FoundError(exit_flag, NoError, 319, MyFile)) RETURN ! !----------------------------------------------------------------------- ! Compute BASIC STATE horizontal mass fluxes (Hz*u/n and Hz*v/m). !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=last_tile(ng),first_tile(ng),-1 CALL set_massflux (ng, tile, iADM) CALL rho_eos (ng, tile, iADM) END DO END DO ! ! Avoid time-stepping if additional delayed IO time-step. ! ADVANCE : IF (ad_advance) THEN ! !----------------------------------------------------------------------- ! Time-step adjoint tracer equations. !----------------------------------------------------------------------- ! ! Compute intermediate BASIC STATE mass fluxes (Huon,Hvom) for use in ! the adjoint horizontal advection (ad_step3d_t) and adjoint vertical ! velocity (ad_omega). ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL reset_massflux (ng, tile, iADM) ! intermediate END DO ! fluxes END DO ! ! Compute basic STATE omega vertical velocity with intermediate mass ! fluxes. Time-step adjoint tracer equations. ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=last_tile(ng),first_tile(ng),-1 CALL omega (ng, tile, iADM) ! BASIC STATE w-velocity CALL ad_step3d_t (ng, tile) END DO END DO ! !----------------------------------------------------------------------- ! Time-step adjoint vertical mixing turbulent equations and passive ! tracer source and sink terms, if applicable. Reinstate original ! BASIC STATE (Huon,Hvom). !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_omega (ng, tile, iADM) CALL set_massflux (ng, tile, iADM) ! BASIC STATE END DO ! mass fluxes END DO ! !----------------------------------------------------------------------- ! Time-step adjoint 3D equations. !----------------------------------------------------------------------- ! ! Reinstate original BASIC STATE omega vertical velocity. Time-step ! 3D adjoint momentum equations and couple with vertically integrated ! equations. ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=last_tile(ng),first_tile(ng),-1 CALL omega (ng, tile, iADM) ! BASIC STATE w-velocity CALL ad_step3d_uv (ng, tile) END DO END DO ! !----------------------------------------------------------------------- ! Adjoint of recompute depths and thicknesses using the new time ! filtered free-surface. This call was moved from "ad_step2d" to here. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=last_tile(ng),first_tile(ng),-1 CALL ad_set_depth (ng, tile, iADM) END DO END DO ! !----------------------------------------------------------------------- ! Solve adjoint vertically integrated primitive equations for the ! free-surface and barotropic momentum components. !----------------------------------------------------------------------- ! LOOP_2D : DO my_iif=MAXVAL(nfast)+1,1,-1 ! ! Corrector step - Apply 2D adjoint time-step corrector scheme. Notice ! ============== that there is not need for a corrector step during ! the auxiliary (nfast+1) time-step. ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) iif(ng)=my_iif IF (iif(ng).lt.(nfast(ng)+1)) THEN DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_step2d (ng, tile) END DO END IF ! ! Set time indices for adjoint predictor step. ! next_indx1=3-indx1(ng) IF (.not.PREDICTOR_2D_STEP(ng)) THEN PREDICTOR_2D_STEP(ng)=.TRUE. !^ knew(ng)=next_indx1 !^ kstp(ng)=3-knew(ng) !^ krhs(ng)=3 !^ kt=knew(ng) ks=kstp(ng) knew(ng)=krhs(ng) kstp(ng)=kt krhs(ng)=ks !^ IF (my_iif.lt.(nfast(ng)+1)) indx1(ng)=next_indx1 END IF END DO ! ! Predictor step - Advance adjoint barotropic equations using 2D ! ============== time-step predictor scheme. No actual time- ! stepping is performed during the auxiliary (nfast+1) time-step. ! It is needed to finalize the fast-time averaging of 2D fields, ! if any, and compute the new time-evolving depths. ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (my_iif.le.(nfast(ng)+1)) THEN DO tile=last_tile(ng),first_tile(ng),-1 CALL ad_step2d (ng, tile) END DO END IF ! ! Set time indices for next adjoint corrector step. The ! PREDICTOR_2D_STEP switch it is assumed to be false before the ! first time-step. ! IF (PREDICTOR_2D_STEP(ng).and. & & my_iif.le.(nfast(ng)+1)) THEN PREDICTOR_2D_STEP(ng)=.FALSE. !^ IF (iif(ng).eq.1) THEN !^ kstp(ng)=indx1(ng) !^ ELSE !^ kstp(ng)=3-indx1(ng) !^ END IF !^ knew(ng)=3 !^ krhs(ng)=indx1(ng) !^ ks=knew(ng) knew(ng)=krhs(ng) krhs(ng)=ks END IF END DO END DO LOOP_2D END IF ADVANCE ! !----------------------------------------------------------------------- ! If appropriate, read observation and model state at observation ! locations. Then, compute adjoint forcing terms due to observations. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) HalfDT=0.5_r8*dt(ng) IF (((time(ng)-HalfDT).le.ObsTime(ng)).and. & & (ObsTime(ng).lt.(time(ng)+HalfDT))) THEN ProcessObs(ng)=.TRUE. CALL obs_read (ng, iADM, backward) ELSE ProcessObs(ng)=.FALSE. END IF ! DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_htobs (ng, tile, iADM) END DO IF (FoundError(exit_flag, NoError, & & 719, MyFile)) RETURN END DO ! !----------------------------------------------------------------------- ! If appropriate, add representer coefficients (Beta hat) impulse ! forcing to adjoint solution. Read next impulse record, if available. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (ProcessObs(ng)) THEN DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_force_dual (ng, tile, kstp(ng), nstp(ng)) END DO END IF END DO ! !----------------------------------------------------------------------- ! Compute adjoint right-hand-side terms for 3D equations. If ! applicable, time-step turbulent mixing schemes. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_rhs3d (ng, tile) END DO END DO ! !----------------------------------------------------------------------- ! Set adjoint free-surface to it time-averaged value. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_set_zeta (ng, tile) END DO END DO ! !----------------------------------------------------------------------- ! Compute adjoint vertical mixing coefficients for momentum and ! tracers. Compute adjoint S-coordinate vertical velocity, ! diagnostically from horizontal mass divergence. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=last_tile(ng),first_tile(ng),-1 CALL ad_omega (ng, tile, iADM) END DO END DO ! !----------------------------------------------------------------------- ! Set adjoint fields for vertical boundary conditions. Process tidal ! forcing, if any. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_set_vbc (ng, tile) END DO END DO ! !----------------------------------------------------------------------- ! Compute adjoint density related fields and horizontal mass fluxes ! (Hz*u/n and Hz*v/m). If applicable, compute and report diagnostics ! and accumulate time-averaged adjoint fields which needs a ! irreversible loop in shared-memory jobs. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) DO tile=first_tile(ng),last_tile(ng),+1 ! irreversible CALL ad_rho_eos (ng, tile, iADM) CALL ad_set_massflux (ng, tile, iADM) CALL ad_diag (ng, tile) END DO END DO IF (FoundError(exit_flag, NoError, 974, MyFile)) RETURN ! !----------------------------------------------------------------------- ! Update 3D time-level indices. !----------------------------------------------------------------------- ! ! The original forward time-stepping indices are advanced as follows: ! ! nstp(ng)=1+MOD(iic(ng)-ntstart(ng),2) ! nnew(ng)=3-nstp(ng) ! nrhs(ng)=nstp(ng) ! ! This yields the only 2 permutations: ! ! nstp nnew nrhs ! 1 2 1 ! 2 1 2 ! ! With nstp=1, nnew=1 and nrhs=2 at time zero, this is equivalent to ! the following: ! ! nstp(ng)=nnew(ng) ! nnew(ng)=nrhs(ng) ! nrhs(ng)=nstp(ng) ! ! The adjoint of this is as follows: ! ! nstp(ng)=nrhs(ng) ! nrhs(ng)=nnew(ng) ! nnew(ng)=nstp(ng) ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (iic(ng).ne.ntend(ng)) THEN nrhs(ng)=nnew(ng) nnew(ng)=nstp(ng) nstp(ng)=nrhs(ng) END IF END DO ! !----------------------------------------------------------------------- ! If not a restart, initialize all time levels and compute other ! initial fields. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (iic(ng).eq.ntend(ng)) THEN ! ! Adjoint of initialize other state variables. ! DO tile=last_tile(ng),first_tile(ng),-1 CALL ad_ini_fields (ng, tile, iADM) END DO ! ! Adjoint of initialize free-surface and compute initial level ! thicknesses and depths. ! DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_set_depth (ng, tile, iADM) CALL ad_ini_zeta (ng, tile, iADM) END DO END IF END DO ! !----------------------------------------------------------------------- ! Gather weak constraint forcing to storage arrays. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (iic(ng).ne.ntstart(ng)) THEN DO tile=first_tile(ng),last_tile(ng),+1 CALL ad_set_depth (ng, tile, iADM) CALL frc_ADgather (ng, tile) END DO END IF END DO ! !----------------------------------------------------------------------- ! If appropriate, write out fields into output NetCDF files. Notice ! that IO data is written in delayed and serial mode. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) CALL ad_output (ng) IF (FoundError(exit_flag, NoError, & & 1178, MyFile)) RETURN END DO ! !----------------------------------------------------------------------- ! Copy storage arrays index 1 into index 2, and clear index 1. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) IF (MOD(iic(ng)-1,nADJ(ng)).eq.0) THEN DO tile=first_tile(ng),last_tile(ng),+1 CALL frc_clear (ng, tile) END DO END IF END DO ! !----------------------------------------------------------------------- ! Set adjoint time indices and time clock. !----------------------------------------------------------------------- ! DO ig=1,GridsInLayer(nl) ng=GridNumber(ig,nl) iic(ng)=iic(ng)-1 step_counter(ng)=step_counter(ng)-1 time(ng)=time(ng)-dt(ng) CALL time_string (time(ng), time_code(ng)) END DO END DO STEP_LOOP END DO NEST_LAYER END DO KERNEL_LOOP ! RETURN END SUBROUTINE ad_main3d