!
! Toon Moene 20180913-15,17-21 - Initial version.
!            20180922          - Fix the domain decomposition and the
!                                Semi-Lagrangian code that depends on it.
!            20181229          - Switch to 0-based indexing for model arrays to
!                                make index conversion local <-> global simpler.
!            20190222          - Remove defined assignment; intrinsic assignment works.
!            20201004          - Reworked to use module procedures instead of type bound ones.
!                                This needs a defined assignment.
!            20201114          - Enforce coarrays to be equal size on every image;
!                                Stop if the number of images and domain size don't allow this.
!            20201128          - Correct synchronisation of images.
!            20210213          - Make subgrid scale physics access the whole atmosphere.
!            20250727          - Set the random seed on each image in such a way that
!                                different implementations will yield the same result numbers.
!
! Compile:
!    gfortran -fcoarray=lib    random-weather.f90 -lcaf_openmpi # Using OpenMPI's Message Passing Implementation
!    gfortran -fcoarray=lib    random-weather.f90 -lcaf_shmem   # Using "Vehreschild" shared memory method
!    gfortran -fcoarray=shared random-weather.f90 -lcaf_shared  # Using "Koenig" shared memory method
!
! Run:
!    echo ' &config / ' | mpirun -np 9 ./a.out
!    (export GFORTRAN_NUM_IMAGES=9; echo ' &config / ' | ./a.out)
!
! A mock "weather" forecasting program that goes through the motions
! of "representing" the weather on a limited area spread over several
! images to show the use of coarrays.
!
! Why are weather forecasting systems multi-million-line, multi-language
! beasts when you can do it with a simple, 500 line, Fortran program ?
!
! Because *a lot* is left out below:
!
! A. Semi-Lagrangian Advection.
!
!    The atmosphere is not a gas in a stationary rectangular box,
!    but resides in a curvilinear, rotating coordinate system, with
!    an uneven surface (mountains).
!
!    It is subject to forces: gravity, and the pseudo forces due to
!    rotation: the centrifugal and Coriolis force.
!
!    (Its mass is cleverly hidden in the surface pressure, along with
!    the - just as implicit - assumption of hydrostatic equilibrium).
!
! B. "Random" sub grid scale variation.
!
!    These adjustments of the rate of change of atmospheric quantities are
!    certainly not "random":
!
!    1. Irradiation by the Sun and returned infrared radiation, and its
!       interaction with (trace-)gasses, clouds and aerosols.
!    2. Conversion between the three phases of water: vapor, liquid and solid,
!       and the resulting heat exchange and release of precipitation.
!    3. Differential heating of land and sea surface.
!    4. Other interaction with the surface: friction, saturation, runoff, evaporation.
!    5. Turbulence.
!
! C. Adjustment of atmospheric quantitites at the boundaries with input from a global
!    model, whose data has to be read in and interpolated to the model's own grid.
!
! D. Adjustment of the weather development in the model due to new information
!    coming in from various observations: land stations, ships, balloons,
!    aircraft, radar and satellites.
!
module constants
!
!       1. Physics
!                        Pascal      Kelvin       m/s         m/s         m/s         kg/kg
real, parameter :: PS = 100000.0, T = 300.0,    U = 0.0,    V = 0.0,    W = 0.0,      Q = 0.002      ! Mean value
real, parameter :: dPSdt = 0.1,  dTdt = 0.002, dUdt = 0.1, dVdt = 0.1, dWdt = 0.001, dQdt = 0.000001 ! Maximum tendencies.
!
!       2. Domain size and discretization; forecast time and step
!                      number    number    number    number       meters           meters        seconds       seconds
integer, parameter :: DNX = 72, DNY = 70, DNZ = 30, BDSIZE = 4, HORSTEP = 10000, VERSTEP = 100, FCLEN = 3600, TIMSTEP = 240
!
end module constants
!
module grid
!
! Partition the work into "slabs" of the computational domain, which is
! defined by the discretization of the partial differential equations
! resulting from the physical laws governing atmospheric motion.
! Note that the code below makes extensive use of the fact that a coarray/coscalar *without*
! coindices refers to its part on the local slab (called an image in Fortran).
!
   integer :: nxglobal[*], nyglobal[*], nzglobal[*]                         ! Extent of the entire grid
   integer :: hgrid_step[*], vgrid_step[*], boundary_size[*], time_step[*]  ! Grid steps (meters), etc.
   integer :: nxslab[*], nyslab[*], npx[*], npy[*]                          ! Extent and number of the horizontal slabs
