#include "ATM_PHYS_OPTIONS.h"

CBOP
C !ROUTINE: ATM_PHYS_DRIVER

C !INTERFACE: ==========================================================
      SUBROUTINE ATM_PHYS_DRIVER(
     I                    myTime, myIter, myThid )

C !DESCRIPTION:
C Calculate custom tendency terms outside k-loop in DO_OCEANIC_PHYS

C !USES: ===============================================================
      use radiation_mod
      use lscale_cond_mod
      use dargan_bettsmiller_mod
      use surface_flux_mod
      use vert_turb_driver_mod
      use vert_diff_mod, only: gcm_vert_diff_down,
     &                         gcm_vert_diff_up,
     &                         surf_diff_type
      use mixed_layer_mod, only: mixed_layer
      use constants_mod, only:  HLv

      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"
#include "ATM_PHYS_PARAMS.h"
#include "ATM_PHYS_VARS.h"

C !INPUT PARAMETERS: ===================================================
C  myTime   :: Current time in simulation
C  myIter   :: Current time-step number
C  myThid   :: my Thread Id number
      _RL     myTime
      INTEGER myIter, myThid

C !OUTPUT PARAMETERS: ==================================================

C !LOCAL VARIABLES: ====================================================
C  bi,bj    :: Tile indices
C  lat2d    :: latitude of grid-cell center          [rad]
C pHalf3d   :: pressure at interface between 2 levels [Pa]
C pFull3d   :: pressure at level center               [Pa]
C zHalf3d   :: height of interface between 2 levels   [m]
C zFull3d   :: height of level center                 [m]
C  t3d      :: absolute temperature                   [K]
C  q3d      :: specific humidity                    [kg/kg]
C  u3d      :: wind speed, 1rst component (X-dir)    [m/s]
C  v3d      :: wind speed, 2nd  component (Y-dir)    [m/s]
      INTEGER bi, bj
      _RL lat2d   (sNx,sNy)
      _RL pHalf3d (sNx,sNy,Nr+1)
      _RL pFull3d (sNx,sNy,Nr)
      _RL zHalf3d (sNx,sNy,Nr+1)
      _RL zFull3d (sNx,sNy,Nr)
      _RL t3d     (sNx,sNy,Nr)
      _RL q3d     (sNx,sNy,Nr)
      _RL u3d     (sNx,sNy,Nr)
      _RL v3d     (sNx,sNy,Nr)
      _RL tdt3d   (sNx,sNy,Nr)
      _RL qdt3d   (sNx,sNy,Nr)
      _RL udt3d   (sNx,sNy,Nr)
      _RL vdt3d   (sNx,sNy,Nr)
CEOP
      _RL s_sw_dwn(sNx,sNy)
      _RL s_lw_dwn(sNx,sNy)
      _RL t_surf  (sNx,sNy)
C-- radiation fields:
      _RL albedo_2d (sNx,sNy)
      _RL ozone_3d  (sNx,sNy,Nr)
      _RL ozoneColmn(sNx,sNy,Nr+1)
      _RL dtrans_3d (sNx,sNy,Nr)
      _RL dtrans_win(sNx,sNy,Nr)
      _RL b_3d      (sNx,sNy,Nr)
      _RL b_win     (sNx,sNy,Nr)
      _RL lw_down_3d(sNx,sNy,Nr+1)
      _RL sw_down_3d(sNx,sNy,Nr+1)
      _RL sw_net_3d (sNx,sNy,Nr+1)
      _RL lw_net_3d (sNx,sNy,Nr+1)
      _RL adj_lw_up (sNx,sNy)
      _RL rad_dt_tg (sNx,sNy,Nr)

      LOGICAL coldT (sNx,sNy)
