UKMO Staging [fb_rtd2, bf_coupling condition, bf_smc_omp] (#254)

* Fb rtd2 (#7)

FB_rtd2: Output boundary conditions to rotated grids

* This feature allows the user to specify the pole of each nested grid output of boundary conditions (nest) file.
  - There is no effect unless the switch RTD is set for compilation.
  - The contents of b.c. output files nestI.ww3 are not affected.
  - A model compiled with combined switches SMC and RTD is not affected.

* Updated manual Sec. "3.4.9 Rotated grids"

* Added regtests for testing input/output BCs in rotated pole context.

* Updated revision of ww3_bound.ftn to 7.11

* Added new ww3_tr1 regtest (rotated pole) to matrix.base. Also, split the SMC/RTD regtests into their own separate switches.

Co-authored-by: ukmo-chris.bunney <christopher.bunney@metoffice.gov.uk>

* Double allocation error of FIELD variable [w3wavemd] (#14)

* Fixed allocation of FIELD variable.
* Removed defunct OMPX switches and FLOMP variable

* Change coupling condition so that it works with all compilers (#12)

Co-authored-by: Carsten Hansen <cha@fcoo.dk>
Co-authored-by: Juan Manuel Castillo Sanchez <48921434+ukmo-juan-castillo@users.noreply.github.com>

manual/num/rotagrid.tex

@@ -4,8 +4,8 @@ \subsubsection{~Rotated grids} \label{sub:num_space_rotagrid}
 
 \noindent
 The rotated grid is a latitude-longitude (lat-lon) grid and is obtained by
-rotating the North Pole to a new position at latitude $\phi_{p}$ and
-longitude $\lambda_{p}$ in the standard latitude-longitude system.  The new
+rotating the North Pole along a longitude $\lambda_{p}$ to a new position at
+latitude $\phi_{p}$ in the standard latitude-longitude system.  The new
 pole position is chosen so that the model domain of interest may be placed
 around the rotated equatorial area for a evenly-spaced lat-lon mesh. For this
 reason the rotated grid is also known as \emph{Equatorial grid}. For instance,
@@ -15,22 +15,41 @@ \subsubsection{~Rotated grids} \label{sub:num_space_rotagrid}
 grid allows a much more evenly spaced lat-lon mesh in the NAEW domain than the
 standard lat-lon grid in the same area. 
 
-In \ws\, the rotated grid is implemented with minimum changes to the original
+In \ws, the rotated grid is implemented with minimum changes to the original
 lat-lon grid. In fact, the rotated grid is treated just like the standard 
 lat-lon grid inside the model. To set up and run a rotated grid model configuration,
 users should choose the regular lat-lon grid along with the {\code RTD} switch.
-The rotated pole position is set using the {\code PLAT} and {\code PLON} variables in the
-{\bf ww3\_grid.inp} namelist {\code ROTD}. Model input files, like wind, current and ice files 
-should be mapped on to the rotated grid. For convenience of nesting in standard 
-lat-lon grid frameworks, boundary conditions provided to and output from the 
+The rotated pole position is set using the {\code PLAT} and {\code PLON} variables
+under the namelist {\code ROTD} in the input file {\bf ww3\_grid.inp}
+(see Sect.~\ref{sec:config011}). If the pole is set as {\code PLAT = 90.0},
+{\code PLON = -180.0}, the grid is treated as a standard lat-lon system.
+
+Model input files, like wind, current and ice files
+should be mapped on to the rotated grid. For convenience of nesting in standard
+lat-lon grid frameworks, boundary conditions data provided to and output from the
 rotated grid use spectra referenced to a standard grid north and standard lat-lon
-grid points values, which are converted into rotated grid lat-lon inside \ws\.
-The list of 2D spectral output locations in {\bf ww3\_shel.inp} are also specified in 
-standard lat-lon. 
+grid points values, which are converted into rotated grid lat-lon inside \ws.
+The list of 2D spectral output locations in {\bf ww3\_shel.inp} or
+{\bf ww3\_shel.nml} are also specified in standard lat-lon.
+
+When nesting from a standard grid to a rotated grid model, both the outer and inner
+grids should be built with the {\code RTD} switch set when compiling executables for
+the models. When nesting from a rotated grid to a standard grid model, the inner
+(standard) grid model does not require to be built with {\code RTD}.
+
+Output of spectra at boundary points to one-way nested inner grids are transferred
+as described in Appendix~\ref{app:nest}. Output b.c.\ are defined in the input
+file (see {\bf ww3\_grid.inp}, Sect.~\ref{sec:config011}) as a sequence of straight
+lines given in coordinates of the inner grid. If the inner grid is rotated, it's
+pole position must be defined as the values of the array elements
+{\code BPLAT(\sl{n})} and {\code BPLON(\sl{n})} under namelist {\code ROTB}
+for the outer grid. The array index {\sl{n}} is the index ({\sl{1:9}}) of the boundary
+conditions file {\file nest{\sl{n}}.ww3}.
 
-Model directional and x-y vector outputs can be converted to a standard grid
-north reference by setting the UNROT variable in the {\bf ww3\_grid.inp} namelist
-ROTD to True. With this set, for point outputs lat-lon locations all directional
+Model directional and x-y vector outputs from a rotated grid can be converted
+to a standard grid north reference by setting the UNROT variable in the
+{\bf ww3\_grid.inp} namelist ROTD to True. With this set, for point outputs
+lat-lon locations all directional
 values such as wind direction, current direction and 2D spectra are converted
 into standard lat-lon orientation. Functions to de-rotate gridded
 fields are applied in {\bf ww3\_ounf}, {\bf ww3\_outf} and {\bf ww3\_grib}. 

model/ftn/w3iobcmd.ftn

@@ -103,6 +103,8 @@
 !/                  (T.J. Campbell, NRL)                ( version 6.02 )
 !/    01-Mar-2018 : Rotate boundary points and directions 
 !/                  of input spectra for rotated grids  ( version 6.02 )
+!/    07-Oct-2019 : RTD option with standard lat-lon
+!/                  grid when nesting to rotated grid   ( version 7.11 )
 !/
 !  1. Purpose :
 !
@@ -230,8 +232,6 @@
 !/RTD  !!    Declare rotation angle and rotated lat/lon variables for
 !/RTD  !!    boundary points.   JGLi12Jun2012
 !/RTD      REAL, ALLOCATABLE      :: Anglbdy(:), ELatbdy(:), ELonbdy(:) 
-      REAL, ALLOCATABLE      :: XBPO2(:), YBPO2(:)
-!/RTD      REAL, ALLOCATABLE      :: ANG(:)
 !/RTD      REAL                   :: Spectr(NK*NTH)
 !/RTD      REAL                   :: XRLIM, YRLIM
       REAL, ALLOCATABLE       :: TMPSPC(:,:)
@@ -311,26 +311,12 @@
 !/T                WRITE (NDST,9002) IFILE, NDSL(IFILE), IDSTRBC,       &
 !/T                       VERBPTBC, NBO(IFILE)-NBO(IFILE-1)
 !
-!
-                 ! Temporary storage for boundary points:
-                 ! The arrays are required to allow rotation of
-                 ! x,y coordinate data under the /RTD option
-                 ALLOCATE( XBPO2(NBO(IFILE)-NBO(IFILE-1)) )
-                 ALLOCATE( YBPO2(NBO(IFILE)-NBO(IFILE-1)) )
-!/RTD                 ALLOCATE( ANG(NBO(IFILE)-NBO(IFILE-1)) )
-                 XBPO2=XBPO(NBO(IFILE-1)+1:NBO(IFILE))
-                 YBPO2=YBPO(NBO(IFILE-1)+1:NBO(IFILE))
-!
-!/RTD                 ! Convert rotated lat/lon back to standard pole grid.
-!/RTD                 CALL W3EQTOLL ( YBPO(NBO(IFILE-1)+1:NBO(IFILE)),   &
-!/RTD                                 XBPO(NBO(IFILE-1)+1:NBO(IFILE)),   &
-!/RTD                                 YBPO2,XBPO2,ANG,POLAT,POLON,       &
-!/RTD                                 NBO(IFILE)-NBO(IFILE-1) )
-!/RTD
-                ! Modified write below to write out XBPO2 and YBPO2,
-                ! instead of XBPO/YBPO
+!/RTD ! By running the ww3_grid program the arrays XBPO, YBPO have been
+!/RTD ! remapped to standard lat-lon and stored in mod_def.*
+!/RTD !
                 WRITE (NDSL(IFILE))                                   &
-                       XBPO2, YBPO2,                                  &
+                       (XBPO(I),I=NBO(IFILE-1)+1,NBO(IFILE)),         &
+                       (YBPO(I),I=NBO(IFILE-1)+1,NBO(IFILE)),         &
                      ((IPBPO(I,J),I=NBO(IFILE-1)+1,NBO(IFILE)),J=1,4),&
                      ((RDBPO(I,J),I=NBO(IFILE-1)+1,NBO(IFILE)),J=1,4)
 !
@@ -340,10 +326,6 @@
 !/T0                         YBPO(I), (IPBPO(I,J),J=1,4),             &
 !/T0                         (RDBPO(I,J),J=1,4)
 !/T0                  END DO
-!
-                ! Deallocate temporary storage for boundary points
-                DEALLOCATE( XBPO2, YBPO2 )
-!/RTD                DEALLOCATE( ANG )
 !
                 END DO
             END IF
@@ -404,6 +386,10 @@
               ((IPBPI(I,J),I=1,NBI),J=1,4),                           &
               ((RDBPI(I,J),I=1,NBI),J=1,4)
 !
+!/RTD ! All boundary conditions position arrays XBPI, YBPI are defined
+!/RTD ! in standard lat/lon coordinates. If Polat = 90. (and Polon = -180.),
+!/RTD ! the b.c. positions don't need to be remapped 
+!/RTD        IF ( Polat < 90. ) THEN
 !/RTD !!   Convert standard into rotated lat/lon. JGLi12Jun2012
 !/RTD          ALLOCATE ( Anglbdy(NBI),  ELatbdy(NBI), ELonbdy(NBI) )
 !/RTD
@@ -433,6 +419,8 @@
 !/RTD          ENDDO
 !/RTD
 !/RTD          DEALLOCATE ( Anglbdy,  ELatbdy, ELonbdy )
+!/RTD
+!/RTD        END IF ! ( Polat < 90. )
 !/RTD
           FLOK   = .TRUE.
           IF (GTYPE .EQ. UNGTYPE) THEN
@@ -541,11 +529,15 @@
 !/DIST                                    CG(1+(IS-1)/NTH,ISEA)
 !/DIST                END DO
 !
-!/RTD              ! Added spectral turning for rotated grid
-!/RTD              ! (rotate back to standard pole)
-!/RTD              Spectr = ABPOS(:,ISOUT)
-!/RTD              CALL W3ACTURN( NTH, NK, -AnglD(ISEA), Spectr )
-!/RTD              ABPOS(:,ISOUT) = Spectr
+!/RTD ! Polat == 90. means the grid is standard lat-lon, and the spectra
+!/RTD ! need not be rotated back 
+!/RTD              IF ( Polat < 90. ) THEN
+!/RTD                ! Added spectral turning for rotated grid
+!/RTD                ! (rotate back to standard pole)
+!/RTD                Spectr = ABPOS(:,ISOUT)
+!/RTD                CALL W3ACTURN( NTH, NK, -AnglD(ISEA), Spectr )
+!/RTD                ABPOS(:,ISOUT) = Spectr
+!/RTD                END IF
 !
               WRITE (NDSL(IFILE)) (ABPOS(IS,ISOUT),IS=1,NSPEC)
 !

model/ftn/w3updtmd.ftn

@@ -54,6 +54,8 @@
 !/                  and add ICEH initialization in W3UICE.
 !/    13-Jan-2016 : Changed initial value of ICEDMAX    ( version 5.08 )
 !/    26-Mar-2018 : Sea-point only Wnd/Cur input.  JGLi ( version 6.04 )
+!/    07-Oct-2019 : RTD option with standard lat-lon
+!/                  grid when nesting to rotated grid   ( version 7.11 )
 !/
 !/    Copyright 2009-2014 National Weather Service (NWS),
 !/       National Oceanic and Atmospheric Administration.  All rights
@@ -903,6 +905,7 @@
 !/    12-Jun-2012 : Add /RTD option or rotated grid option. 
 !/                  (Jian-Guo Li)                       ( version 4.06 )
 !/    06-Jun-2018 : Add DEBUGIOBC/SETUP/DEBUGW3ULEV     ( version 6.04 )
+!/    13-Jun-2019 : Rotation only if POLAT<90 (C.Hansen)( version 7.11 ) 
 !/
 !  1. Purpose :
 !
@@ -955,7 +958,7 @@
 !/ ------------------------------------------------------------------- /
       USE W3GDATMD, ONLY: NSPEC, MAPWN, SIG2, DDEN
 !/RTD !!   Use rotation angle and action conversion sub.  JGLi12Jun2012
-!/RTD      USE W3GDATMD, ONLY: NK, NTH, NSPEC, AnglD 
+!/RTD      USE W3GDATMD, ONLY: NK, NTH, NSPEC, AnglD, PoLat
 !/RTD      USE W3SERVMD, ONLY: W3ACTURN
       USE W3ADATMD, ONLY: CG
       USE W3ODATMD, ONLY: NBI, ABPI0, ABPIN, ISBPI, IPBPI, RDBPI,     &
@@ -1022,10 +1025,13 @@
           END DO
 !
 !/RTD !!  Rotate the spectra if model is on rotated grid.  JGLi12Jun2012
-!/RTD        Spectr = BBPIN(:,IBI)
-!/RTD        AnglBP = AnglD(ISEA)
-!/RTD        CALL  W3ACTURN( NTH, NK, AnglBP, Spectr )
-!/RTD        BBPIN(:,IBI) = Spectr
+!/RTD !!  PoLat == 90. if the grid is standard lat/lon (C. Hansen 20190613)
+!/RTD        IF ( PoLat < 90. ) THEN
+!/RTD          Spectr = BBPIN(:,IBI)
+!/RTD          AnglBP = AnglD(ISEA)
+!/RTD          CALL  W3ACTURN( NTH, NK, AnglBP, Spectr )
+!/RTD          BBPIN(:,IBI) = Spectr
+!/RTD        END IF
 !/RTD
 !
         END DO

model/ftn/w3wavemd.ftn

@@ -80,6 +80,8 @@
 !/                  (M. Accensi & F. Ardhuin, IFREMER)
 !/    27-Aug-2015 : Update for ICEH, ICEF               ( version 5.08 )
 !/    14-Sep-2018 : Remove PALM implementation          ( version 6.06 )
+!/    15-Sep-2020 : Bugfix FIELD allocation. Remove     ( version 7.11 )
+!/                  defunct OMPX switches.
 !/
 !/    Copyright 2009-2014 National Weather Service (NWS),
 !/       National Oceanic and Atmospheric Administration.  All rights
@@ -154,7 +156,6 @@
 !       !/DIST  Id.
 !       !/MPI   Id.
 !       !/OMPG  Id.
-!       !/OMPX  Id.
 !
 !       !/PR1   First order propagation schemes.
 !       !/PR2   ULTIMATE QUICKEST scheme.
@@ -456,8 +457,7 @@
       LOGICAL                 :: FLACT, FLZERO, FLFRST, FLMAP, TSTAMP,&
                                  SKIP_O, FLAG_O, FLDDIR, READBC,      &
                                  FLAG0 = .FALSE., FLOUTG, FLPFLD,     &
-                                 FLPART, LOCAL, FLOUTG2,              &
-                                 FLOMP = .FALSE.
+                                 FLPART, LOCAL, FLOUTG2
 !
 !Li   Logical variable to control regular gird lines in conflict with SMC option.
 !AR   SMC option is in conflict with lofical variables for regular grid ... chicken ... egg ... stuff 
@@ -497,8 +497,6 @@
 !
 ! 0.a Set pointers to data structure
 !
-!/OMPX      FLOMP  = .TRUE.
-!
 !/COU      SCREEN   =  333
 !
 !/DEBUGINIT      WRITE(740+IAPROC,*) 'W3WAVE, step 1'
@@ -576,8 +574,7 @@
           FCUT  = SIG(NK) * TPIINV
         END IF
 !
-!!Li  ALLOCATE ( FIELD(1-NY:NY*(NX+2)) )
-      IF( RGLGRD .AND. .NOT. FLOMP ) ALLOCATE ( FIELD(1-NY:NY*(NX+2)) )
+      IF( RGLGRD ) ALLOCATE ( FIELD(1-NY:NY*(NX+2)) )
 !
 !/SMC !!Li   Otherwise use sea point only field
 !/SMC       ALLOCATE ( FIELD(NCel) )
@@ -1728,17 +1725,11 @@
 !/MPI         CALL MPI_STARTALL (NRQSG1, IRQSG1(1,2), IERR_MPI)
 !/MPI       END IF
 !
-!!/OMPX/!$OMP PARALLEL PRIVATE (ISPEC,FIELD)
-!
 !/DEBUGRUN        WRITE(740+IAPROC,*) 'W3WAVE, step 6.14'
 !/DEBUGRUN        FLUSH(740+IAPROC)
-
-            IF ( FLOMP ) ALLOCATE ( FIELD(1-NY:NY*(NX+2)) )
 !
 ! Initialize FIELD variable 
              FIELD = 0.
-!
-!!/OMPX/!$OMP DO SCHEDULE (DYNAMIC,1)
 !
             DO ISPEC=1, NSPEC
               IF ( IAPPRO(ISPEC) .EQ. IAPROC ) THEN
@@ -1785,14 +1776,7 @@
 !/MEMCHECK       write(740+IAPROC,*) 'memcheck_____:', 'WW3_WAVE TIME LOOP 17'
 !/MEMCHECK       call getMallocInfo(mallinfos)
 !/MEMCHECK       call printMallInfo(IAPROC,mallInfos)
-
-
-!!/OMPX/!$OMP END DO
-
-              IF ( FLOMP ) DEALLOCATE ( FIELD )
-
-!!/OMPX/!$OMP END PARALLEL
-
+!
 !Li   Initialise IK IX IY in case ARC option is not used to avoid warnings.
               IK=1
               IX=1
@@ -2400,7 +2384,7 @@
 !/OASIS ! Send variables to atmospheric or ocean circulation or ice model
 !/OASIS !
 !/OASIS                              IF (DTOUT(7).NE.0) THEN
-!/OASIS                                IF ( (MOD(ID_OASIS_TIME/DTOUT(7),1.0)  .LT. 1.E-7 ) .AND. &
+!/OASIS                                IF ( (MOD(ID_OASIS_TIME,NINT(DTOUT(7))) .EQ. 0 ) .AND. &
 !/OASIS                                     (DSEC21 (TIME00, TIME) .GT. 0.0) ) THEN
 !/OASACM                                  CALL SND_FIELDS_TO_ATMOS()
 !/OASOCM                                  CALL SND_FIELDS_TO_OCEAN()
@@ -2609,8 +2593,7 @@
 
       IF ( IAPROC .EQ. NAPLOG ) WRITE (NDSO,902)
 !
-      IF ( .NOT. FLOMP ) DEALLOCATE ( FIELD )
-!
+      DEALLOCATE(FIELD)
       DEALLOCATE(TAUWX, TAUWY)
 !
 !/MEMCHECK       write(740+IAPROC,*) 'memcheck_____:', 'WW3_WAVE END W3WAVE'

model/ftn/ww3_bound.ftn

@@ -12,7 +12,7 @@
 !/    28-Aug-2012 : adaptation from SHOM/Ifremer code   ( version 4.08 )
 !/    01-Nov-2012 : Bug correction for NKI != NK        ( version 4.08 )
 !/    20-Oct-2016 : Error statement updates             ( version 5.15 ) 
-!/    21-Jul-2020 : Support rotated pole grid           ( version 7.06 )
+!/    21-Jul-2020 : Support rotated pole grid           ( version 7.11 )
 !/                                  Chris Bunney, UKMO.
 !/
 !/    Copyright 2012-2012 National Weather Service (NWS),