!
end module grid
!
module types
!
! Declaration of the types '{global|local}_model_state' and their module procedures
!
   type global_model_state
      real, allocatable :: ps(:,:)  [:,:] ! Surface pressure
      real, allocatable :: t(:,:,:) [:,:] ! Temperature
      real, allocatable :: u(:,:,:) [:,:] ! U-component of wind (along parallels)
      real, allocatable :: v(:,:,:) [:,:] ! V-component of wind (along meridians)
      real, allocatable :: w(:,:,:) [:,:] ! W-component of wind (vertical)
      real, allocatable :: q(:,:,:) [:,:] ! Specific humidity
   end type global_model_state
!
   type local_model_state
      real, allocatable :: ps(:,:)  ! Surface pressure
      real, allocatable :: t(:,:,:) ! Temperature
      real, allocatable :: u(:,:,:) ! U-component of wind (along parallels)
      real, allocatable :: v(:,:,:) ! V-component of wind (along meridians)
      real, allocatable :: w(:,:,:) ! W-component of wind (vertical)
      real, allocatable :: q(:,:,:) ! Specific humidity
   contains
      final :: dealloc              ! Deallocate a local model state that goes out of scope
   end type local_model_state
!
   interface operator(+)
      module procedure add                ! Standard
   end interface                          !  Overloading
   interface operator(*)                  !   Of
      module procedure int_mult           !    Operators
   end interface                          !     And
   interface assignment(=)                !      Assignment
      module procedure copy_local_to_global
   end interface
   interface operator(.phys.)             ! Operator for
      module procedure phys               !  subgrid scale
   end interface                          !   physics
   interface operator(.sla.)              ! Operator for
      module procedure sla                !  Semi-Lagrangian
   end interface                          !   Advection
!
contains
!
! Implementation of the module procedures on types '{global|local}_model_state'.
!
   subroutine dealloc(lms)
      type(local_model_state) :: lms
      if (allocated(lms % ps)) deallocate(lms % ps)
      if (allocated(lms % t )) deallocate(lms % t)
      if (allocated(lms % u )) deallocate(lms % u)
      if (allocated(lms % v )) deallocate(lms % v)
      if (allocated(lms % w )) deallocate(lms % w)
      if (allocated(lms % q )) deallocate(lms % q)
   end subroutine dealloc
!
   subroutine global_alloc(gms, nx, ny, nz, npx)
      type(global_model_state) :: gms
      integer, intent(in) :: nx, ny, nz, npx
      allocate( gms % ps(0:nx-1, 0:ny-1)        [0:npx-1, 0:*] )
      allocate( gms % t(0:nx-1, 0:ny-1, 0:nz-1) [0:npx-1, 0:*] )
      allocate( gms % u(0:nx-1, 0:ny-1, 0:nz-1) [0:npx-1, 0:*] )
      allocate( gms % v(0:nx-1, 0:ny-1, 0:nz-1) [0:npx-1, 0:*] )
      allocate( gms % w(0:nx-1, 0:ny-1, 0:nz-1) [0:npx-1, 0:*] )
      allocate( gms % q(0:nx-1, 0:ny-1, 0:nz-1) [0:npx-1, 0:*] )
   end subroutine global_alloc