C-- output from convection & LSC :
      _RL t3d_tmp    (sNx,sNy,Nr)
      _RL q3d_tmp    (sNx,sNy,Nr)
      _RL cond_dt_tg (sNx,sNy,Nr)
      _RL cond_dt_qg (sNx,sNy,Nr)
      _RL rain2d     (sNx,sNy)
      _RL snow2d     (sNx,sNy)
      _RL q_ref      (sNx,sNy,Nr)
      _RL t_ref      (sNx,sNy,Nr)
      _RL bmflag     (sNx,sNy)
      _RL klzbs      (sNx,sNy)
      _RL cape       (sNx,sNy)
      _RL cin        (sNx,sNy)
      _RL invtau_bm_t(sNx,sNy)
      _RL invtau_bm_q(sNx,sNy)
      _RL capeflag   (sNx,sNy)
C-- Input/Output for surface flux:
      _RL  q_surf(sNx,sNy)
      _RL  u_surf(sNx,sNy), v_surf(sNx,sNy)
      _RL  rough_mom(sNx,sNy), rough_heat(sNx,sNy)
      _RL  rough_moist(sNx,sNy), gust(sNx,sNy)
      _RL  flux_t(sNx,sNy), flux_q(sNx,sNy), flux_r(sNx,sNy)
      _RL  flux_u(sNx,sNy), flux_v(sNx,sNy)
      _RL  drag_m(sNx,sNy), drag_t(sNx,sNy), drag_q(sNx,sNy)
      _RL  w_atm(sNx,sNy)
      _RL  ustar(sNx,sNy), bstar(sNx,sNy), qstar(sNx,sNy)
      _RL  dhdt_surf(sNx,sNy), dedt_surf(sNx,sNy), dedq_surf(sNx,sNy)
      _RL  drdt_surf(sNx,sNy)
      _RL  dhdt_atm(sNx,sNy), dedq_atm(sNx,sNy), dtaudv_atm(sNx,sNy)
      LOGICAL land(sNx,sNy), avail(sNx,sNy)
C-- Input for turb:
      _RL fracland(sNx,sNy)
      _RL rough(sNx,sNy)
C-- Output from turbulence driver:
      _RL  diff_t(sNx,sNy,Nr), diff_m(sNx,sNy,Nr)
      _RL diff_dt_tg (sNx,sNy,Nr)
      _RL diff_dt_qg (sNx,sNy,Nr)
      _RL diss_heat  (sNx,sNy,Nr)
c     TYPE(surf_diff_type) :: tri_surf ! used by gcm_vert_diff
      _RL tri_surf_dtmass(sNx,sNy)
      _RL tri_surf_dflux_t(sNx,sNy), tri_surf_dflux_q(sNx,sNy)
      _RL tri_surf_delta_t(sNx,sNy), tri_surf_delta_q(sNx,sNy)
      _RL e_global(sNx,sNy,Nr-1)
      _RL f_t_global(sNx,sNy,Nr-1), f_q_global(sNx,sNy,Nr-1)
C-- Mixed Layer fields:
      _RL ocean_qflux(sNx,sNy)
      _RL mixLayDepth(sNx,sNy)
      _RL delta_t_surf(sNx,sNy)
C--
      _RL dpFac(sNx,sNy)
      _RL conv_T2theta
      INTEGER k, kc
