C $Header: /u/gcmpack/MITgcm/pkg/atm_phys/atm_phys_driver.F,v 1.10 2017/08/12 22:12:48 jmc Exp $
C $Name: $
#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 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-- 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/q_ref
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,
O 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, 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