!
   subroutine local_alloc(lms, nx, ny, nz)
      type(local_model_state) :: lms
      integer, intent(in) :: nx, ny, nz
      allocate( lms % ps(0:nx-1, 0:ny-1)        )
      allocate( lms % t(0:nx-1, 0:ny-1, 0:nz-1) )
      allocate( lms % u(0:nx-1, 0:ny-1, 0:nz-1) )
      allocate( lms % v(0:nx-1, 0:ny-1, 0:nz-1) )
      allocate( lms % w(0:nx-1, 0:ny-1, 0:nz-1) )
      allocate( lms % q(0:nx-1, 0:ny-1, 0:nz-1) )
   end subroutine local_alloc
!
   subroutine global_init(gms, ps, t, u, v, w, q)
      type(global_model_state), intent(inout) :: gms
      real, intent(in) :: ps, t, u, v, w, q
      gms % ps = ps
      gms % t  = t 
      gms % u  = u 
      gms % v  = v 
      gms % w  = w 
      gms % q  = q 
   end subroutine global_init
!
   subroutine copy_local_to_global(gms, lms)
      type(global_model_state), intent(inout) :: gms
      type(local_model_state),  intent(in)    :: lms
      sync all ! Before overwriting the global model state,
               ! ensure all its uses are complete.
      gms % ps = lms % ps
      gms % t  = lms % t
      gms % u  = lms % u
      gms % v  = lms % v
      gms % w  = lms % w
      gms % q  = lms % q
   end subroutine copy_local_to_global
!
   function phys(gms)
      use constants, only : dPSdt, dTdt, dUdt, dVdt, dWdt, dQdt
      type(local_model_state)              :: phys
      type(global_model_state), intent(in) :: gms
      integer :: nx, ny, nz
      nx = size(gms % t, 1)
      ny = size(gms % t, 2)
      nz = size(gms % t, 3)
      call local_alloc(phys, nx, ny, nz)
!
! Compute the subgrid scale physical tendencies (random fraction of maximum, for now).
! Note that the full model atmosphere is the argument to this function, so it is possible
! to compute tendencies on the current image due to physical processes *outside* it
! (slant solar radiation or slant precipitation, three dimensional turbulence).
! See the "semi-langrangian" function for how slant access is done.
!
      call random_number(phys % ps); phys % ps = (1.0 - 2.0 * phys % ps) * dPSdt ! Pa/s
      call random_number(phys % t ); phys % t  = (1.0 - 2.0 * phys % t ) * dTdt  ! K/s
      call random_number(phys % u ); phys % u  = (1.0 - 2.0 * phys % u ) * dUdt  ! m/s^2
      call random_number(phys % v ); phys % v  = (1.0 - 2.0 * phys % v ) * dVdt  ! m/s^2
      call random_number(phys % w ); phys % w  = (1.0 - 2.0 * phys % w ) * dWdt  ! m/s^2
      call random_number(phys % q ); phys % q  = (1.0 - 2.0 * phys % q ) * dQdt  ! kg/kg/s
   end function phys
!
   function add(lms1, lms2)
      type(local_model_state)             :: add
      type(local_model_state), intent(in) :: lms1, lms2
      integer :: nx, ny, nz
      nx  = size(lms1 % t, 1)
      ny  = size(lms1 % t, 2)
      nz  = size(lms1 % t, 3)
      call local_alloc(add, nx, ny, nz)
      add % ps = lms1 % ps + lms2 % ps
      add % t  = lms1 % t  + lms2 % t 
      add % u  = lms1 % u  + lms2 % u 
      add % v  = lms1 % v  + lms2 % v 
      add % w  = lms1 % w  + lms2 % w 
      add % q  = lms1 % q  + lms2 % q 
   end function add