c     INTEGER ioUnit
c     _RS     dummyRS(1)
#ifdef COMPONENT_MODULE
      INTEGER i, j
      _RL taux(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL tauy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#endif /* COMPONENT_MODULE */

C--   Loops on tile indices bi,bj
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

C--   Get surface temperature
        t_surf(:,:) = atmPhys_SST(1:sNx,1:sNy,bi,bj)
#ifdef COMPONENT_MODULE
        IF ( useCoupler ) THEN
C--   take surface data from the ocean component
C     to replace MxL fields (if use sea-ice) or directly atm_phys SST
          CALL ATM_APPLY_IMPORT(
     I             maskInC,
     U             t_surf, taux,
c    U             sst1(1,myThid), oice1,
     I             myTime, myIter, bi, bj, myThid )
        ENDIF
#endif /* COMPONENT_MODULE */
#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
C--   fill-up Atm_Phys diagnostics for state variables
         CALL DIAGNOSTICS_FILL( t_surf, 'AtPh_SST',
     &                          0, 1, 3, bi, bj, myThid )
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */
c       IF ( myIter.EQ.nIter0 ) THEN
c        ioUnit = 0
c        CALL MDS_WRITEVEC_LOC(
c    I                     'SST.check', writeBinaryPrec, ioUnit,
c    I                     'RL', sNx*sNy, t_surf, dummyRS,
c    I                     bi, bj, 1, myIter, myThid )
c       ENDIF
        ocean_qflux(:,:) = atmPhys_Qflx(1:sNx,1:sNy,bi,bj)
        mixLayDepth(:,:) = atmPhys_MxLD(1:sNx,1:sNy,bi,bj)
        albedo_2d(:,:) = atmPhys_Albedo(1:sNx,1:sNy,bi,bj)
C--   copy atmPhys_Ozone to ozone_3d for rad scheme
        DO k=1,Nr
          kc = Nr-k+1
          ozone_3d(:,:,k) = atmPhys_Ozone(1:sNx,1:sNy,kc,bi,bj)
        ENDDO

C--   Get grid and dynamical fields from main model common blocks
        CALL ATM_PHYS_DYN2PHYS(
     O                    lat2d, pHalf3d, pFull3d,
     O                    zHalf3d, zFull3d,
     O                    t3d, q3d, u3d, v3d,
     I                    bi, bj, myTime, myIter, myThid )

C--   initialise
        coldT(:,:) = .FALSE.
        EmPmR(:,:,bi,bj) = 0.
        rain2d(:,:) = 0.
        snow2d(:,:) = 0.

C--   initialise tendency
        tdt3d = 0.
        qdt3d = 0.
        udt3d = 0.
        vdt3d = 0.
        cond_dt_tg = 0.
        cond_dt_qg = 0.

        IF (lwet_convection) THEN
          CALL DARGAN_BETTSMILLER(
     I                   deltaT, t3d, q3d, pFull3d, pHalf3d, coldT,
     O                   rain2d, snow2d, cond_dt_tg, cond_dt_qg,
     O                   q_ref,      bmflag,
     O                   klzbs,        cape,
     O                   cin,         t_ref,
     O                   invtau_bm_t,     invtau_bm_q,
     O                   capeflag,
     I                   bi,bj,myIter,myThid )

C-    store updated state to pass to LSC (addition from POG)
          t3d_tmp = t3d + cond_dt_tg
          q3d_tmp = q3d + cond_dt_qg

          cond_dt_tg = cond_dt_tg / deltaT
          cond_dt_qg = cond_dt_qg / deltaT
          rain2d = rain2d / deltaT
          EmPmR(1:sNx,1:sNy,bi,bj) = -rain2d(:,:)

          tdt3d = tdt3d + cond_dt_tg
          qdt3d = qdt3d + cond_dt_qg

#ifdef ALLOW_DIAGNOSTICS
          IF ( useDiagnostics ) THEN
           CALL DIAGNOSTICS_FILL( rain2d , 'AtPhCnvP',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( cape   , 'AtPhCAPE',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( cin    , 'AtPhCnIn',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( klzbs  , 'AtPhKlzb',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( bmflag , 'AtPhConv',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( invtau_bm_t, 'AtPhRlxT',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( invtau_bm_q, 'AtPhRlxQ',
     &                            0, 1, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( t_ref  , 'AtPh_Trf',
     &                           -1,Nr, 3, bi, bj, myThid )
           CALL DIAGNOSTICS_FILL( q_ref  , 'AtPh_Qrf',
     &                           -1,Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( cond_dt_tg, 'AtPhdTcv',
     &                           -1,Nr, 3, bi, bj, myThid )
c          snow2d = snow2d / deltaT
c          CALL DIAGNOSTICS_FILL( snow2d , 'SDIAG1  ',
c    &                            0, 1, 3, bi, bj, myThid )
          ENDIF
#endif /* ALLOW_DIAGNOSTICS */
        ELSE
          t3d_tmp = t3d
          q3d_tmp = q3d
        ENDIF

        cond_dt_tg = 0.
        cond_dt_qg = 0.
        rain2d(:,:) = 0.
        CALL LSCALE_COND(
     I              t3d_tmp, q3d_tmp, pFull3d, pHalf3d, coldT,
     O              rain2d, snow2d, cond_dt_tg, cond_dt_qg, q_ref,
     I              myThid )
        cond_dt_tg = cond_dt_tg / deltaT
        cond_dt_qg = cond_dt_qg / deltaT
        rain2d = rain2d / deltaT
        EmPmR(1:sNx,1:sNy,bi,bj) = EmPmR(1:sNx,1:sNy,bi,bj)
     &                           - rain2d(:,:)

        tdt3d = tdt3d + cond_dt_tg
        qdt3d = qdt3d + cond_dt_qg

#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
           CALL DIAGNOSTICS_FILL( rain2d , 'AtPhLscP',
     &                            0, 1, 3, bi, bj, myThid )
C     Re-use "q_ref" array to calculate Relative Humidity (in %)
           q_ref = 100. _d 0 * q3d_tmp/MAX( q_ref, 1. _d -12 )
           CALL DIAGNOSTICS_FILL( q_ref , 'RELHUM  ',
     &                           -1, Nr, 3, bi, bj, myThid )
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

        IF ( two_stream ) THEN
          CALL RADIATION_DOWN(
     I                   sNx,sNy, myTime, lat2d, pHalf3d, t3d, q3d,
     I                   albedo_2d, ozone_3d,
     O                   ozoneColmn, s_sw_dwn, s_lw_dwn,
     O                   dtrans_3d, dtrans_win, b_3d, b_win,
     O                   lw_down_3d, sw_down_3d, myThid )

c         write(errorMessageUnit,'(A,I5,I3,4F9.3)')
c    &      'it,bj,sw,lw:',myIter,bj,
c    &      sw_down_3d(1,1,15), sw_down_3d(1,sNy,15),
c    &      lw_down_3d(1,1,15), lw_down_3d(1,sNy,15)
        ENDIF

        IF (.TRUE.) THEN
         land = .false.
         avail = .true.
         rough_mom = roughness_mom
         rough_heat = roughness_heat
         rough_moist = roughness_moist
         gust = 1.0
         u_surf = 0.
         v_surf = 0.
         CALL SURFACE_FLUX(
     I            t3d(:,:,Nr), q3d(:,:,Nr), u3d(:,:,Nr), v3d(:,:,Nr),
     I            pFull3d(:,:,Nr), zFull3d(:,:,Nr), pHalf3d(:,:,Nr+1),
     I            t_surf, t_surf,
     U            q_surf,
     I            u_surf, v_surf,
     I            rough_mom, rough_heat, rough_moist, gust,
     O            flux_t, flux_q, flux_r, flux_u, flux_v,
     O            drag_m, drag_t, drag_q, w_atm,
     O            ustar, bstar, qstar,
     O            dhdt_surf, dedt_surf, dedq_surf, drdt_surf,
     O            dhdt_atm, dedq_atm, dtaudv_atm,
     I            deltaT, land(:,:), avail(:,:), myThid  )
        ENDIF

        IF ( two_stream ) THEN
          rad_dt_tg = tdt3d
          CALL RADIATION_UP(
     I                   sNx,sNy, myTime, lat2d, pHalf3d, t_surf, t3d,
     U                   tdt3d, lw_net_3d, sw_net_3d,
     I                   albedo_2d, ozoneColmn, dtrans_3d, dtrans_win,
     I                   b_3d, b_win, lw_down_3d, sw_down_3d, myThid )

c         write(errorMessageUnit,'(A,I5,I3,4F9.3)')
c    &      'it,bj,t,tdt:',myIter,bj,
c    &      t3d(1,1,15), t3d(1,sNy,15),
c    &      tdt3d(1,1,15)*86400., tdt3d(1,sNy,15)*86400.
          rad_dt_tg = tdt3d - rad_dt_tg
        ELSE
          rad_dt_tg = 0.
        ENDIF

        IF (turb) THEN
         fracland = 0.0
         rough = roughness_mom
         CALL VERT_TURB_DRIVER(  1, 1, myTime, myTime+deltaT, deltaT,
     I             fracland(:,:), pHalf3d, pFull3d, zHalf3d, zFull3d,
     I             ustar, bstar, rough,
     I             u3d, v3d, t3d, q3d,
     O             diff_t(:,:,:), diff_m(:,:,:), gust(:,:),
     I             myThid )
        ENDIF

        diff_dt_tg = tdt3d
        diff_dt_qg = qdt3d

        CALL GCM_VERT_DIFF_DOWN( 1, 1, deltaT,
     I           u3d, v3d, t3d, q3d,
     I           diff_m(:,:,:), diff_t(:,:,:),
     I           pHalf3d, pFull3d, zFull3d,
     U           flux_u(:,:),  flux_v(:,:),  dtaudv_atm,
     U           udt3d, vdt3d, tdt3d,
     I           qdt3d,
     O           diss_heat(:,:,:),
     U           tri_surf_dtmass,
     U           tri_surf_dflux_t, tri_surf_dflux_q,
     U           tri_surf_delta_t, tri_surf_delta_q,
     O           e_global, f_t_global, f_q_global,
     I           myThid )

        CALL MIXED_LAYER(
     I                    myTime,
     U                    t_surf(:,:),
     U                    flux_t(:,:), flux_q(:,:), flux_r(:,:),
     I                    deltaT,
     I                    s_sw_dwn(:,:), s_lw_dwn(:,:),
     U                    tri_surf_dtmass,
     U                    tri_surf_dflux_t, tri_surf_dflux_q,
     U                    tri_surf_delta_t, tri_surf_delta_q,
     I                    dhdt_surf(:,:), dedt_surf(:,:),
     I                    dedq_surf(:,:), drdt_surf(:,:),
     I                    dhdt_atm(:,:),  dedq_atm(:,:),
     I                    ocean_qflux, mixLayDepth,
     O                    delta_t_surf(:,:),
     I                    myThid )

        CALL GCM_VERT_DIFF_UP ( 1, 1, deltat,
     I           tri_surf_delta_t, tri_surf_delta_q,
     I           e_global, f_t_global, f_q_global,
     O           tdt3d, qdt3d,
     I           myThid )

        diff_dt_tg = tdt3d - diff_dt_tg
        diff_dt_qg = qdt3d - diff_dt_qg

C-- update Upward LW (assuming adjustment of surf LW emission goes directly to space)
        adj_lw_up = drdt_surf*delta_t_surf
        DO k=1,Nr+1
         lw_net_3d(:,:,k) = lw_net_3d(:,:,k) + adj_lw_up
        ENDDO

        DO k=1,Nr
          kc = Nr-k+1
          conv_T2theta = (atm_po/rC(kc))**atm_kappa
#ifdef NONLIN_FRSURF
          IF ( select_rStar.GE.1 ) THEN
           atmPhys_dT(1:sNx,1:sNy,kc,bi,bj) = tdt3d(:,:,k)
     &                   * conv_T2theta/pStarFacK(1:sNx,1:sNy,bi,bj)
          ELSE
#else /* NONLIN_FRSURF */
          IF ( .TRUE. ) THEN
#endif /* NONLIN_FRSURF */
           atmPhys_dT(1:sNx,1:sNy,kc,bi,bj) = tdt3d(:,:,k)
     &                   * conv_T2theta
          ENDIF
          atmPhys_dQ(1:sNx,1:sNy,kc,bi,bj)  = qdt3d(:,:,k)
C-    Note: multiply A-grid tend. by hFacC (<-> dpFac) here and
C     C-grid average will be divided by hFacW,S when applied to Dynamics
          dpFac(:,:) = ( pHalf3d(:,:,k+1) - pHalf3d(:,:,k)
     &                 )*recip_drF(kc)
          atmPhys_dU(1:sNx,1:sNy,kc,bi,bj) = udt3d(:,:,k)*dpFac(:,:)
          atmPhys_dV(1:sNx,1:sNy,kc,bi,bj) = vdt3d(:,:,k)*dpFac(:,:)
        ENDDO

C--   Update SST
        IF ( atmPhys_stepSST ) THEN
          atmPhys_SST(1:sNx,1:sNy,bi,bj) = t_surf(:,:)
        ENDIF

        EmPmR(1:sNx,1:sNy,bi,bj) = EmPmR(1:sNx,1:sNy,bi,bj)
     &                           + flux_q(:,:)
        Qnet(1:sNx,1:sNy,bi,bj) = flux_t(:,:) + flux_r(:,:)
     &                          - s_lw_dwn(:,:) - s_sw_dwn(:,:)
     &                          + flux_q(:,:)*HLv
        Qsw (1:sNx,1:sNy,bi,bj) = -s_sw_dwn(:,:)
#ifdef COMPONENT_MODULE
C-    surface wind-stress (+ = down) at grid cell center (A-grid):
        taux(1:sNx,1:sNy,bi,bj) = -flux_u
        tauy(1:sNx,1:sNy,bi,bj) = -flux_v
#endif /* COMPONENT_MODULE */

#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
          CALL DIAGNOSTICS_FILL( atmPhys_dT , 'AtPhdTdt',
     &                           0, Nr, 1, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( atmPhys_dQ , 'AtPhdQdt',
     &                           0, Nr, 1, bi, bj, myThid )
c         CALL DIAGNOSTICS_FILL( atmPhys_dU , 'AtPhdUdt',
c    &                           0, Nr, 1, bi, bj, myThid )
c         CALL DIAGNOSTICS_FILL( atmPhys_dV , 'AtPhdVdt',
c    &                           0, Nr, 1, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( udt3d , 'AtPhdUdt',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( vdt3d , 'AtPhdVdt',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( diff_t , 'AtPhDifT',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( diff_m , 'AtPhDifM',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( rad_dt_tg , 'AtPhdTrd',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( sw_net_3d , 'AtPhNSR ',
     &                           -1, Nr+1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( lw_net_3d , 'AtPhNLR ',
     &                           -1, Nr+1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( sw_down_3d, 'AtPhDSR ',
     &                           -1, Nr+1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( lw_down_3d, 'AtPhDLR ',
     &                           -1, Nr+1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( sw_down_3d, 'AtPhInSR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( sw_net_3d , 'AtPhNTSR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( lw_net_3d , 'AtPhOLR ',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( sw_down_3d(:,:,Nr+1),'AtPhDSSR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( s_sw_dwn,  'AtPhNSSR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( s_lw_dwn,  'AtPhDSLR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( flux_r ,   'AtPhUSLR',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( flux_t ,   'AtPhSens',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( flux_q ,   'AtPhEvap',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( flux_u , 'AtPhTauX',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( flux_v , 'AtPhTauY',
     &                           0, 1, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( cond_dt_tg, 'AtPhdTlc',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( diff_dt_tg, 'AtPhdtTg',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( diff_dt_qg, 'AtPhdtQg',
     &                           -1, Nr, 3, bi, bj, myThid )
          CALL DIAGNOSTICS_FILL( diss_heat,  'AtPhDisH',
     &                           -1, Nr, 3, bi, bj, myThid )
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C--   end bi,bj loops.
       ENDDO
      ENDDO

      CALL EXCH_UV_AGRID_3D_RL( atmPhys_dU, atmPhys_dV,
     &                          .TRUE., Nr, myThid )
#ifdef COMPONENT_MODULE
C-    surface wind-stress at grid cell center (A-grid):
      IF ( useCoupler ) THEN
        CALL EXCH_UV_AGRID_3D_RL( taux, tauy,
     &                           .TRUE., 1, myThid )
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=2-OLy,sNy+OLy
          DO i=2-OLx,sNx+OLx
            fu(i,j,bi,bj) = halfRL
     &                    *( taux(i-1,j,bi,bj) + taux(i,j,bi,bj) )
            fv(i,j,bi,bj) = halfRL
     &                    *( tauy(i,j-1,bi,bj) + tauy(i,j,bi,bj) )
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL ATM_STORE_MY_DATA( myTime, myIter, myThid )
      ENDIF
#endif /* COMPONENT_MODULE */

      RETURN
      END
