MODULE normalization_mod ! !git $Id$ !svn $Id: normalization.F 1189 2023-08-15 21:26:58Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! This routine computes the background covariance, B, normalization ! ! factors using the exact or approximated approach. These factors ! ! ensure that the diagonal elements of B are equal to unity. Notice ! ! that in applications with land/sea masking, it will produce large ! ! changes in the covariance structures near the boundary. ! ! ! ! The exact method is very expensive: the normalization factors are ! ! computed by perturbing each model grid cell with a delta function ! ! scaled by the area (2D factors) or volume (3D factors), and then ! ! convolving with the squared-root adjoint and tangent diffusion ! ! operators. ! ! ! ! The approximated method is cheaper: the normalization factors are ! ! computed using the randomization approach of Fisher and Courtier ! ! (1995). The factors are initialized with randon numbers having a ! ! uniform distribution (drawn from a normal distribution with zero ! ! mean and unity variance). Then, they scaled the inverse, squared ! ! root cell area (2D factors) or volume (3D factors) and convoluted ! ! with the squared-root adjoint and tangent diffuse operators over ! ! a specified number of iterations, Nrandom. ! ! ! ! References: ! ! ! ! Fisher, M. and. P. Courtier, 1995: Estimating the covariance ! ! matrices of analysis and forecast error in variational data ! ! assimilation, ECMWF Technical Memo N. 220, ECMWF, Reading, ! ! UK. ! ! ! ! www.ecmwf.int/publications/library/ecpublications/_pdf/tm/ ! ! 001-300/tm220.pdf ! ! ! ! Weaver, A. and P. Courtier, 2001: Correlation modeling on the ! ! sphere using a generalized diffusion equation, Q.J.R. Meteo. ! ! Soc, 127, 1815-1846. ! ! ! !======================================================================! ! USE mod_kinds USE mod_param USE mod_parallel USE mod_grid USE mod_fourdvar USE mod_iounits USE mod_ncparam USE mod_netcdf USE mod_mixing USE mod_ocean USE mod_scalars USE mod_stepping ! USE ad_conv_2d_mod USE ad_conv_3d_mod USE bc_2d_mod USE bc_3d_mod USE mp_exchange_mod USE set_depth_mod USE tl_conv_2d_mod USE tl_conv_3d_mod USE white_noise_mod ! USE distribute_mod, ONLY : mp_reduce USE nf_fwrite2d_mod, ONLY : nf_fwrite2d USE nf_fwrite3d_mod, ONLY : nf_fwrite3d USE strings_mod, ONLY : FoundError ! implicit none ! PUBLIC :: normalization PRIVATE :: normalization_tile PRIVATE :: randomization_tile PRIVATE :: wrt_norm2d_nf90 PRIVATE :: wrt_norm3d_nf90 ! CONTAINS ! !*********************************************************************** SUBROUTINE normalization (ng, tile, ifac) !*********************************************************************** ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, ifac ! ! Local variable declarations. ! integer :: IminS, ImaxS, JminS, JmaxS integer :: LBi, UBi, LBj, UBj, LBij, UBij ! ! Set horizontal starting and ending indices for automatic private ! storage arrays. ! IminS=BOUNDS(ng)%Istr(tile)-3 ImaxS=BOUNDS(ng)%Iend(tile)+3 JminS=BOUNDS(ng)%Jstr(tile)-3 JmaxS=BOUNDS(ng)%Jend(tile)+3 ! ! Determine array lower and upper bounds in the I- and J-directions. ! LBi=BOUNDS(ng)%LBi(tile) UBi=BOUNDS(ng)%UBi(tile) LBj=BOUNDS(ng)%LBj(tile) UBj=BOUNDS(ng)%UBj(tile) ! ! Set array lower and upper bounds for MIN(I,J) directions and ! MAX(I,J) directions. ! LBij=BOUNDS(ng)%LBij UBij=BOUNDS(ng)%UBij ! ! Compute background error covariance normalization factors using ! the very expensive exact method. ! IF (Nmethod(ng).eq.0) THEN CALL normalization_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), ifac, & & GRID(ng) % pm, & & GRID(ng) % om_p, & & GRID(ng) % om_r, & & GRID(ng) % om_u, & & GRID(ng) % om_v, & & GRID(ng) % pn, & & GRID(ng) % on_p, & & GRID(ng) % on_r, & & GRID(ng) % on_u, & & GRID(ng) % on_v, & & GRID(ng) % pmon_p, & & GRID(ng) % pmon_r, & & GRID(ng) % pmon_u, & & GRID(ng) % pnom_p, & & GRID(ng) % pnom_r, & & GRID(ng) % pnom_v, & & GRID(ng) % pmask, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & & GRID(ng) % h, & & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & & MIXING(ng) % Kh, & & MIXING(ng) % Kv, & & OCEAN(ng) % b_t, & & OCEAN(ng) % b_u, & & OCEAN(ng) % b_v, & & OCEAN(ng) % b_zeta, & & OCEAN(ng) % b_ubar, & & OCEAN(ng) % b_vbar) ! ! Compute background error covariance normalization factors using ! the approximated randomization method. ! ELSE IF (Nmethod(ng).eq.1) THEN CALL randomization_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), ifac, & & GRID(ng) % pm, & & GRID(ng) % om_p, & & GRID(ng) % om_r, & & GRID(ng) % om_u, & & GRID(ng) % om_v, & & GRID(ng) % pn, & & GRID(ng) % on_p, & & GRID(ng) % on_r, & & GRID(ng) % on_u, & & GRID(ng) % on_v, & & GRID(ng) % pmon_p, & & GRID(ng) % pmon_r, & & GRID(ng) % pmon_u, & & GRID(ng) % pnom_p, & & GRID(ng) % pnom_r, & & GRID(ng) % pnom_v, & & GRID(ng) % pmask, & & GRID(ng) % rmask, & & GRID(ng) % umask, & & GRID(ng) % vmask, & & GRID(ng) % h, & & GRID(ng) % Hz, & & GRID(ng) % z_r, & & GRID(ng) % z_w, & & MIXING(ng) % Kh, & & MIXING(ng) % Kv, & & OCEAN(ng) % b_t, & & OCEAN(ng) % b_u, & & OCEAN(ng) % b_v, & & OCEAN(ng) % b_zeta, & & OCEAN(ng) % b_ubar, & & OCEAN(ng) % b_vbar) END IF ! RETURN END SUBROUTINE normalization ! !*********************************************************************** SUBROUTINE normalization_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, ifac, & & pm, om_p, om_r, om_u, om_v, & & pn, on_p, on_r, on_u, on_v, & & pmon_p, pmon_r, pmon_u, & & pnom_p, pnom_r, pnom_v, & & pmask, rmask, & & umask, vmask, & & h, & & Hz, z_r, z_w, & & Kh, & & Kv, & & VnormR, VnormU, VnormV, & & HnormR, HnormU, HnormV) !*********************************************************************** ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew, ifac ! real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: om_p(LBi:,LBj:) real(r8), intent(in) :: om_r(LBi:,LBj:) real(r8), intent(in) :: om_u(LBi:,LBj:) real(r8), intent(in) :: om_v(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) real(r8), intent(in) :: on_p(LBi:,LBj:) real(r8), intent(in) :: on_r(LBi:,LBj:) real(r8), intent(in) :: on_u(LBi:,LBj:) real(r8), intent(in) :: on_v(LBi:,LBj:) real(r8), intent(in) :: pmon_p(LBi:,LBj:) real(r8), intent(in) :: pmon_r(LBi:,LBj:) real(r8), intent(in) :: pmon_u(LBi:,LBj:) real(r8), intent(in) :: pnom_p(LBi:,LBj:) real(r8), intent(in) :: pnom_r(LBi:,LBj:) real(r8), intent(in) :: pnom_v(LBi:,LBj:) real(r8), intent(in) :: pmask(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) real(r8), intent(in) :: Kh(LBi:,LBj:) real(r8), intent(in) :: Kv(LBi:,LBj:,0:) real(r8), intent(inout) :: h(LBi:,LBj:) real(r8), intent(out) :: VnormR(LBi:,LBj:,:,:,:) real(r8), intent(out) :: VnormU(LBi:,LBj:,:,:) real(r8), intent(out) :: VnormV(LBi:,LBj:,:,:) real(r8), intent(out) :: HnormR(LBi:,LBj:,:) real(r8), intent(out) :: HnormU(LBi:,LBj:,:) real(r8), intent(out) :: HnormV(LBi:,LBj:,:) real(r8), intent(out) :: Hz(LBi:,LBj:,:) real(r8), intent(out) :: z_r(LBi:,LBj:,:) real(r8), intent(out) :: z_w(LBi:,LBj:,0:) ! ! Local variable declarations. ! logical :: Ldiffer, Lsame ! integer :: Imin, Imax, Jmin, Jmax integer :: i, ic, ifile, is, j, jc, rec integer :: NSUB integer :: UBt, itrc, k, kc, ntrc real(dp) :: my_time real(r8) :: cff, compute real(r8) :: my_dot, Gdotp ! real(r8), dimension(LBi:UBi,LBj:UBj) :: A2d real(r8), dimension(LBi:UBi,LBj:UBj) :: Hscale real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: A3d real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: Vscale ! character (len=3 ) :: op_handle character (len=40 ) :: Text character (len=256) :: ncname character (len=*), parameter :: MyFile = & & "ROMS/Utility/normalization.F"//", normalization_tile" ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU integer :: Iend, IendB, IendP, IendR, IendT integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV integer :: Jend, JendB, JendP, JendR, JendT integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1 integer :: Iendp1, Iendp2, Iendp2i, Iendp3 integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1 integer :: Jendp1, Jendp2, Jendp2i, Jendp3 ! Istr =BOUNDS(ng) % Istr (tile) IstrB =BOUNDS(ng) % IstrB (tile) IstrM =BOUNDS(ng) % IstrM (tile) IstrP =BOUNDS(ng) % IstrP (tile) IstrR =BOUNDS(ng) % IstrR (tile) IstrT =BOUNDS(ng) % IstrT (tile) IstrU =BOUNDS(ng) % IstrU (tile) Iend =BOUNDS(ng) % Iend (tile) IendB =BOUNDS(ng) % IendB (tile) IendP =BOUNDS(ng) % IendP (tile) IendR =BOUNDS(ng) % IendR (tile) IendT =BOUNDS(ng) % IendT (tile) Jstr =BOUNDS(ng) % Jstr (tile) JstrB =BOUNDS(ng) % JstrB (tile) JstrM =BOUNDS(ng) % JstrM (tile) JstrP =BOUNDS(ng) % JstrP (tile) JstrR =BOUNDS(ng) % JstrR (tile) JstrT =BOUNDS(ng) % JstrT (tile) JstrV =BOUNDS(ng) % JstrV (tile) Jend =BOUNDS(ng) % Jend (tile) JendB =BOUNDS(ng) % JendB (tile) JendP =BOUNDS(ng) % JendP (tile) JendR =BOUNDS(ng) % JendR (tile) JendT =BOUNDS(ng) % JendT (tile) ! Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3 Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2 Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1 IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2 IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1 Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1 Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2 Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3 Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3 Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2 Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1 JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2 JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1 Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1 Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2 Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3 ! SourceFile=MyFile my_time=tdays(ng)*day2sec ! !----------------------------------------------------------------------- ! Compute time invariant depths (use zero free-surface). !----------------------------------------------------------------------- ! DO i=LBi,UBi DO j=LBj,UBj A2d(i,j)=0.0_r8 END DO END DO CALL set_depth_tile (ng, tile, iNLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & h, & & A2d, & & Hz, z_r, z_w) ! !----------------------------------------------------------------------- ! Compute initial conditions and model error covariance, B, ! normalization factors using the exact method. It involves ! computing the filter variance (convolution) at each point ! independenly. That is, each point is perturbed with a delta ! function, scaled by the inverse squared root of the area (2D) ! or volume (3D), and then convoluted. !----------------------------------------------------------------------- ! IF (Master) WRITE (stdout,10) FILE_LOOP : DO ifile=1,NSA IF (LwrtNRM(ifile,ng)) THEN IF (ifile.eq.1) THEN Text='initial conditions' ELSE IF (ifile.eq.2) THEN Text='model' END IF ! ! Set time record index to write in normalization NetCDF file. ! ncname=NRM(ifile,ng)%name NRM(ifile,ng)%Rindex=NRM(ifile,ng)%Rindex+1 NRM(ifile,ng)%Nrec=NRM(ifile,ng)%Nrec+1 ! ! Write out model time (s). ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL netcdf_put_fvar (ng, iTLM, ncname, & & Vname(1,idtime), my_time, & & start = (/NRM(ifile,ng)%Rindex/), & & total = (/1/), & & ncid = NRM(ifile,ng)%ncid, & & varid = NRM(ifile,ng)%Vid(idtime)) END SELECT IF (FoundError(exit_flag, NoError, 618, MyFile)) RETURN ! ! 2D norm at RHO-points. ! IF (Cnorm(ifile,isFsur)) THEN Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at RHO-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrT,IendT Hscale(i,j)=1.0_r8/SQRT(om_r(i,j)*on_r(i,j)) END DO END DO DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (rmask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO j=LBj,UBj DO i=LBi,UBi A2d(i,j)=0.0_r8 END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A2d(ic,jc)=1.0_r8 END IF CALL ad_conv_r2d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isFsur)/ifac, & & DTsizeH(ifile,isFsur), & & Kh, & & pm, pn, pmon_u, pnom_v, & & rmask, umask, vmask, & & A2d) DO j=JstrT,JendT DO i=IstrT,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO ! my_dot=0.0_r8 DO j=JstrT,JendT DO i=IstrT,IendT my_dot=my_dot+A2d(i,j)*A2d(i,j) END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (R2_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (R2_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN HnormR(ic,jc,ifile)=cff END IF END DO END DO CALL dabc_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormR(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormR(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idFsur, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idFsur), & & NRM(ifile,ng)%Rindex, & & rmask, & & HnormR(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 766, MyFile)) RETURN END IF ! ! 2D norm at U-points. ! IF (Cnorm(ifile,isUbar)) THEN IF (EWperiodic(ng)) THEN Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) ELSE Imin=2 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) END IF IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at U-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrP,IendT Hscale(i,j)=1.0_r8/SQRT(om_u(i,j)*on_u(i,j)) END DO END DO DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (umask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO j=LBj,UBj DO i=LBi,UBi A2d(i,j)=0.0_r8 END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A2d(ic,jc)=1.0_r8 END IF CALL ad_conv_u2d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isUbar)/ifac, & & DTsizeH(ifile,isUbar), & & Kh, & & pm, pn, pmon_r, pnom_p, & & umask, pmask, & & A2d) DO j=JstrT,JendT DO i=IstrP,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO ! my_dot=0.0_r8 DO j=JstrT,JendT DO i=IstrP,IendT my_dot=my_dot+A2d(i,j)*A2d(i,j) END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (U2_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (U2_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN HnormU(ic,jc,ifile)=cff END IF END DO END DO CALL dabc_u2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormU(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormU(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idUbar, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idUbar), & & NRM(ifile,ng)%Rindex, & & umask, & & HnormU(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 923, MyFile)) RETURN END IF ! ! 2D norm at V-points. ! IF (Cnorm(ifile,isVbar)) THEN IF (NSperiodic(ng)) THEN Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) ELSE Imin=1 Imax=Lm(ng) Jmin=2 Jmax=Mm(ng) END IF IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at V-points' CALL my_flush (stdout) END IF DO j=JstrP,JendT DO i=IstrT,IendT Hscale(i,j)=1.0_r8/SQRT(om_v(i,j)*on_v(i,j)) END DO END DO DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (vmask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO j=LBj,UBj DO i=LBi,UBi A2d(i,j)=0.0_r8 END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A2d(ic,jc)=1.0_r8 END IF CALL ad_conv_v2d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isVbar)/ifac, & & DTsizeH(ifile,isVbar), & & Kh, & & pm, pn, pmon_p, pnom_r, & & vmask, pmask, & & A2d) DO j=JstrP,JendT DO i=IstrT,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO ! my_dot=0.0_r8 DO j=JstrP,JendT DO i=IstrT,IendT my_dot=my_dot+A2d(i,j)*A2d(i,j) END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (V2_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (V2_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN HnormV(ic,jc,ifile)=cff END IF END DO END DO CALL dabc_v2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormV(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormV(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idVbar, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idVbar), & & NRM(ifile,ng)%Rindex, & & vmask, & & HnormV(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 1079, MyFile)) RETURN END IF ! ! 3D norm U-points. ! IF (Cnorm(ifile,isUvel)) THEN IF (EWperiodic(ng)) THEN Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) ELSE Imin=2 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) END IF IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at U-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrP,IendT cff=om_u(i,j)*on_u(i,j)*0.5_r8 DO k=1,N(ng) Vscale(i,j,k)=1.0_r8/SQRT(cff*(Hz(i-1,j,k)+Hz(i,j,k))) END DO END DO END DO DO kc=1,N(ng) DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (umask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO k=1,N(ng) DO j=LBj,UBj DO i=LBi,UBi A3d(i,j,k)=0.0_r8 END DO END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A3d(ic,jc,kc)=1.0_r8 END IF CALL ad_conv_u3d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isUvel)/ifac, & & NVsteps(ifile,isUvel)/ifac, & & DTsizeH(ifile,isUvel), & & DTsizeV(ifile,isUvel), & & Kh, Kv, & & pm, pn, & & pmon_r, pnom_p, & & umask, pmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO ! my_dot=0.0_r8 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT my_dot=my_dot+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (R3_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (R3_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN VnormU(ic,jc,kc,ifile)=cff END IF END DO END DO END DO CALL dabc_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormU(:,:,:,ifile)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormU(:,:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idUvel, NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idUvel), & & NRM(ifile,ng)%Rindex, & & umask, & & VnormU(:,:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 1261, MyFile)) RETURN END IF ! ! 3D norm at V-points. ! IF (Cnorm(ifile,isVvel)) THEN IF (NSperiodic(ng)) THEN Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) ELSE Imin=1 Imax=Lm(ng) Jmin=2 Jmax=Mm(ng) END IF IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at V-points' CALL my_flush (stdout) END IF DO j=JstrP,JendT DO i=IstrT,IendT cff=om_v(i,j)*on_v(i,j)*0.5_r8 DO k=1,N(ng) Vscale(i,j,k)=1.0_r8/SQRT(cff*(Hz(i,j-1,k)+Hz(i,j,k))) END DO END DO END DO DO kc=1,N(ng) DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (vmask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO k=1,N(ng) DO j=LBj,UBj DO i=LBi,UBi A3d(i,j,k)=0.0_r8 END DO END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A3d(ic,jc,kc)=1.0_r8 END IF CALL ad_conv_v3d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isVvel)/ifac, & & NVsteps(ifile,isVvel)/ifac, & & DTsizeH(ifile,isVvel), & & DTsizeV(ifile,isVvel), & & Kh, Kv, & & pm, pn, & & pmon_p, pnom_r, & & vmask, pmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=JstrP,JendT DO i=IstrT,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO ! my_dot=0.0_r8 DO k=1,N(ng) DO j=JstrP,JendT DO i=IstrT,IendT my_dot=my_dot+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (V3_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (V3_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN VnormV(ic,jc,kc,ifile)=cff END IF END DO END DO END DO CALL dabc_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormV(:,:,:,ifile)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormV(:,:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idVvel, NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idVvel), & & NRM(ifile,ng)%Rindex, & & vmask, & & VnormV(:,:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 1441, MyFile)) RETURN END IF ! ! 3D norm at RHO-points. ! IF (Master) THEN Lsame=.FALSE. DO itrc=1,NT(ng) is=isTvar(itrc) IF (Cnorm(ifile,is)) Lsame=.TRUE. END DO IF (Lsame) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at RHO-points' CALL my_flush (stdout) END IF END IF ! ! Check if the decorrelation scales for all the tracers are different. ! If not, just compute the normalization factors for the first tracer ! and assign the same value to the rest. Recall that this computation ! is very expensive. ! Ldiffer=.FALSE. Imin=1 Imax=Lm(ng) Jmin=1 Jmax=Mm(ng) DO itrc=2,NT(ng) IF ((Hdecay(ifile,isTvar(itrc ),ng).ne. & & Hdecay(ifile,isTvar(itrc-1),ng)).or. & & (Vdecay(ifile,isTvar(itrc ),ng).ne. & & Vdecay(ifile,isTvar(itrc-1),ng))) THEN Ldiffer=.TRUE. END IF END DO IF (.not.Ldiffer) THEN Lsame=.TRUE. UBt=1 ELSE Lsame=.FALSE. UBt=NT(ng) END IF ! DO j=JstrT,JendT DO i=IstrT,IendT cff=om_r(i,j)*on_r(i,j) DO k=1,N(ng) Vscale(i,j,k)=1.0_r8/SQRT(cff*Hz(i,j,k)) END DO END DO END DO DO itrc=1,UBt is=isTvar(itrc) IF (Cnorm(ifile,is)) THEN DO kc=1,N(ng) DO jc=Jmin,Jmax DO ic=Imin,Imax compute=0.0_r8 IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (rmask(ic,jc).gt.0) compute=1.0_r8 END IF CALL mp_reduce (ng, iTLM, 1, compute, 'SUM') IF (compute.gt.0.0_r8) THEN DO k=1,N(ng) DO j=LBj,UBj DO i=LBi,UBi A3d(i,j,k)=0.0_r8 END DO END DO END DO IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN A3d(ic,jc,kc)=1.0_r8 END IF CALL ad_conv_r3d_tile (ng, tile, iADM, & & LBi, UBi, LBj, UBj, & & 1, N(ng), & & IminS, ImaxS, JminS, JmaxS,& & NghostPoints, & & NHsteps(ifile,is)/ifac, & & NVsteps(ifile,is)/ifac, & & DTsizeH(ifile,is), & & DTsizeV(ifile,is), & & Kh, Kv, & & pm, pn, & & pmon_u, pnom_v, & & rmask, umask, vmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO ! my_dot=0.0_r8 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT my_dot=my_dot+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO ! ! Perform parallel global reduction operation: dot product. ! NSUB=1 ! distributed-memory !$OMP CRITICAL (R3_DOT) IF (tile_count.eq.0) THEN Gdotp=my_dot ELSE Gdotp=Gdotp+my_dot END IF tile_count=tile_count+1 IF (tile_count.eq.NSUB) THEN tile_count=0 op_handle='SUM' CALL mp_reduce (ng, iTLM, 1, Gdotp, op_handle) cff=1.0_r8/SQRT(Gdotp) END IF !$OMP END CRITICAL (R3_DOT) ELSE cff=0.0_r8 END IF IF (((Jstr.le.jc).and.(jc.le.Jend)).and. & & ((Istr.le.ic).and.(ic.le.Iend))) THEN IF (Lsame) THEN DO ntrc=1,NT(ng) VnormR(ic,jc,kc,ifile,ntrc)=cff END DO ELSE VnormR(ic,jc,kc,ifile,itrc)=cff END IF END IF END DO END DO END DO END IF END DO DO itrc=1,NT(ng) is=isTvar(itrc) IF (Cnorm(ifile,is)) THEN CALL dabc_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormR(:,:,:,ifile,itrc)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormR(:,:,:,ifile,itrc)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idTvar(itrc), & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idTvar(itrc)), & & NRM(ifile,ng)%Rindex, & & rmask, & & VnormR(:,:,:,ifile,itrc)) END SELECT IF (FoundError(exit_flag, NoError, & & 1656, MyFile)) RETURN END IF END DO END IF END DO FILE_LOOP ! IF (Master) THEN WRITE (stdout,30) END IF 10 FORMAT (/,' Error Covariance Normalization Factors: ', & & 'Exact Method',/) 20 FORMAT (4x,'Computing',1x,a,1x,a) 30 FORMAT (/) ! RETURN END SUBROUTINE normalization_tile ! !*********************************************************************** SUBROUTINE randomization_tile (ng, tile, & & LBi, UBi, LBj, UBj, LBij, UBij, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, ifac, & & pm, om_p, om_r, om_u, om_v, & & pn, on_p, on_r, on_u, on_v, & & pmon_p, pmon_r, pmon_u, & & pnom_p, pnom_r, pnom_v, & & pmask, rmask, & & umask, vmask, & & h, & & Hz, z_r, z_w, & & Kh, & & Kv, & & VnormR, VnormU, VnormV, & & HnormR, HnormU, HnormV) !*********************************************************************** ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew, ifac ! real(r8), intent(in) :: pm(LBi:,LBj:) real(r8), intent(in) :: om_p(LBi:,LBj:) real(r8), intent(in) :: om_r(LBi:,LBj:) real(r8), intent(in) :: om_u(LBi:,LBj:) real(r8), intent(in) :: om_v(LBi:,LBj:) real(r8), intent(in) :: pn(LBi:,LBj:) real(r8), intent(in) :: on_p(LBi:,LBj:) real(r8), intent(in) :: on_r(LBi:,LBj:) real(r8), intent(in) :: on_u(LBi:,LBj:) real(r8), intent(in) :: on_v(LBi:,LBj:) real(r8), intent(in) :: pmon_p(LBi:,LBj:) real(r8), intent(in) :: pmon_r(LBi:,LBj:) real(r8), intent(in) :: pmon_u(LBi:,LBj:) real(r8), intent(in) :: pnom_p(LBi:,LBj:) real(r8), intent(in) :: pnom_r(LBi:,LBj:) real(r8), intent(in) :: pnom_v(LBi:,LBj:) real(r8), intent(in) :: pmask(LBi:,LBj:) real(r8), intent(in) :: rmask(LBi:,LBj:) real(r8), intent(in) :: umask(LBi:,LBj:) real(r8), intent(in) :: vmask(LBi:,LBj:) real(r8), intent(in) :: Kh(LBi:,LBj:) real(r8), intent(in) :: Kv(LBi:,LBj:,0:) real(r8), intent(inout) :: h(LBi:,LBj:) real(r8), intent(out) :: VnormR(LBi:,LBj:,:,:,:) real(r8), intent(out) :: VnormU(LBi:,LBj:,:,:) real(r8), intent(out) :: VnormV(LBi:,LBj:,:,:) real(r8), intent(out) :: HnormR(LBi:,LBj:,:) real(r8), intent(out) :: HnormU(LBi:,LBj:,:) real(r8), intent(out) :: HnormV(LBi:,LBj:,:) real(r8), intent(out) :: Hz(LBi:,LBj:,:) real(r8), intent(out) :: z_r(LBi:,LBj:,:) real(r8), intent(out) :: z_w(LBi:,LBj:,0:) ! ! Local variable declarations. ! logical :: Ldiffer, Lsame ! integer :: i, ifile, is, iter, j, rec integer :: UBt, itrc, k integer :: start(4), total(4) ! real(dp) :: my_time real(r8) :: Aavg, Amax, Amin, Asqr, FacAvg, FacSqr real(r8) :: cff, val ! real(r8), dimension(LBi:UBi,LBj:UBj) :: A2d real(r8), dimension(LBi:UBi,LBj:UBj) :: A2davg real(r8), dimension(LBi:UBi,LBj:UBj) :: A2dsqr real(r8), dimension(LBi:UBi,LBj:UBj) :: Hscale real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: A3d real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: A3davg real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: A3dsqr real(r8), dimension(LBi:UBi,LBj:UBj,1:N(ng)) :: Vscale ! character (len=40 ) :: Text character (len=256) :: ncname character (len=*), parameter :: MyFile = & & "ROMS/Utility/normalization.F"//", randomization_tile" ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU integer :: Iend, IendB, IendP, IendR, IendT integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV integer :: Jend, JendB, JendP, JendR, JendT integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1 integer :: Iendp1, Iendp2, Iendp2i, Iendp3 integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1 integer :: Jendp1, Jendp2, Jendp2i, Jendp3 ! Istr =BOUNDS(ng) % Istr (tile) IstrB =BOUNDS(ng) % IstrB (tile) IstrM =BOUNDS(ng) % IstrM (tile) IstrP =BOUNDS(ng) % IstrP (tile) IstrR =BOUNDS(ng) % IstrR (tile) IstrT =BOUNDS(ng) % IstrT (tile) IstrU =BOUNDS(ng) % IstrU (tile) Iend =BOUNDS(ng) % Iend (tile) IendB =BOUNDS(ng) % IendB (tile) IendP =BOUNDS(ng) % IendP (tile) IendR =BOUNDS(ng) % IendR (tile) IendT =BOUNDS(ng) % IendT (tile) Jstr =BOUNDS(ng) % Jstr (tile) JstrB =BOUNDS(ng) % JstrB (tile) JstrM =BOUNDS(ng) % JstrM (tile) JstrP =BOUNDS(ng) % JstrP (tile) JstrR =BOUNDS(ng) % JstrR (tile) JstrT =BOUNDS(ng) % JstrT (tile) JstrV =BOUNDS(ng) % JstrV (tile) Jend =BOUNDS(ng) % Jend (tile) JendB =BOUNDS(ng) % JendB (tile) JendP =BOUNDS(ng) % JendP (tile) JendR =BOUNDS(ng) % JendR (tile) JendT =BOUNDS(ng) % JendT (tile) ! Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3 Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2 Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1 IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2 IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1 Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1 Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2 Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3 Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3 Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2 Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1 JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2 JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1 Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1 Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2 Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3 ! SourceFile=MyFile my_time=tdays(ng)*day2sec ! !----------------------------------------------------------------------- ! Compute time invariant depths (use zero free-surface). !----------------------------------------------------------------------- ! DO i=LBi,UBi DO j=LBj,UBj A2d(i,j)=0.0_r8 END DO END DO CALL set_depth_tile (ng, tile, iNLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & h, & & A2d, & & Hz, z_r, z_w) ! !----------------------------------------------------------------------- ! Compute initial conditions and model erro covariance, B, ! normalization factors using the randomization approach of Fisher ! and Courtier (1995). These factors ensure that the diagonal ! elements of B are equal to unity. Notice that in applications ! with land/sea masking, the boundary conditions will produce ! large changes in the covariance structures near the boundary. ! ! Initialize factors with randon numbers ("white-noise") having an ! uniform distribution (zero mean and unity variance). Then, scale ! by the inverse squared root area (2D) or volume (3D) and "color" ! with the diffusion operator. Iterate this step over a specified ! number of ensamble members, Nrandom. !----------------------------------------------------------------------- ! IF (Master) WRITE (stdout,10) FILE_LOOP : DO ifile=1,NSA IF (LwrtNRM(ifile,ng)) THEN IF (ifile.eq.1) THEN Text='initial conditions' ELSE IF (ifile.eq.2) THEN Text='model' END IF ! ! Set randomization summation factors. ! FacAvg=1.0_r8/REAL(Nrandom,r8) FacSqr=SQRT(REAL(Nrandom,r8)) ! ! Set time record index to write in normalization NetCDF file. ! ncname=NRM(ifile,ng)%name NRM(ifile,ng)%Rindex=NRM(ifile,ng)%Rindex+1 NRM(ifile,ng)%Nrec=NRM(ifile,ng)%Nrec+1 ! ! Write out model time (s). ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL netcdf_put_fvar (ng, iTLM, ncname, & & Vname(1,idtime), my_time, & & start = (/NRM(ifile,ng)%Rindex/), & & total = (/1/), & & ncid = NRM(ifile,ng)%ncid, & & varid = NRM(ifile,ng)%Vid(idtime)) END SELECT IF (FoundError(exit_flag, NoError, 3709, MyFile)) RETURN ! ! 2D norm at RHO-points. ! IF (Cnorm(ifile,isFsur)) THEN IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at RHO-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrT,IendT A2davg(i,j)=0.0_r8 A2dsqr(i,j)=0.0_r8 Hscale(i,j)=1.0_r8/SQRT(om_r(i,j)*on_r(i,j)) END DO END DO DO iter=1,Nrandom CALL white_noise2d (ng, iTLM, r2dvar, Rscheme(ng), & & IstrR, IendR, JstrR, JendR, & & LBi, UBi, LBj, UBj, & & Amin, Amax, A2d) DO j=JstrT,JendT DO i=IstrT,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO CALL tl_conv_r2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isFsur)/ifac, & & DTsizeH(ifile,isFsur), & & Kh, & & pm, pn, pmon_u, pnom_v, & & rmask, umask, vmask, & & A2d) DO j=Jstr,Jend DO i=Istr,Iend A2davg(i,j)=A2davg(i,j)+A2d(i,j) A2dsqr(i,j)=A2dsqr(i,j)+A2d(i,j)*A2d(i,j) END DO END DO END DO DO j=Jstr,Jend DO i=Istr,Iend Aavg=FacAvg*A2davg(i,j) Asqr=FacAvg*A2dsqr(i,j) IF (rmask(i,j).gt.0.0_r8) THEN HnormR(i,j,ifile)=1.0_r8/SQRT(Asqr) ELSE HnormR(i,j,ifile)=0.0_r8 END IF END DO END DO CALL dabc_r2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormR(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormR(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idFsur, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idFsur), & & NRM(ifile,ng)%Rindex, & & rmask, & & HnormR(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 3813, MyFile)) RETURN END IF ! ! 2D norm at U-points. ! IF (Cnorm(ifile,isUbar)) THEN IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at U-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrP,IendT A2davg(i,j)=0.0_r8 A2dsqr(i,j)=0.0_r8 Hscale(i,j)=1.0_r8/SQRT(om_u(i,j)*on_u(i,j)) END DO END DO DO iter=1,Nrandom CALL white_noise2d (ng, iTLM, u2dvar, Rscheme(ng), & & Istr, IendR, JstrR, JendR, & & LBi, UBi, LBj, UBj, & & Amin, Amax, A2d) DO j=JstrT,JendT DO i=IstrP,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO CALL tl_conv_u2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isUbar)/ifac, & & DTsizeH(ifile,isUbar), & & Kh, & & pm, pn, pmon_r, pnom_p, & & umask, pmask, & & A2d) DO j=Jstr,Jend DO i=IstrU,Iend A2davg(i,j)=A2davg(i,j)+A2d(i,j) A2dsqr(i,j)=A2dsqr(i,j)+A2d(i,j)*A2d(i,j) END DO END DO END DO DO j=Jstr,Jend DO i=IstrU,Iend Aavg=FacAvg*A2davg(i,j) Asqr=FacAvg*A2dsqr(i,j) IF (umask(i,j).gt.0.0_r8) THEN HnormU(i,j,ifile)=1.0_r8/SQRT(Asqr) ELSE HnormU(i,j,ifile)=0.0_r8 END IF END DO END DO CALL dabc_u2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormU(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormU(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idUbar, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idUbar), & & NRM(ifile,ng)%Rindex, & & umask, & & HnormU(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 3918, MyFile)) RETURN END IF ! ! 2D norm at V-points. ! IF (Cnorm(ifile,isVbar)) THEN IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '2D normalization factors at V-points' CALL my_flush (stdout) END IF DO j=JstrP,JendT DO i=IstrT,IendT A2davg(i,j)=0.0_r8 A2dsqr(i,j)=0.0_r8 Hscale(i,j)=1.0_r8/SQRT(om_v(i,j)*on_v(i,j)) END DO END DO DO iter=1,Nrandom CALL white_noise2d (ng, iTLM, v2dvar, Rscheme(ng), & & IstrR, IendR, Jstr, JendR, & & LBi, UBi, LBj, UBj, & & Amin, Amax, A2d) DO j=JstrP,JendT DO i=IstrT,IendT A2d(i,j)=A2d(i,j)*Hscale(i,j) END DO END DO CALL tl_conv_v2d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isVbar)/ifac, & & DTsizeH(ifile,isVbar), & & Kh, & & pm, pn, pmon_p, pnom_r, & & vmask, pmask, & & A2d) DO j=JstrV,Jend DO i=Istr,Iend A2davg(i,j)=A2davg(i,j)+A2d(i,j) A2dsqr(i,j)=A2dsqr(i,j)+A2d(i,j)*A2d(i,j) END DO END DO END DO DO j=JstrV,Jend DO i=Istr,Iend Aavg=FacAvg*A2davg(i,j) Asqr=FacAvg*A2dsqr(i,j) IF (vmask(i,j).gt.0.0_r8) THEN HnormV(i,j,ifile)=1.0_r8/SQRT(Asqr) ELSE HnormV(i,j,ifile)=0.0_r8 END IF END DO END DO CALL dabc_v2d_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & HnormV(:,:,ifile)) CALL mp_exchange2d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & HnormV(:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm2d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, idVbar, & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idVbar), & & NRM(ifile,ng)%Rindex, & & vmask, & & HnormV(:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 4023, MyFile)) RETURN END IF ! ! 3D norm U-points. ! IF (Cnorm(ifile,isUvel)) THEN IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at U-points' CALL my_flush (stdout) END IF DO j=JstrT,JendT DO i=IstrP,IendT val=om_u(i,j)*on_u(i,j)*0.5_r8 DO k=1,N(ng) A3davg(i,j,k)=0.0_r8 A3dsqr(i,j,k)=0.0_r8 Vscale(i,j,k)=1.0_r8/SQRT(val*(Hz(i-1,j,k)+Hz(i,j,k))) END DO END DO END DO DO iter=1,Nrandom CALL white_noise3d (ng, iTLM, u3dvar, Rscheme(ng), & & Istr, IendR, JstrR, JendR, & & LBi, UBi, LBj, UBj, 1, N(ng), & & Amin, Amax, A3d) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO CALL tl_conv_u3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isUvel)/ifac, & & NVsteps(ifile,isUvel)/ifac, & & DTsizeH(ifile,isUvel), & & DTsizeV(ifile,isUvel), & & Kh, Kv, & & pm, pn, & & pmon_r, pnom_p, & & umask, pmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=Jstr,Jend DO i=IstrU,Iend A3davg(i,j,k)=A3davg(i,j,k)+A3d(i,j,k) A3dsqr(i,j,k)=A3dsqr(i,j,k)+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO END DO DO k=1,N(ng) DO j=Jstr,Jend DO i=IstrU,Iend Aavg=FacAvg*A3davg(i,j,k) Asqr=FacAvg*A3dsqr(i,j,k) IF (umask(i,j).gt.0.0_r8) THEN VnormU(i,j,k,ifile)=1.0_r8/SQRT(Asqr) ELSE VnormU(i,j,k,ifile)=0.0_r8 END IF END DO END DO END DO CALL dabc_u3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormU(:,:,:,ifile)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormU(:,:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idUvel, NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idUvel), & & NRM(ifile,ng)%Rindex, & & umask, & & VnormU(:,:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 4151, MyFile)) RETURN END IF ! ! 3D norm at V-points. ! IF (Cnorm(ifile,isVvel)) THEN IF (Master) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at V-points' CALL my_flush (stdout) END IF DO j=JstrP,JendT DO i=IstrT,IendT val=om_v(i,j)*on_v(i,j)*0.5_r8 DO k=1,N(ng) A3davg(i,j,k)=0.0_r8 A3dsqr(i,j,k)=0.0_r8 Vscale(i,j,k)=1.0_r8/SQRT(val*(Hz(i,j-1,k)+Hz(i,j,k))) END DO END DO END DO DO iter=1,Nrandom CALL white_noise3d (ng, iTLM, v3dvar, Rscheme(ng), & & IstrR, IendR, Jstr, JendR, & & LBi, UBi, LBj, UBj, 1, N(ng), & & Amin, Amax, A3d) DO k=1,N(ng) DO j=JstrP,JendT DO i=IstrT,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO CALL tl_conv_v3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,isVvel)/ifac, & & NVsteps(ifile,isVvel)/ifac, & & DTsizeH(ifile,isVvel), & & DTsizeV(ifile,isVvel), & & Kh, Kv, & & pm, pn, & & pmon_p, pnom_r, & & vmask, pmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=JstrV,Jend DO i=Istr,Iend A3davg(i,j,k)=A3davg(i,j,k)+A3d(i,j,k) A3dsqr(i,j,k)=A3dsqr(i,j,k)+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO END DO DO k=1,N(ng) DO j=JstrV,Jend DO i=Istr,Iend Aavg=FacAvg*A3davg(i,j,k) Asqr=FacAvg*A3dsqr(i,j,k) IF (vmask(i,j).gt.0.0_r8) THEN VnormV(i,j,k,ifile)=1.0_r8/SQRT(Asqr) ELSE VnormV(i,j,k,ifile)=0.0_r8 END IF END DO END DO END DO CALL dabc_v3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormV(:,:,:,ifile)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormV(:,:,:,ifile)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idVvel, NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idVvel), & & NRM(ifile,ng)%Rindex, & & vmask, & & VnormV(:,:,:,ifile)) END SELECT IF (FoundError(exit_flag, NoError, 4277, MyFile)) RETURN END IF ! ! 3D norm at RHO-points. ! IF (Master) THEN Lsame=.FALSE. DO itrc=1,NT(ng) is=isTvar(itrc) IF (Cnorm(ifile,is)) Lsame=.TRUE. END DO IF (Lsame) THEN WRITE (stdout,20) TRIM(Text), & & '3D normalization factors at RHO-points' CALL my_flush (stdout) END IF END IF ! ! Check if the decorrelation scales for all the tracers are different. ! If not, just compute the normalization factors for the first tracer ! and assign the same value to the rest. Recall that this computation ! is very expensive. ! Ldiffer=.FALSE. DO itrc=2,NT(ng) IF ((Hdecay(ifile,isTvar(itrc ),ng).ne. & & Hdecay(ifile,isTvar(itrc-1),ng)).or. & & (Vdecay(ifile,isTvar(itrc ),ng).ne. & & Vdecay(ifile,isTvar(itrc-1),ng))) THEN Ldiffer=.TRUE. END IF END DO IF (.not.Ldiffer) THEN Lsame=.TRUE. UBt=1 ELSE Lsame=.FALSE. UBt=NT(ng) END IF ! DO j=JstrT,JendT DO i=IstrT,IendT val=om_r(i,j)*on_r(i,j) DO k=1,N(ng) Vscale(i,j,k)=1.0_r8/SQRT(val*Hz(i,j,k)) END DO END DO END DO DO itrc=1,UBt is=isTvar(itrc) IF (Cnorm(ifile,is)) THEN DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT A3davg(i,j,k)=0.0_r8 A3dsqr(i,j,k)=0.0_r8 END DO END DO END DO DO iter=1,Nrandom CALL white_noise3d (ng, iTLM, r3dvar, Rscheme(ng), & & IstrR, IendR, JstrR, JendR, & & LBi, UBi, LBj, UBj, 1, N(ng), & & Amin, Amax, A3d) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT A3d(i,j,k)=A3d(i,j,k)*Vscale(i,j,k) END DO END DO END DO CALL tl_conv_r3d_tile (ng, tile, iTLM, & & LBi, UBi, LBj, UBj, 1, N(ng), & & IminS, ImaxS, JminS, JmaxS, & & NghostPoints, & & NHsteps(ifile,is)/ifac, & & NVsteps(ifile,is)/ifac, & & DTsizeH(ifile,is), & & DTsizeV(ifile,is), & & Kh, Kv, & & pm, pn, & & pmon_u, pnom_v, & & rmask, umask, vmask, & & Hz, z_r, & & A3d) DO k=1,N(ng) DO j=Jstr,Jend DO i=Istr,Iend A3davg(i,j,k)=A3davg(i,j,k)+A3d(i,j,k) A3dsqr(i,j,k)=A3dsqr(i,j,k)+A3d(i,j,k)*A3d(i,j,k) END DO END DO END DO END DO DO k=1,N(ng) DO j=Jstr,Jend DO i=Istr,Iend Aavg=FacAvg*A3davg(i,j,k) Asqr=FacAvg*A3dsqr(i,j,k) IF (rmask(i,j).gt.0.0_r8) THEN VnormR(i,j,k,ifile,itrc)=1.0_r8/SQRT(Asqr) ELSE VnormR(i,j,k,ifile,itrc)=0.0_r8 END IF END DO END DO END DO END IF END DO IF (Lsame) THEN DO itrc=2,NT(ng) DO k=1,N(ng) DO j=Jstr,Jend DO i=Istr,Iend VnormR(i,j,k,ifile,itrc)=VnormR(i,j,k,ifile,1) END DO END DO END DO END DO END IF DO itrc=1,NT(ng) is=isTvar(itrc) IF (Cnorm(ifile,is)) THEN CALL dabc_r3d_tile (ng, tile, & & LBi, UBi, LBj, UBj, 1, N(ng), & & VnormR(:,:,:,ifile,itrc)) CALL mp_exchange3d (ng, tile, iTLM, 1, & & LBi, UBi, LBj, UBj, 1, N(ng), & & NghostPoints, & & EWperiodic(ng), NSperiodic(ng), & & VnormR(:,:,:,ifile,itrc)) ! SELECT CASE (NRM(ifile,ng)%IOtype) CASE (io_nf90) CALL wrt_norm3d_nf90 (ng, tile, iTLM, ncname, & & LBi, UBi, LBj, UBj, 1, N(ng), & & idTvar(itrc), & & NRM(ifile,ng)%ncid, & & NRM(ifile,ng)%Vid(idTvar(itrc)), & & NRM(ifile,ng)%Rindex, & & rmask, & & VnormR(:,:,:,ifile,itrc)) END SELECT IF (FoundError(exit_flag, NoError, & & 4456, MyFile)) RETURN END IF END DO END IF END DO FILE_LOOP ! IF (Master) THEN WRITE (stdout,30) END IF 10 FORMAT (/,' Error Covariance Factors: Randomization Method',/) 20 FORMAT (4x,'Computing',1x,a,1x,a) 30 FORMAT (/) ! RETURN END SUBROUTINE randomization_tile ! !*********************************************************************** SUBROUTINE wrt_norm2d_nf90 (ng, tile, model, ncname, & & LBi, UBi, LBj, UBj, & & ifield, ncid, ncvarid, tindex, & & Amask, & & A) !*********************************************************************** ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: ifield, ncid, ncvarid, tindex character (len=*), intent(in) :: ncname ! real(r8), intent(in) :: Amask(LBi:,LBj:) real(r8), intent(in) :: A(LBi:,LBj:) ! ! Local variable declarations. ! integer :: gfactor, gtype, status ! real(dp) :: scale ! character (len=*), parameter :: MyFile = & & "ROMS/Utility/normalization.F"//", wrt_norm2d_nf90" ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU integer :: Iend, IendB, IendP, IendR, IendT integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV integer :: Jend, JendB, JendP, JendR, JendT integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1 integer :: Iendp1, Iendp2, Iendp2i, Iendp3 integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1 integer :: Jendp1, Jendp2, Jendp2i, Jendp3 ! Istr =BOUNDS(ng) % Istr (tile) IstrB =BOUNDS(ng) % IstrB (tile) IstrM =BOUNDS(ng) % IstrM (tile) IstrP =BOUNDS(ng) % IstrP (tile) IstrR =BOUNDS(ng) % IstrR (tile) IstrT =BOUNDS(ng) % IstrT (tile) IstrU =BOUNDS(ng) % IstrU (tile) Iend =BOUNDS(ng) % Iend (tile) IendB =BOUNDS(ng) % IendB (tile) IendP =BOUNDS(ng) % IendP (tile) IendR =BOUNDS(ng) % IendR (tile) IendT =BOUNDS(ng) % IendT (tile) Jstr =BOUNDS(ng) % Jstr (tile) JstrB =BOUNDS(ng) % JstrB (tile) JstrM =BOUNDS(ng) % JstrM (tile) JstrP =BOUNDS(ng) % JstrP (tile) JstrR =BOUNDS(ng) % JstrR (tile) JstrT =BOUNDS(ng) % JstrT (tile) JstrV =BOUNDS(ng) % JstrV (tile) Jend =BOUNDS(ng) % Jend (tile) JendB =BOUNDS(ng) % JendB (tile) JendP =BOUNDS(ng) % JendP (tile) JendR =BOUNDS(ng) % JendR (tile) JendT =BOUNDS(ng) % JendT (tile) ! Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3 Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2 Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1 IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2 IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1 Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1 Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2 Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3 Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3 Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2 Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1 JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2 JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1 Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1 Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2 Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3 ! !----------------------------------------------------------------------- ! Write out requested 2D field into normalization NetCDF file. Since ! the computation of normalization coefficients is a very expensive ! computation, synchronize file to disk. !----------------------------------------------------------------------- ! IF (FoundError(exit_flag, NoError, 6004, MyFile)) RETURN ! ! Set grid type factor to write full (gfactor=1) fields or water ! points (gfactor=-1) fields only. ! gfactor=1 ! ! Write out 2D normalization field. ! gtype=gfactor*Iinfo(1,ifield,ng) scale=1.0_dp status=nf_fwrite2d(ng, model, ncid, ifield, & & ncvarid, tindex, gtype, & & LBi, UBi, LBj, UBj, scale, & & Amask, & & A) IF (FoundError(status, nf90_noerr, 6026, MyFile)) THEN IF (Master) THEN WRITE (stdout,10) TRIM(Vname(1,ifield)), tindex END IF exit_flag=3 ioerror=status RETURN END IF ! ! Synchronize normalization NetCDF file to disk to allow other ! processes to access data immediately after it is written. ! CALL netcdf_sync (ng, model, ncname, ncid) IF (FoundError(exit_flag, NoError, 6039, MyFile)) RETURN IF (Master) WRITE (stdout,20) TRIM(Vname(1,ifield)), tindex CALL my_flush (stdout) ! 10 FORMAT (/,' WRT_NORM2D_NF90 - error while writing variable: ',a, & & /,19x 'into normalization NetCDF file for time record: ', & & i0) 20 FORMAT (7x,'wrote ',a, t21,'normalization factors into record ', & & i0) ! RETURN END SUBROUTINE wrt_norm2d_nf90 ! !*********************************************************************** SUBROUTINE wrt_norm3d_nf90 (ng, tile, model, ncname, & & LBi, UBi, LBj, UBj, LBk, UBk, & & ifield, ncid, ncvarid, tindex, & & Amask, & & A) !*********************************************************************** ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk integer, intent(in) :: ifield, ncid, ncvarid, tindex character (len=*), intent(in) :: ncname ! real(r8), intent(in) :: Amask(LBi:,LBj:) real(r8), intent(in) :: A(LBi:,LBj:,LBk:) ! ! Local variable declarations. ! integer :: gfactor, gtype, status ! real(dp) :: scale ! character (len=*), parameter :: MyFile = & & "ROMS/Utility/normalization.F"//", wrt_norm3d_nf90" ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU integer :: Iend, IendB, IendP, IendR, IendT integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV integer :: Jend, JendB, JendP, JendR, JendT integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1 integer :: Iendp1, Iendp2, Iendp2i, Iendp3 integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1 integer :: Jendp1, Jendp2, Jendp2i, Jendp3 ! Istr =BOUNDS(ng) % Istr (tile) IstrB =BOUNDS(ng) % IstrB (tile) IstrM =BOUNDS(ng) % IstrM (tile) IstrP =BOUNDS(ng) % IstrP (tile) IstrR =BOUNDS(ng) % IstrR (tile) IstrT =BOUNDS(ng) % IstrT (tile) IstrU =BOUNDS(ng) % IstrU (tile) Iend =BOUNDS(ng) % Iend (tile) IendB =BOUNDS(ng) % IendB (tile) IendP =BOUNDS(ng) % IendP (tile) IendR =BOUNDS(ng) % IendR (tile) IendT =BOUNDS(ng) % IendT (tile) Jstr =BOUNDS(ng) % Jstr (tile) JstrB =BOUNDS(ng) % JstrB (tile) JstrM =BOUNDS(ng) % JstrM (tile) JstrP =BOUNDS(ng) % JstrP (tile) JstrR =BOUNDS(ng) % JstrR (tile) JstrT =BOUNDS(ng) % JstrT (tile) JstrV =BOUNDS(ng) % JstrV (tile) Jend =BOUNDS(ng) % Jend (tile) JendB =BOUNDS(ng) % JendB (tile) JendP =BOUNDS(ng) % JendP (tile) JendR =BOUNDS(ng) % JendR (tile) JendT =BOUNDS(ng) % JendT (tile) ! Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3 Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2 Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1 IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2 IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1 Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1 Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2 Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3 Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3 Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2 Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1 JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2 JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1 Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1 Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2 Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3 ! !----------------------------------------------------------------------- ! Write out requested 3D field into normalization NetCDF file. Since ! the computation of normalization coefficients is a very expensive ! computation, synchronize file to disk. !----------------------------------------------------------------------- ! IF (FoundError(exit_flag, NoError, 6195, MyFile)) RETURN ! ! Set grid type factor to write full (gfactor=1) fields or water ! points (gfactor=-1) fields only. ! gfactor=1 ! ! Write out 3D normalization field. ! gtype=gfactor*Iinfo(1,ifield,ng) scale=1.0_dp status=nf_fwrite3d(ng, model, ncid, ifield, & & ncvarid, tindex, gtype, & & LBi, UBi, LBj, UBj, LBk, UBk, scale, & & Amask, & & A) IF (FoundError(status, nf90_noerr, 6217, MyFile)) THEN IF (Master) THEN WRITE (stdout,10) TRIM(Vname(1,ifield)), tindex END IF exit_flag=3 ioerror=status RETURN END IF ! ! Synchronize normalization NetCDF file to disk to allow other ! processes to access data immediately after it is written. ! CALL netcdf_sync (ng, model, ncname, ncid) IF (FoundError(exit_flag, NoError, 6230, MyFile)) RETURN IF (Master) WRITE (stdout,20) TRIM(Vname(1,ifield)), tindex CALL my_flush (stdout) ! 10 FORMAT (/,' WRT_NORM3D_NF90 - error while writing variable: ',a, & & /,19x,'into normalization NetCDF file for time record: ', & & i0) 20 FORMAT (7x,'wrote ',a, t21,'normalization factors into record ', & & i0) ! RETURN END SUBROUTINE wrt_norm3d_nf90 END MODULE normalization_mod