!
   function int_mult(lms, ifactor)
      type(local_model_state)             :: int_mult
      type(local_model_state), intent(in) :: lms
      integer, intent(in) :: ifactor
      integer :: nx, ny, nz
      nx  = size(lms % t, 1)
      ny  = size(lms % t, 2)
      nz  = size(lms % t, 3)
      call local_alloc(int_mult, nx, ny, nz)
      int_mult % ps = lms % ps * ifactor
      int_mult % t  = lms % t  * ifactor 
      int_mult % u  = lms % u  * ifactor
      int_mult % v  = lms % v  * ifactor
      int_mult % w  = lms % w  * ifactor
      int_mult % q  = lms % q  * ifactor
   end function int_mult
!
   function sla(gms)
      use grid
      type(local_model_state)              :: sla
      type(global_model_state), intent(in) :: gms
      integer :: nx, jx, ny, jy, nz, jz, ico(2), ioff(3), idep(3)
      nx  = size(gms % t, 1)
      ny  = size(gms % t, 2)
      nz  = size(gms % t, 3)
      ico = this_image(gms % ps) ! The coindices of the model coarrays on this image.
      call local_alloc(sla, nx, ny, nz)
!
!        Semi-Lagrangian advection looks upstream (i.e., up*wind* - to the departure point)
!        at the atmospheric quantities to determine their values in the arrival point.
!        The real thing is more complicated than just transporting the values, but it
!        gives an idea of the (co)index computations and *communications* involved.
!
      do jz = 0, nz-1
         do jy = 0, ny-1
            do jx = 0, nx-1
!
!        Compute the two horizontal global indices of the departure point based on the local upwind speed.
!
               ioff(1) = nint(gms % u(jx, jy, jz) * time_step / hgrid_step)       ! Upwind X distance in grid cells ...
               idep(1) = max(0, min(nxglobal-1, jx + nxslab * ico(1) - ioff(1)))  ! determines X departure global index.
               ioff(2) = nint(gms % v(jx, jy, jz) * time_step / hgrid_step)       ! Upwind Y distance in grid cells ...
               idep(2) = max(0, min(nyglobal-1, jy + nyslab * ico(2) - ioff(2)))  ! determines Y departure global index.
               ioff(3) = nint(gms % w(jx, jy, jz) * time_step / vgrid_step)       ! Upwind Z distance in grid cells ...
               idep(3) = max(0, min(nzglobal-1, jz                   - ioff(3)))  ! determines Z departure global index.
!
!        The departure point is accessed by splitting its two horizontal global indices into two local ones and two coindices.
!
               if (jz == 0) &
                  sla % ps(jx, jy) = gms % ps(mod(idep(1),nxslab), mod(idep(2),nyslab)         ) [idep(1)/nxslab, idep(2)/nyslab]
               sla % t(jx, jy, jz) = gms % t (mod(idep(1),nxslab), mod(idep(2),nyslab), idep(3)) [idep(1)/nxslab, idep(2)/nyslab]
               sla % u(jx, jy, jz) = gms % u (mod(idep(1),nxslab), mod(idep(2),nyslab), idep(3)) [idep(1)/nxslab, idep(2)/nyslab]
               sla % v(jx, jy, jz) = gms % v (mod(idep(1),nxslab), mod(idep(2),nyslab), idep(3)) [idep(1)/nxslab, idep(2)/nyslab]
               sla % w(jx, jy, jz) = gms % w (mod(idep(1),nxslab), mod(idep(2),nyslab), idep(3)) [idep(1)/nxslab, idep(2)/nyslab]
               sla % q(jx, jy, jz) = gms % q (mod(idep(1),nxslab), mod(idep(2),nyslab), idep(3)) [idep(1)/nxslab, idep(2)/nyslab]
            enddo
         enddo
      enddo
   end function sla
