MODULE roms_kernel_mod ! !git $Id$ !svn $Id: optobs_roms.h 1166 2023-05-17 20:11:58Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group W. G. Zhang ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! ROMS/TOMS Optimal Observation Driver: ! ! ! ! This driver computes the adjoint sensitivity of a function or ! ! index, J, in terms of space and/or time integrals of the model ! ! state, S(zeta,u,v,T,...). Small changes, dS, in S will lead to ! ! changes dJ in J: ! ! ! ! dJ = (dJ/dzeta) dzeta + (dJ/du) du + (dJ/dv) dv + (dJ/dt) dT ... ! ! ! ! and ! ! ! ! dJ/dS = transpose(R) S ! ! ! ! where transpose(R) is the adjoint propagator. It implies that ! ! the sensitivity for ALL variables, parameters, and space-time ! ! points can be computed from a single integration of the adjoint ! ! model. ! ! ! ! These routines control the initialization, time-stepping, and ! ! finalization of ROMS/TOMS model following ESMF conventions: ! ! ! ! ROMS_initialize ! ! ROMS_run ! ! ROMS_finalize ! ! ! ! Reference: ! ! ! ! ! !======================================================================= ! USE mod_param USE mod_parallel USE mod_arrays USE mod_fourdvar USE mod_iounits USE mod_ncparam USE mod_ocean USE mod_scalars USE mod_stepping ! USE close_io_mod, ONLY : close_inp, close_out USE convolve_mod, ONLY : error_covariance USE def_norm_mod, ONLY : def_norm USE inp_par_mod, ONLY : inp_par USE get_state_mod, ONLY : get_state #ifdef MCT_LIB # ifdef ATM_COUPLING USE mct_coupler_mod, ONLY : initialize_ocn2atm_coupling # endif # ifdef WAV_COUPLING USE mct_coupler_mod, ONLY : initialize_ocn2wav_coupling # endif #endif USE normalization_mod, ONLY : normalization USE strings_mod, ONLY : FoundError USE tl_def_ini_mod, ONLY : tl_def_ini USE wrt_rst_mod, ONLY : wrt_rst #if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC USE zeta_balance_mod, ONLY : balance_ref, biconj #endif ! implicit none ! PUBLIC :: ROMS_initialize PUBLIC :: ROMS_run PUBLIC :: ROMS_finalize ! CONTAINS ! SUBROUTINE ROMS_initialize (first, mpiCOMM) ! !======================================================================= ! ! ! This routine allocates and initializes ROMS/TOMS state variables ! ! and internal and external parameters. ! ! ! !======================================================================= ! ! Imported variable declarations. ! logical, intent(inout) :: first ! integer, intent(in), optional :: mpiCOMM ! ! Local variable declarations. ! logical :: allocate_vars = .TRUE. ! #ifdef DISTRIBUTE integer :: MyError, MySize #endif integer :: STDrec, Tindex integer :: chunk_size, ng, thread #ifdef _OPENMP integer :: my_threadnum #endif ! character (len=*), parameter :: MyFile = & & __FILE__//", ROMS_initialize" #ifdef DISTRIBUTE ! !----------------------------------------------------------------------- ! Set distribute-memory (mpi) world communictor. !----------------------------------------------------------------------- ! IF (PRESENT(mpiCOMM)) THEN OCN_COMM_WORLD=mpiCOMM ELSE OCN_COMM_WORLD=MPI_COMM_WORLD END IF CALL mpi_comm_rank (OCN_COMM_WORLD, MyRank, MyError) CALL mpi_comm_size (OCN_COMM_WORLD, MySize, MyError) #endif ! !----------------------------------------------------------------------- ! On first pass, initialize model parameters a variables for all ! nested/composed grids. Notice that the logical switch "first" ! is used to allow multiple calls to this routine during ensemble ! configurations. !----------------------------------------------------------------------- ! IF (first) THEN first=.FALSE. ! ! Initialize parallel control switches. These scalars switches are ! independent from standard input parameters. ! CALL initialize_parallel ! ! Read in model tunable parameters from standard input. Allocate and ! initialize variables in several modules after the number of nested ! grids and dimension parameters are known. ! CALL inp_par (iADM) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN ! ! Set domain decomposition tile partition range. This range is ! computed only once since the "first_tile" and "last_tile" values ! are private for each parallel thread/node. ! #if defined _OPENMP MyThread=my_threadnum() #elif defined DISTRIBUTE MyThread=MyRank #else MyThread=0 #endif DO ng=1,Ngrids chunk_size=(NtileX(ng)*NtileE(ng)+numthreads-1)/numthreads first_tile(ng)=MyThread*chunk_size last_tile (ng)=first_tile(ng)+chunk_size-1 END DO ! ! Initialize internal wall clocks. Notice that the timings does not ! includes processing standard input because several parameters are ! needed to allocate clock variables. ! IF (Master) THEN WRITE (stdout,10) 10 FORMAT (/,' Process Information:',/) END IF ! DO ng=1,Ngrids DO thread=THREAD_RANGE CALL wclock_on (ng, iADM, 0, __LINE__, MyFile) END DO END DO ! ! Allocate and initialize all model state arrays. ! CALL ROMS_allocate_arrays (allocate_vars) CALL ROMS_initialize_arrays IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END IF #if defined MCT_LIB && (defined ATM_COUPLING || defined WAV_COUPLING) ! !----------------------------------------------------------------------- ! Initialize coupling streams between model(s). !----------------------------------------------------------------------- ! DO ng=1,Ngrids # ifdef ATM_COUPLING CALL initialize_ocn2atm_coupling (ng, MyRank) # endif # ifdef WAV_COUPLING CALL initialize_ocn2wav_coupling (ng, MyRank) # endif END DO #endif ! !----------------------------------------------------------------------- ! Read in standard deviation factors for error covariance. !----------------------------------------------------------------------- ! ! Initial conditions standard deviation. They are loaded in Tindex=1 ! of the e_var(...,Tindex) state variables. ! STDrec=1 Tindex=1 DO ng=1,Ngrids CALL get_state (ng, 10, 10, STD(1,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO ! ! Model error standard deviation. They are loaded in Tindex=2 ! of the e_var(...,Tindex) state variables. ! STDrec=1 Tindex=2 DO ng=1,Ngrids IF (NSA.eq.2) THEN CALL get_state (ng, 11, 11, STD(2,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END IF END DO #ifdef ADJUST_BOUNDARY ! ! Open boundary conditions standard deviation. ! STDrec=1 Tindex=1 DO ng=1,Ngrids CALL get_state (ng, 12, 12, STD(3,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO #endif #if defined ADJUST_WSTRESS || defined ADJUST_STFLUX ! ! Surface forcing standard deviation. ! STDrec=1 Tindex=1 DO ng=1,Ngrids CALL get_state (ng, 13, 13, STD(4,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO #endif ! RETURN END SUBROUTINE ROMS_initialize ! SUBROUTINE ROMS_run (RunInterval) ! !======================================================================= ! ! ! This routine computes the adjoint sensitivity analysis, dJ/dS, ! ! to the specified functional J. The sensitivity masking arrays ! ! Rscope, Uscope, and Vscope are used to evaluate the functional ! ! in the desired spatial area. ! ! ! !======================================================================= ! ! Imported variable declarations ! real(dp), intent(in) :: RunInterval ! seconds ! ! Local variable declarations. ! logical :: Lposterior ! integer :: ng, tile integer :: Lbck, Lini, Rec, Rec1, Rec2 integer :: NRMrec ! real (r8) :: str_day, end_day ! character (len=6) :: driver character (len=*), parameter :: MyFile = & & __FILE__//", ROMS_run" ! !======================================================================= ! Run model for all nested grids, if any. !======================================================================= ! ! Initialize relevant parameters. ! DO ng=1,Ngrids Lold(ng)=1 ! old minimization time index Lnew(ng)=2 ! new minimization time index LreadFWD(ng)=.TRUE. END DO Lini=1 ! NLM initial conditions record in INI Lbck=1 ! background record in INI Rec1=1 Rec2=2 driver='optobs' #ifdef FORWARD_FLUXES ! ! Set the BLK structure to contain the nonlinear model surface fluxes ! needed by the tangent linear and adjoint models. Also, set switches ! to process that structure in routine "check_multifile". Notice that ! it is possible to split the solution into multiple NetCDF files to ! reduce their size. ! ! The switch LreadFRC is deactivated because all the atmospheric ! forcing, including shortwave radiation, is read from the NLM ! surface fluxes or is assigned during ESM coupling. Such fluxes ! are available from the QCK structure. There is no need for reading ! and processing from the FRC structure input forcing-files. ! CALL edit_multifile ('QCK2BLK') IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN DO ng=1,Ngrids LreadBLK(ng)=.TRUE. LreadFRC(ng)=.FALSE. END DO #endif #ifdef BALANCE_OPERATOR ! ! Set NLM model background trajectory to process in the balance ! operator. ! DO ng=1,Ngrids INI(ng)%name=FWD(ng)%name END DO # ifdef ZETA_ELLIPTIC ! ! Compute the reference zeta and biconjugate gradient arrays ! required for the balance of free surface. ! IF (balance(isFsur)) THEN DO ng=1,Ngrids CALL get_state (ng, iNLM, 2, INI(ng), Lini, Lini) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN ! DO tile=first_tile(ng),last_tile(ng),+1 CALL balance_ref (ng, tile, Lini) CALL biconj (ng, tile, iNLM, Lini) END DO wrtZetaRef(ng)=.TRUE. END DO END IF # endif #endif ! !----------------------------------------------------------------------- ! Initialize adjoint model and define sensitivity functional. !----------------------------------------------------------------------- ! DO ng=1,Ngrids Lstiffness=.FALSE. CALL ad_initial (ng) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO ! ! Activate adjoint output. ! DO ng=1,Ngrids LdefADJ(ng)=.TRUE. LwrtADJ(ng)=.TRUE. LcycleADJ(ng)=.FALSE. END DO ! ! Define tangent linear initial conditions file. ! DO ng=1,Ngrids LdefITL(ng)=.TRUE. CALL tl_def_ini (ng) LdefITL(ng)=.FALSE. IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO ! !----------------------------------------------------------------------- ! Compute or read in background-error correlations normalization ! factors. !----------------------------------------------------------------------- ! ! If computing, write out factors to NetCDF. This is an expensive ! computation and needs to be computed once for a particular ! application grid. ! DO ng=1,Ngrids IF (ANY(LwrtNRM(:,ng))) THEN CALL def_norm (ng, iNLM, 1) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN IF (NSA.eq.2) THEN CALL def_norm (ng, iNLM, 2) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END IF #ifdef ADJUST_BOUNDARY CALL def_norm (ng, iNLM, 3) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN #endif #if defined ADJUST_WSTRESS || defined ADJUST_STFLUX CALL def_norm (ng, iNLM, 4) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN #endif DO tile=first_tile(ng),last_tile(ng),+1 CALL normalization (ng, tile, 2) END DO LdefNRM(1:4,ng)=.FALSE. LwrtNRM(1:4,ng)=.FALSE. ELSE NRMrec=1 CALL get_state (ng, 14, 14, NRM(1,ng), NRMrec, 1) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN IF (NSA.eq.2) THEN CALL get_state (ng, 15, 15, NRM(2,ng), NRMrec, 2) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END IF #ifdef ADJUST_BOUNDARY CALL get_state (ng, 16, 16, NRM(3,ng), NRMrec, 1) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN #endif #if defined ADJUST_WSTRESS || defined ADJUST_STFLUX CALL get_state (ng, 17, 17, NRM(3,ng), NRMrec, 1) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN #endif END IF END DO ! !------------------------------------------------------------------------ ! Compute the gradient or index, dJ/dS, of the sensitivity functional. !------------------------------------------------------------------------ ! DO ng=1,Ngrids str_day=tdays(ng) end_day=str_day-ntimes(ng)*dt(ng)*sec2day IF ((DstrS(ng).eq.0.0_r8).and.(DendS(ng).eq.0.0_r8)) THEN DstrS(ng)=end_day DendS(ng)=str_day END IF IF (Master) THEN WRITE (stdout,10) 'AD', DendS(ng), DstrS(ng) END IF IF ((DstrS(ng).gt.str_day).or.(DstrS(ng).lt.end_day)) THEN IF (Master) WRITE (stdout,20) 'DstrS = ', DstrS(ng), & & end_day, str_day exit_flag=7 RETURN END IF IF ((DendS(ng).gt.str_day).or.(DendS(ng).lt.end_day)) THEN IF (Master) WRITE (stdout,20) 'DendS = ', DendS(ng), & & end_day, str_day exit_flag=7 RETURN END IF END DO ! ! Time-step adjoint model. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,30) 'AD', ng, ntstart(ng), ntend(ng) END IF END DO ! #ifdef SOLVE3D CALL ad_main3d (RunInterval) #else CALL ad_main2d (RunInterval) #endif IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN ! ! Convolve adjoint trajectory with error covariances and write ! TLM initial conditions. ! Lposterior=.FALSE. CALL error_covariance (iTLM, driver, -1, -1, & & Lbck, Lini, Lold, Lnew, & & Rec1, Rec2, Lposterior) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN ! !----------------------------------------------------------------------- ! Run tangent linear model for all nested grids, if any. !----------------------------------------------------------------------- ! ! Initialize tangent linear model. ! DO ng=1,Ngrids CALL tl_initial (ng) IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN END DO ! ! Activate tangent linear output. ! DO ng=1,Ngrids LdefTLM(ng)=.TRUE. LwrtTLM(ng)=.TRUE. LcycleTLM(ng)=.FALSE. END DO ! ! Time-step tangent linear model. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,30) 'TL', ng, ntstart(ng), ntend(ng) END IF END DO ! #ifdef SOLVE3D CALL tl_main3d (RunInterval) #else CALL tl_main2d (RunInterval) #endif IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN ! 10 FORMAT (/,'AD ROMS/TOMS: adjoint forcing time range: ', & & f12.4,' - ',f12.4 ,/) 20 FORMAT (/,' Out of range adjoint forcing time, ',a,f12.4,/, & & ' It must be between ',f12.4,' and ',f12.4) 30 FORMAT (/,1x,a,1x,'ROMS/TOMS: started time-stepping:', & & ' (Grid: ',i2.2,' TimeSteps: ',i0,' - ',i0,')',/) ! RETURN END SUBROUTINE ROMS_run ! SUBROUTINE ROMS_finalize ! !======================================================================= ! ! ! This routine terminates ROMS/TOMS nonlinear and adjoint models ! ! execution. ! ! ! !======================================================================= ! ! Local variable declarations. ! integer :: Fcount, ng, thread ! character (len=*), parameter :: MyFile = & & __FILE__//", ROMS_finalize" ! !----------------------------------------------------------------------- ! If blowing-up, save latest model state into RESTART NetCDF file. !----------------------------------------------------------------------- ! ! If cycling restart records, write solution into the next record. ! DO ng=1,Ngrids IF (LwrtRST(ng).and.(exit_flag.eq.1)) THEN IF (Master) WRITE (stdout,10) 10 FORMAT (/,' Blowing-up: Saving latest model state into ', & & ' RESTART file',/) Fcount=RST(ng)%load IF (LcycleRST(ng).and.(RST(ng)%Nrec(Fcount).ge.2)) THEN RST(ng)%Rindex=2 LcycleRST(ng)=.FALSE. END IF blowup=exit_flag exit_flag=NoError #ifdef DISTRIBUTE CALL wrt_rst (ng, MyRank) #else CALL wrt_rst (ng, -1) #endif END IF END DO ! !----------------------------------------------------------------------- ! Stop model and time profiling clocks, report memory requirements, and ! close output NetCDF files. !----------------------------------------------------------------------- ! ! Stop time clocks. ! IF (Master) THEN WRITE (stdout,20) 20 FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/) END IF ! DO ng=1,Ngrids DO thread=THREAD_RANGE CALL wclock_off (ng, iADM, 0, __LINE__, MyFile) END DO END DO ! ! Report dynamic memory and automatic memory requirements. ! CALL memory ! ! Close IO files. ! DO ng=1,Ngrids CALL close_inp (ng, iADM) END DO CALL close_out ! RETURN END SUBROUTINE ROMS_finalize END MODULE roms_kernel_mod