fabm_update Subroutine

      subroutine fabm_update(joff, js, je, is, ie)
         !! update the environment and calculate the source/sink terms -
         !! is called with the same argument list as mom() calls tracer()
         !! i.e. the specification of the active UVic window - typically
         !! the full domain on modern hardware
      integer, intent(in) :: joff
         !! offset row in global window
      integer, intent(in) :: js
         !! start row
      integer, intent(in) :: je
         !! end row
      integer, intent(in) :: is
         !! start column
      integer, intent(in) :: ie
         !! end column

         integer :: i,j,k,n
            ! local loop counters
         real :: wloc
         real :: flux

#ifdef DEBUG
        !print*, 'fabm_update:',joff,js,je,is,ie
        !print*, 'fabm_update:',t(53,1,53,itemp,1),src(53,1,53,3:6)
        !print*, 'fabm_update:',t(53,1,53,itemp,1),t(53,1,53,3:6,0)
#endif
         surface_flux = 0._rke
         surface_sms = 0._rke
         pelagic_sms = 0._rke
         bottom_flux = 0._rke
         bottom_sms = 0._rke
         src = 0._rke

         ! t(:,:,:,var,0) is updated in loadmw() in mom()
         ! this is done before the call to tracer() - and thus
         ! data are ready here
         call update_data(joff)

#if defined O_fabm_check_state
         if (nsurface > 0) then
            do j=js,je
               call model%check_surface_state(is,ie,j,repair,
     &                                        valid_surf)
            end do
         end if
         if (npelagic > 0) then
            do j=js,je
               do k=1,km
                  call model%check_interior_state(is,ie,k,j,repair,
     &                                            valid_pel)
               end do
            end do
         end if
         if (nbottom > 0) then
            do j=js,je
               call model%check_bottom_state(is,ie,j,repair,
     &                                       valid_bott)
            end do
         end if
#endif

         call model%prepare_inputs()

         ! update the surface
         if (nsurface > 0) then
            do j=js,je
               call model%get_surface_sources(is,ie,j,
     &              surface_flux(is:ie,j,:),surface_sms(is:ie,j,:))
            end do
         end if

         ! update the pelagic
         do j=js,je
            pelagic_sms = 0._rke
            do k=1,km
               call model%get_interior_sources(is,ie,k,j,
     &              pelagic_sms(is:ie,k,1,:))
            end do
            src(is:ie,:,j,offset:) = pelagic_sms(is:ie,:,1,:)
         end do

         ! update the bottom
         if (nbottom > 0) then
            do j=js,je
               call model%get_bottom_sources(is,ie,j,
     &              bottom_flux(is:ie,j,:),bottom_sms(is:ie,j,:))
            end do
         end if

         ! fold the surface and bottom flux terms
         do j=js,je
            do i=is,ie
               if (kmt(i,j) > 0) then
                  k=1 ! surface
                  src(i,k,j,offset:)=src(i,k,j,offset:)+
     &                          surface_flux(i,j,:)/dz(i,k,j)
                  k=kmt(i,j) ! bottom
                  src(i,k,j,offset:)=src(i,k,j,offset:)+
     &                          bottom_flux(i,j,:)/dz(i,k,j)
               end if
            end do
         end do

         ! vertical velocities
         do j=js,je
            do k=1,km
               call model%get_vertical_movement(is,ie,k,j,w(is:ie,k,:))
            end do
            ! do vertical advection - first-order upstream
            do n=1,npelagic
               if ( .not. (any(w(is:ie,:,n) /= 0.0_rke))) cycle
               do i=is,ie
                  do k=1,kmt(i,j)-1
                     wloc = -0.5_rke*(w(i,k,n) + w(i,k+1,n))
                     if (wloc > 0.0_rke) then
                        flux = wloc*t(i,k  ,j,n,0)
                     else
                        flux = wloc*t(i,k+1,j,n,0)
                     end if
                     t(i,k  ,j,n,0) = t(i,k  ,j,n,0) - flux/dz(i,k  ,j)
                     t(i,k+1,j,n,0) = t(i,k+1,j,n,0) + flux/dz(i,k+1,j)
                  end do
               end do
            end do
         end do

         call model%finalize_outputs()
      end subroutine fabm_update