!
   subroutine bndrel(gms, ps, t, u, v, w, q)
      use grid
      type(global_model_state), intent(inout) :: gms
      real,                     intent(in)    :: ps, t, u, v, w, q
      real    :: alpha
      integer :: nx, jx, ny, jy, nz, jz, ibegin, iend, ico(2), uco(2)
      nx = size(gms % t, 1)
      ny = size(gms % t, 2)
      nz = size(gms % t, 3)
      ico = this_image(gms % ps) ! The coindices of the model coarrays on this image.
      uco = ucobound(gms % ps)
!
!        Boundary relaxation is done to keep the "weather" inside the limited domain
!        current with the developments outside it. In this "mock" weather forecasting
!        program, we relax to fixed mean values.
!
!        To prevent the relaxation to be doubly performed at the corners,
!        the following convention is used:
!
!               North
!             |-------
!        West |      | East
!             -------|
!               South
!
!        The origin of the grid is in the South-West corner.
!
!  1. Relax (part of) the South boundary, if present in this slab
      if (ico(2) == 0) then                ! The South boundary (y=0)
         if (ico(1) == uco(1)) then        ! at the East side (x=max)
            iend = nx - 1 - boundary_size
         else                              ! ... or if not ...
            iend = nx - 1
         endif
         do jz = 0, nz-1
            do jy = 0, boundary_size-1
               alpha = real(boundary_size - jy) / boundary_size
               do jx = 0, iend
                  if (jz == 0) &
                     gms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * gms % ps(jx, jy)
                  gms % t(jx, jy, jz) = alpha * t  + (1.0 - alpha) * gms % t(jx, jy, jz)
                  gms % u(jx, jy, jz) = alpha * u  + (1.0 - alpha) * gms % u(jx, jy, jz)
                  gms % v(jx, jy, jz) = alpha * v  + (1.0 - alpha) * gms % v(jx, jy, jz)
                  gms % w(jx, jy, jz) = alpha * w  + (1.0 - alpha) * gms % w(jx, jy, jz)
                  gms % q(jx, jy, jz) = alpha * q  + (1.0 - alpha) * gms % q(jx, jy, jz)
               enddo
            enddo
         enddo
      endif
!  2. Relax (part of) the North boundary, if present in this slab
      if (ico(2) == uco(2)) then           ! The North boundary (y=max)
         if (ico(1) == 0) then             ! at the West side (x=0)
            ibegin = boundary_size
         else                              ! ... or if not ...
            ibegin = 0
         endif
         do jz = 0, nz-1
            do jy = ny-1, ny-1 - boundary_size
               alpha = real(jy - ny-1 + boundary_size) / boundary_size
               do jx = ibegin, nx-1
                  if (jz == 0) &
                     gms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * gms % ps(jx, jy)
                  gms % t(jx, jy, jz) = alpha * t  + (1.0 - alpha) * gms % t(jx, jy, jz)
                  gms % u(jx, jy, jz) = alpha * u  + (1.0 - alpha) * gms % u(jx, jy, jz)
                  gms % v(jx, jy, jz) = alpha * v  + (1.0 - alpha) * gms % v(jx, jy, jz)
                  gms % w(jx, jy, jz) = alpha * w  + (1.0 - alpha) * gms % w(jx, jy, jz)
                  gms % q(jx, jy, jz) = alpha * q  + (1.0 - alpha) * gms % q(jx, jy, jz)
               enddo
            enddo
         enddo
      endif
!  3. Relax (part of) the West boundary, if present in this slab
      if (ico(1) == 0) then                ! The West boundary (x=0)
         if (ico(2) == 0) then             ! at the South side (y=0)
            ibegin = boundary_size
         else                              ! ... or if not ...
            ibegin = 0
         endif
         do jz = 0, nz-1
            do jy = ibegin, ny-1
               do jx = 0, boundary_size-1
                  alpha = real(boundary_size - jx) / boundary_size
                  if (jz == 0) &
                     gms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * gms % ps(jx, jy)
                  gms % t(jx, jy, jz) = alpha * t  + (1.0 - alpha) * gms % t(jx, jy, jz)
                  gms % u(jx, jy, jz) = alpha * u  + (1.0 - alpha) * gms % u(jx, jy, jz)
                  gms % v(jx, jy, jz) = alpha * v  + (1.0 - alpha) * gms % v(jx, jy, jz)
                  gms % w(jx, jy, jz) = alpha * w  + (1.0 - alpha) * gms % w(jx, jy, jz)
                  gms % q(jx, jy, jz) = alpha * q  + (1.0 - alpha) * gms % q(jx, jy, jz)
               enddo
            enddo
         enddo
      endif
!  4. Relax (part of) the East boundary, if present in this slab
      if (ico(1) == uco(1)) then           ! The East boundary (x=max)
         if (ico(2) == uco(2)) then        ! at the North end (y=max)
            iend = ny - boundary_size - 1
         else                              ! ... or if not ...
            iend = ny - 1
         endif
         do jz = 0, nz-1
            do jy = 0, iend
               do jx = nx-1, nx-1 - boundary_size
                  alpha = real(jx - nx-1 + boundary_size) / boundary_size
                  if (jz == 0) &
                     gms % ps(jx, jy) = alpha * ps + (1.0 - alpha) * gms % ps(jx, jy)
                  gms % t(jx, jy, jz) = alpha * t  + (1.0 - alpha) * gms % t(jx, jy, jz)
                  gms % u(jx, jy, jz) = alpha * u  + (1.0 - alpha) * gms % u(jx, jy, jz)
                  gms % v(jx, jy, jz) = alpha * v  + (1.0 - alpha) * gms % v(jx, jy, jz)
                  gms % w(jx, jy, jz) = alpha * w  + (1.0 - alpha) * gms % w(jx, jy, jz)
                  gms % q(jx, jy, jz) = alpha * q  + (1.0 - alpha) * gms % q(jx, jy, jz)
               enddo
            enddo
         enddo
      endif
   end subroutine bndrel
!
end module types
!
program random_weather
!
! 0. Declarations.
!
!    0.1. Modules
! 
use constants, only : &
   DNX, DNY, DNZ, HORSTEP, VERSTEP, BDSIZE, FCLEN, TIMSTEP, &
   PS, T, U, V, W, Q
use grid
use types
!
implicit none
!
!    0.2. Configuration variables
!
integer :: forecast_length[*]
namelist /config/ nxglobal, nyglobal, nzglobal, boundary_size, hgrid_step, vgrid_step, forecast_length, time_step
!
!    0.3. Model state variables
!
type(global_model_state) :: atmosphere
!
!    0.4. Local variables
!
integer :: i, time, np1, np2, npxy, nseeds
integer, allocatable :: iseeds(:)
!
! 1. Determine the number of images and compute the domain decomposition.
!
if (this_image() == 1) then
!
!    1.1. First, set the domain definition defaults.
!
   nxglobal = DNX; nyglobal = DNY; nzglobal = DNZ; hgrid_step = HORSTEP; vgrid_step = VERSTEP
   boundary_size = BDSIZE; forecast_length = FCLEN; time_step = TIMSTEP
!
!    1.2. Read the domain definition, overriding the defaults.
!
   read(*,config)
!
!    1.3. Check whether equal sized slabs are possible with this choice of
!         number of images and domain size.
!
   npxy = num_images()
   if (mod(nxglobal * nyglobal, npxy) /= 0) then
      print*, 'Number of images and domain size do not allow an equal-sized partitioning.'
      stop 1
   endif
!
!    1.4. Compute the domain decomposition. Note: only horizontal decomposition.
!         The decomposition should have an "aspect ratio" close to that of the grid
!         to get maximum load balancing, i.e.: npx / npy = grid_x / grid_y.
!         Multiply both sides with number_of_images (npx * npy) -> npx^2 = (grid_x/grid_y)*n_o_i.
!
   np1 = int(sqrt(real(nxglobal*npxy/nyglobal)))
   do
      np2 = npxy / np1
      if (np1 * np2 == npxy) then
         if (mod(nxglobal, np1) == 0 .and. mod(nyglobal, np2) == 0) then
            npx = np1; npy = np2; exit
         endif
         if (mod(nxglobal, np2) == 0 .and. mod(nyglobal, np1) == 0) then
            npx = np2; npy = np1; exit
         endif
      endif
      np1 = np1 + 1
   enddo
   nxslab = nxglobal / npx
   nyslab = nyglobal / npy
!
!    1.5. Check that the size of the boundary relaxation zone doesn't
!         exceed the size of the slabs.
!
   if      (boundary_size > nxslab) then
      print*,'Size of boundary zone larger than x side of slabs'
      stop 1
   else if (boundary_size > nyslab) then
      print*,'Size of boundary zone larger than y side of slabs'
      stop 1
   endif
!
!    1.6. Then distribute the resulting definition to the other images.
!
!     Why won't this work as nxglobal[2:] = nxglobal ?
   do i = 2, npxy
      nxglobal[i] = nxglobal; nyglobal[i] = nyglobal; nzglobal[i] = nzglobal; boundary_size[i] = boundary_size
      hgrid_step[i] = hgrid_step; vgrid_step[i] = vgrid_step; forecast_length[i] = forecast_length; time_step[i] = time_step
      nxslab[i] = nxslab; nyslab[i] = nyslab; npx[i] = npx; npy[i] = npy
   enddo
endif
!
!    1.7. Ensure all images have the domain definition information before proceeding.
!
sync all
!
!    1.8. Ensure repeatable random numbers that are distinct on each image.
!
call random_seed(size = nseeds)
allocate(iseeds(nseeds))
iseeds = this_image()
call random_seed(put = iseeds)
!
print '(5(a,i4),a)', &
   'Decomposition information on image',this_image(),' : there are', &
   npx,' *',npy,' slabs; the slabs are',nxslab,' *',nyslab,' grid cells in size.'
!
! 2. Allocate and "compute" the initial state of the model "atmosphere".
!
!    2.1. The initial values of the atmospheric quantities.
!
call global_alloc(atmosphere, nxslab, nyslab, nzglobal, npx)
call global_init (atmosphere, PS, T, U, V, W, Q)
!
!    2.2. Ensure atmosphere is initialised on all images before entering the loop.
!
sync all
!
! 3. "Integrate" the model in time.
!
do time = 0, forecast_length, time_step
!
!    3.1. Perform "semi-lagrangian" advection (i.e., downstream)
!         and add the increments due to subgrid scale physics.
!
   atmosphere = .sla. atmosphere + .phys. atmosphere * time_step
!
!    Technical note: All images are in sync when the copy to "atmosphere" is excuted by
!                    the defined assignment procedure copy_local_to_global due to
!                    the "sync all" in it assuring that.
!                    Because the subsequent statements in this loop only update the
!                    local part of atmosphere, no further synchronisation is necessary
!                    until all updates are done (after boundary relaxation).
!
!    3.2. Perform boundary relaxation towards mean values (to keep this simple).
!
   call bndrel(atmosphere, PS, T, U, V, W, Q)
   sync all
!
!    3.3. Print some diagnostic model output.
!
   print*, 'Time',time,' Image',this_image(),         &
      ' PS=', atmosphere % ps(nxslab/2, nyslab/2),    &
      ' T= ', atmosphere % t (nxslab/2, nyslab/2, 0), &
      ' U= ', atmosphere % u (nxslab/2, nyslab/2, 0), &
      ' V= ', atmosphere % v (nxslab/2, nyslab/2, 0), &
      ' W= ', atmosphere % w (nxslab/2, nyslab/2, 0), &
      ' Q= ', atmosphere % q (nxslab/2, nyslab/2, 0)
!
end do
!
! 4. Complete the program.
!
end program random_weather