diff --git a/physics/GFS_rrtmg_setup.F90 b/physics/GFS_rrtmg_setup.F90
index 384d5252d..80c033bf1 100644
--- a/physics/GFS_rrtmg_setup.F90
+++ b/physics/GFS_rrtmg_setup.F90
@@ -171,8 +171,6 @@ subroutine GFS_rrtmg_setup_init ( si, levr, ictm, isol, solar_file, ico2, &
       ! Initialize the CCPP error handling variables
       errmsg = ''
       errflg = 0
-
-      if (is_initialized) return
       
       if (do_RRTMGP) then
         write(errmsg,'(*(a))') "Logic error: do_RRTMGP must be set to .false."
@@ -213,6 +211,8 @@ subroutine GFS_rrtmg_setup_init ( si, levr, ictm, isol, solar_file, ico2, &
          print *, 'lextop=',lextop, ' ltp=',ltp
       endif
 
+      if (is_initialized) return
+
       ! Call initialization routines
       call sol_init ( me, isol, solar_file, con_solr_2008,con_solr_2002,&
            con_pi )
diff --git a/physics/GFS_rrtmgp_setup.F90 b/physics/GFS_rrtmgp_setup.F90
index 76db14279..3e90cc96b 100644
--- a/physics/GFS_rrtmgp_setup.F90
+++ b/physics/GFS_rrtmgp_setup.F90
@@ -74,8 +74,6 @@ subroutine GFS_rrtmgp_setup_init(do_RRTMGP, imp_physics, imp_physics_fer_hires,
     ! Initialize the CCPP error handling variables
     errmsg = ''
     errflg = 0
-    
-    if (is_initialized) return
 
     ! Consistency checks
     if (.not. do_RRTMGP) then
@@ -125,6 +123,8 @@ subroutine GFS_rrtmgp_setup_init(do_RRTMGP, imp_physics, imp_physics_fer_hires,
     iyear0 = 0
     monthd = 0
 
+    if (is_initialized) return
+
     ! Call initialization routines..
     call sol_init ( me, isol, solar_file, con_solr_2008, con_solr_2002, con_pi )
     call aer_init ( levr, me, iaermdl, iaerflg, lalw1bd, aeros_file, con_pi, con_t0c,    &
diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90
index c60247cf6..eecc5493c 100644
--- a/physics/module_sf_mynn.F90
+++ b/physics/module_sf_mynn.F90
@@ -106,6 +106,7 @@ MODULE module_sf_mynn
 
   REAL(kind_phys), DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, &
                                      psih_stab,psih_unstab
+!$acc declare create(psim_stab, psim_unstab, psih_stab, psih_unstab)
 
 CONTAINS
 
@@ -371,6 +372,20 @@ SUBROUTINE SFCLAY_mynn(                           &
       errflg = 0
       errmsg = ''
 
+!$acc enter data copyin( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D,       &
+!$acc                       pattern_spp_sfc)
+
+!$acc enter data copyin( UST_WAT(:), UST_LND(:),  UST_ICE(:),     &
+!$acc                    MOL(:),     QFLX(:),     HFLX(:),        &
+!$acc                    QSFC(:),    QSFC_WAT(:), QSFC_LND(:),    &
+!$acc                    QSFC_ICE(:))
+
+!$acc enter data create( dz8w1d(:),  dz2w1d(:),   U1D(:),         &
+!$acc                    V1D(:),     U1D2(:),     V1D2(:),        &
+!$acc                    QV1D(:),    QC1D(:),     P1D(:),         &
+!$acc                    T1D(:),     rstoch1D(:), qstar(:))
+
+
       IF (debug_code >= 1) THEN
         write(*,*)"======= printing of constants:"
         write(*,*)"cp=",    cp," g=",     grav
@@ -382,6 +397,10 @@ SUBROUTINE SFCLAY_mynn(                           &
       itf=ite !MIN0(ite,ide-1)
       ktf=kte !MIN0(kte,kde-1)
 
+!$acc parallel loop present(dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D,       &
+!$acc                       pattern_spp_sfc,dz8w1d,dz2w1d,U1D,    &
+!$acc                       V1D,U1D2,V1D2,QV1D,QC1D,P1D,T1D,      &
+!$acc                       rstoch1D,qstar)
       DO i=its,ite
          dz8w1d(I) = dz8w(i,kts)
          dz2w1d(I) = dz8w(i,kts+1)
@@ -403,6 +422,9 @@ SUBROUTINE SFCLAY_mynn(                           &
       ENDDO
 
       IF (itimestep==1 .AND. iter==1) THEN
+!$acc parallel loop present(U1D,V1D,UST_WAT,UST_LND,UST_ICE,MOL,  &
+!$acc                       QFLX,HFLX,QV3D,QSFC,QSFC_WAT,         &
+!$acc                       QSFC_LND,QSFC_ICE)
          DO i=its,ite
             IF (.not. flag_restart) THEN
                !Everything here is used before calculated
@@ -432,6 +454,9 @@ SUBROUTINE SFCLAY_mynn(                           &
          ENDDO
       ENDIF
 
+!$acc exit data delete( dz8w,U3D,V3D,QV3D,QC3D,P3D,T3D,         &
+!$acc                   pattern_spp_sfc, QC1D)
+
       CALL SFCLAY1D_mynn(flag_iter,                             &
            J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,                       &
            U1D2,V1D2,dz2w1d,                                    &
@@ -471,6 +496,16 @@ SUBROUTINE SFCLAY_mynn(                           &
            its,ite, jts,jte, kts,kte,                           &
            errmsg, errflg                                       )
 
+!$acc exit data copyout( UST_WAT(:), UST_LND(:),  UST_ICE(:),   &
+!$acc                    MOL(:),     QFLX(:),     HFLX(:),      &
+!$acc                    QSFC(:),    QSFC_WAT(:), QSFC_LND(:),  &
+!$acc                    QSFC_ICE(:))
+
+!$acc exit data delete( dz8w1d(:),   dz2w1d(:),   U1D(:),       &
+!$acc                   V1D(:),      U1D2(:),     V1D2(:),      &
+!$acc                   QV1D(:),     T1D(:),      P1D(:),       &
+!$acc                   rstoch1D(:), qstar(:))
+
     END SUBROUTINE SFCLAY_MYNN
 
 !-------------------------------------------------------------------
@@ -629,6 +664,22 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
       character(len=*), intent(inout) :: errmsg
       integer,          intent(inout) :: errflg
 
+! Local fixed-size errmsg character array for error messages on accelerator
+! devices distinct from the host (e.g. GPUs). Necessary since OpenACC does
+! not support assumed-size (len=*) arrays like errmsg. Additional
+! device_errflg integer to denote when device_errmsg needs to be synced
+! with errmsg.
+      character(len=512) :: device_errmsg
+      integer            :: device_errflg
+
+! Special versions of the fixed-size errmsg character array for error messages
+! on the device and it's errflag counterpart. These are necessary to ensure
+! the return statements at lines 1417 and 2030 are executed only for this
+! special case, and not any and all error messages set on the device.
+      character(len=512) :: device_special_errmsg
+      integer            :: device_special_errflg
+
+
 !----------------------------------------------------------------
 ! LOCAL VARS
 !----------------------------------------------------------------
@@ -678,7 +729,65 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
       ! Initialize error-handling
       errflg = 0
       errmsg = ''
+      device_errflg = errflg
+      device_errmsg = errmsg
+      device_special_errflg = errflg
+      device_special_errmsg = errmsg
 !-------------------------------------------------------------------
+!$acc update device(psim_stab, psim_unstab, psih_stab, psih_unstab)
+
+!$acc enter data create( ZA,     ZA2,    THV1D,  TH1D,   TC1D,   TV1D,  &
+!$acc                    RHO1D,  QVSH,   PSIH2,  PSIM10, PSIH10, WSPDI, &
+!$acc                    GOVRTH, PSFC,   THCON,                         &
+!$acc                    zratio_lnd,   zratio_ice,   zratio_wat,        &
+!$acc                    TSK_lnd,      TSK_ice,      TSK_wat,           &
+!$acc                    THSK_lnd,     THSK_ice,     THSK_wat,          &
+!$acc                    THVSK_lnd,    THVSK_ice,    THVSK_wat,         &
+!$acc                    GZ1OZ0_lnd,   GZ1OZ0_ice,   GZ1OZ0_wat,        &
+!$acc                    GZ1OZt_lnd,   GZ1OZt_ice,   GZ1OZt_wat,        &
+!$acc                    GZ2OZ0_lnd,   GZ2OZ0_ice,   GZ2OZ0_wat,        &
+!$acc                    GZ2OZt_lnd,   GZ2OZt_ice,   GZ2OZt_wat,        &
+!$acc                    GZ10OZ0_lnd,  GZ10OZ0_ice,  GZ10OZ0_wat,       &
+!$acc                    GZ10OZt_lnd,  GZ10OZt_ice,  GZ10OZt_wat,       &
+!$acc                    ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat,      &
+!$acc                    ZT_lnd,       ZT_ice,       ZT_wat,            &
+!$acc                    ZQ_lnd,       ZQ_ice,       ZQ_wat,            &
+!$acc                    PSIQ_lnd,     PSIQ_ice,     PSIQ_wat,          &
+!$acc                    PSIQ2_lnd,    PSIQ2_ice,    PSIQ2_wat,         &
+!$acc                    QSFCMR_lnd,   QSFCMR_ice,   QSFCMR_wat )
+
+!$acc enter data copyin(flag_iter, dry, wet, icy, CPM, MAVAIL, &
+!$acc                   QFX, FLHC, FLQC, CHS, CH, CHS2, CQS2, USTM, &
+!$acc                   HFX, LH, wstar, qstar, PBLH, ZOL, MOL, RMOL, &
+!$acc                   T2, TH2, Q2, QV1D, PSFCPA, &
+!$acc                   WSPD, U10, V10, U1D, V1D, U1D2, V1D2, &
+!$acc                   T1D, P1D, rstoch1D, sigmaf, &
+!$acc                   shdmax, vegtype, z0pert, ztpert, dx, QGH, &
+!$acc                   dz2w1d, dz8w1d, &
+!$acc                   stress_wat, stress_lnd, stress_ice, &
+!$acc                   rb_wat, rb_lnd, rb_ice, &
+!$acc                   tskin_wat, tskin_lnd, tskin_ice, &
+!$acc                   tsurf_wat, tsurf_lnd, tsurf_ice, &
+!$acc                   psim, psih, &
+!$acc                   UST_wat, UST_lnd, UST_ice, &
+!$acc                   ZNT_wat, ZNT_lnd, ZNT_ice, &
+!$acc                   QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
+!$acc                   QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
+!$acc                   HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
+!$acc                   PSIX_wat, PSIX_lnd, PSIX_ice, &
+!$acc                   PSIX10_wat, PSIX10_lnd, PSIX10_ice, &
+!$acc                   PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
+!$acc                   PSIT_lnd, PSIT_wat, PSIT_ice, &
+!$acc                   ch_lnd, ch_wat, ch_ice, &
+!$acc                   cm_lnd, cm_wat, cm_ice, &
+!$acc                   snowh_lnd, snowh_wat, snowh_ice, &
+!$acc                   device_errmsg, device_errflg, &
+!$acc                   device_special_errmsg, device_special_errflg)
+
+!$acc parallel loop present(PSFCPA, PSFC, QSFC, T1D, flag_iter, tsurf_lnd,    &
+!$acc                       QSFC_wat, QSFCMR_wat, wet, TSK_wat, tskin_wat,    &
+!$acc                       QSFC_lnd, QSFCMR_lnd, dry, TSK_lnd, tskin_lnd,    &
+!$acc                       QSFC_ice, QSFCMR_ice, icy, TSK_ice, tskin_ice)
       DO I=its,ite
 
          ! PSFC ( in cmb) is used later in saturation checks
@@ -791,6 +900,10 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
        endif ! flag_iter
       ENDDO
 
+!$acc serial present(pblh, PSFCPA, dz8w1d, qflx, hflx,                       &
+!$acc      dry, tskin_lnd, tsurf_lnd, qsfc_lnd, znt_lnd, ust_lnd, snowh_lnd, &
+!$acc      icy, tskin_ice, tsurf_ice, qsfc_ice, znt_ice, ust_ice, snowh_ice, &
+!$acc      wet, tskin_wat, tsurf_wat, qsfc_wat, znt_wat, ust_wat, snowh_wat)
       IF (debug_code >= 1) THEN
         write(0,*)"ITIMESTEP=",ITIMESTEP," iter=",iter
         DO I=its,ite
@@ -815,7 +928,12 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
            ENDIF
         ENDDO
       ENDIF
+!$acc end serial
 
+!$acc parallel loop present(PSFC, PSFCPA, QVSH, QV1D, THCON, flag_iter,        &
+!$acc       dry, tskin_lnd, TSK_lnd, tsurf_lnd, THSK_lnd, THVSK_lnd, qsfc_lnd, &
+!$acc       icy, tskin_ice, TSK_ice, tsurf_ice, THSK_ice, THVSK_ice, qsfc_ice, &
+!$acc       wet, tskin_wat, TSK_wat, tsurf_wat, THSK_wat, THVSK_wat, qsfc_wat)
       DO I=its,ite
          ! PSFC ( in cmb) is used later in saturation checks
          PSFC(I)=PSFCPA(I)/1000.
@@ -853,18 +971,21 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
         endif ! flag_iter
       ENDDO
 
+!$acc parallel loop present(TH1D, T1D, P1D, TC1D)
       DO I=its,ite
          ! CONVERT LOWEST LAYER TEMPERATURE TO POTENTIAL TEMPERATURE:
          TH1D(I)=T1D(I)*(100000./P1D(I))**ROVCP  !(Theta, K)
          TC1D(I)=T1D(I)-273.15                   !(T, Celsius)
       ENDDO
 
+!$acc parallel loop present(THV1D, TH1D, QVSH, TV1D, T1D)
       DO I=its,ite
          ! CONVERT TO VIRTUAL TEMPERATURE
          THV1D(I)=TH1D(I)*(1.+EP1*QVSH(I))             !(K)
          TV1D(I)=T1D(I)*(1.+EP1*QVSH(I))               !(K)
       ENDDO
 
+!$acc parallel loop present(RHO1D, P1D, TV1D, TH1D, ZA, ZA2, dz2w1d, dz8w1d, GOVRTH)
       DO I=its,ite
          RHO1D(I)=P1D(I)/(Rd*TV1D(I))     !now using value calculated in sfc driver
          ZA(I)=0.5*dz8w1d(I)              !height of first half-sigma level
@@ -873,11 +994,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
       ENDDO
 
       !tgs - should QFX and HFX be separate for land, ice and water?
+!$acc parallel loop present(QFX, QFLX, RHO1D, HFX, HFLX)
       DO I=its,ite
          QFX(i)=QFLX(i)*RHO1D(I)
          HFX(i)=HFLX(i)*RHO1D(I)*cp
       ENDDO
 
+!$acc serial present(THV1D, TV1D, RHO1D, GOVRTH, &
+!$acc                dry, tsk_lnd, thvsk_lnd,    &
+!$acc                icy, tsk_ice, thvsk_ice,    &
+!$acc                wet, tsk_wat, thvsk_wat)
       IF (debug_code ==2) THEN
         !write(*,*)"ITIMESTEP=",ITIMESTEP
         DO I=its,ite
@@ -890,7 +1016,9 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
           write(*,*)"RHO1D=", RHO1D(i)," GOVRTH=",GOVRTH(i)
         ENDDO
       ENDIF
+!$acc end serial
 
+!$acc parallel loop present(T1D,P1D,QGH,QV1D,CPM)
       DO I=its,ite
          ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP
          ! Q2SAT = QGH IN LSM
@@ -908,6 +1036,10 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
          CPM(I)=CP*(1.+0.84*QV1D(I))
       ENDDO
 
+!$acc serial present(QGH,                       &
+!$acc                wet, QSFC_wat, QSFCMR_wat, &
+!$acc                dry, QSFC_lnd, QSFCMR_lnd, &
+!$acc                icy, QSFC_ice, QSFCMR_ice)
       IF (debug_code == 2) THEN
          write(*,*)"ITIMESTEP=",ITIMESTEP
          DO I=its,ite
@@ -925,7 +1057,13 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
             endif
          ENDDO
       ENDIF
+!$acc end serial
 
+!$acc parallel loop present(flag_iter,U1D,V1D,WSPD,wet,dry,icy,    &
+!$acc                    THV1D,THVSK_wat,THVSK_lnd,THVSK_ice,   &
+!$acc                    hfx,RHO1D,qfx,WSTAR,pblh,dx,GOVRTH,ZA, &
+!$acc                    TSK_wat,TSK_lnd,TSK_ice,               &
+!$acc                    rb_wat,rb_lnd,rb_ice)
       DO I=its,ite
         if( flag_iter(i) ) then
          ! DH* 20200401 - note. A weird bug in Intel 18 on hera prevents using the
@@ -1067,6 +1205,35 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
 !--------------------------------------------------------------------
 !--------------------------------------------------------------------
 
+!$acc parallel loop present(flag_iter, PSFCPA, dz8w1d, pblh, &
+!$acc                    device_errmsg, device_errflg, &
+!$acc                    device_special_errmsg, device_special_errflg, &
+!$acc                    wet, dry, icy, &
+!$acc                    ZT_wat, ZT_lnd, ZT_ice, &
+!$acc                    ZNT_wat, ZNT_lnd, ZNT_ice, &
+!$acc                    ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc                    UST_wat, UST_lnd, UST_ice, &
+!$acc                    ZQ_wat, ZQ_lnd, ZQ_ice, &
+!$acc                    snowh_wat, snowh_lnd, snowh_ice, &
+!$acc                    THVSK_wat, THVSK_lnd, THVSK_ice, &
+!$acc                    tskin_wat, tskin_lnd, tskin_ice, &
+!$acc                    tsurf_wat, tsurf_lnd, tsurf_ice, &
+!$acc                    qsfc_wat, qsfc_lnd, qsfc_ice, &
+!$acc                    GZ1OZ0_wat, GZ1OZt_wat, GZ2OZ0_wat, GZ2OZt_wat, GZ10OZ0_wat, GZ10OZt_wat, &
+!$acc                    GZ1OZ0_lnd, GZ1OZt_lnd, GZ2OZ0_lnd, GZ2OZt_lnd, GZ10OZ0_lnd, GZ10OZt_lnd, &
+!$acc                    GZ1OZ0_ice, GZ1OZt_ice, GZ2OZ0_ice, GZ2OZt_ice, GZ10OZ0_ice, GZ10OZt_ice, &
+!$acc                    zratio_wat, zratio_lnd, zratio_ice, &
+!$acc                    stress_wat, stress_lnd, stress_ice, &
+!$acc                    rb_wat, rb_lnd, rb_ice, &
+!$acc                    qflx, qflx_lnd, &
+!$acc                    hflx, hflx_lnd, &
+!$acc                    psim, psih, psim10, psih10, psih2, &
+!$acc                    psix_wat, psix10_wat, psit_wat, psit2_wat, psiq_wat, psiq2_wat, &
+!$acc                    psix_lnd, psix10_lnd, psit_lnd, psit2_lnd, psiq_lnd, psiq2_lnd, &
+!$acc                    psix_ice, psix10_ice, psit_ice, psit2_ice, psiq_ice, psiq2_ice, &
+!$acc                    WSPD, WSPDI, U1D, V1D, TC1D, THV1D, rstoch1D, USTM, ZA, ZOL, QVSH, &
+!$acc                    shdmax, vegtype, z0pert, ztpert, mol, rmol, wstar, qstar, sigmaf)
+
  DO I=its,ite
    if( flag_iter(i) ) then
 
@@ -1082,10 +1249,12 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
        if (sfc_z0_type >= 0) then ! Avoid calculation is using wave model
           ! CALCULATE z0 (znt)
           !--------------------------------------
+
           IF (debug_code == 2) THEN
             write(*,*)"=============Input to ZNT over water:"
             write(*,*)"u*:",UST_wat(i)," wspd=",WSPD(i)," visc=",visc," za=",ZA(I)
           ENDIF
+
           IF ( PRESENT(ISFTCFLX) ) THEN
              IF ( ISFTCFLX .EQ. 0 ) THEN
                 IF (COARE_OPT .EQ. 3.0) THEN
@@ -1170,7 +1339,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
              ENDIF
           ELSEIF ( ISFTCFLX .EQ. 4 ) THEN
              !GFS zt formulation
-             CALL GFS_zt_wat(ZT_wat(i),ZNTstoch_wat(i),restar,WSPD(i),ZA(i),sfc_z0_type,errmsg,errflg)
+             CALL GFS_zt_wat(ZT_wat(i),ZNTstoch_wat(i),restar,WSPD(i),ZA(i),sfc_z0_type,device_errmsg,device_errflg)
              ZQ_wat(i)=ZT_wat(i)
           ENDIF
        ELSE
@@ -1183,6 +1352,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
                                     rstoch1D(i),spp_sfc)
           ENDIF
        ENDIF
+
        IF (debug_code > 1) THEN
          write(*,*)"=============Output ZT & ZQ over water:"
          write(*,*)"ZT:",ZT_wat(i)," ZQ:",ZQ_wat(i)
@@ -1230,9 +1400,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
              ELSEIF ( IZ0TLND .EQ. 2 ) THEN
                 ! DH note - at this point, qstar is either not initialized
                 ! or initialized to zero, but certainly not set correctly
-                errmsg = 'Logic error: qstar is not set correctly when calling Yang_2008'
-                errflg = 1
+                device_special_errmsg = 'Logic error: qstar is not set correctly when calling Yang_2008'
+                device_special_errflg = 1
+#ifndef _OPENACC
+! Necessary since OpenACC does not support branching in parallel code
+! Must sync errmsg and errflg with device_errmsg and device_errflg, respectively
+! so that proper error message and error flag codes are returned.
+                errmsg = device_special_errmsg
+                errflg = device_special_errflg
                 return
+#endif
                 CALL Yang_2008(ZNTSTOCH_lnd(i),ZT_lnd(i),ZQ_lnd(i),UST_lnd(i),MOL(I),&
                               qstar(I),restar,visc)
              ELSEIF ( IZ0TLND .EQ. 3 ) THEN
@@ -1249,6 +1426,7 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
                          UST_lnd(I),KARMAN,1.0_kind_phys,0,spp_sfc,rstoch1D(i))
           ENDIF
        ENDIF
+
        IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN
          write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i
          write(0,*)" ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i)
@@ -1258,7 +1436,6 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
          " dz=",dz8w1d(i)," qflx=",qflx_lnd(i)," hflx=",hflx_lnd(i)," hpbl=",pblh(i)
        ENDIF
 
-
        GZ1OZ0_lnd(I)= LOG((ZA(I)+ZNTstoch_lnd(I))/ZNTstoch_lnd(I))
        GZ1OZt_lnd(I)= LOG((ZA(I)+ZNTstoch_lnd(i))/ZT_lnd(i))
        GZ2OZ0_lnd(I)= LOG((2.0+ZNTstoch_lnd(I))/ZNTstoch_lnd(I))
@@ -1821,6 +1998,26 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
    endif ! flag_iter
  ENDDO   ! end i-loop
 
+#ifdef _OPENACC
+! Necessary since OpenACC does not support branching in parallel code.
+! Must sync host errflg, errmsg to determine if return must be triggered
+! and correct error message and error flag code returned.
+! This code is being executed on the HOST side only, pulling data from DEVICE.
+!$acc exit data copyout(device_special_errflg, device_special_errmsg)
+ IF (device_special_errflg /= 0) THEN
+    errflg = device_special_errflg
+    errmsg = device_special_errmsg
+    return
+ ENDIF
+#endif
+
+!$acc serial present(wet, dry, icy,                  &
+!$acc        PSIM, PSIH, CPM, RHO1D, ZOL, wspd, MOL, &
+!$acc        wstar, qstar, THV1D, HFX, MAVAIL, QVSH, &
+!$acc         THVSK_wat,    THVSK_lnd,    THVSK_ice, &
+!$acc           UST_wat,      UST_lnd,      UST_ice, &
+!$acc      ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc            zt_wat,       zt_lnd,       zt_ice)
  IF (debug_code == 2) THEN
     DO I=its,ite
        IF(wet(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(wet)"
@@ -1841,10 +2038,29 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
        write(*,*)"============================================="
     ENDDO ! end i-loop
  ENDIF
+!$acc end serial
 
    !----------------------------------------------------------
    !  COMPUTE SURFACE HEAT AND MOISTURE FLUXES
    !----------------------------------------------------------
+!$acc parallel loop present(flag_iter, dry, wet, icy, &
+!$acc                    QFX, HFX, FLHC, FLQC, LH, CHS, CH, CHS2, CQS2, &
+!$acc                    RHO1D, MAVAIL, USTM, &
+!$acc                    UST_lnd, UST_wat, UST_ice, &
+!$acc                    PSIQ_lnd, PSIT_lnd, PSIX_lnd, &
+!$acc                    PSIQ_wat, PSIT_wat, PSIX_wat, &
+!$acc                    PSIQ_ice, PSIT_ice, PSIX_ice, &
+!$acc                    PSIQ2_lnd, PSIT2_lnd, &
+!$acc                    PSIQ2_wat, PSIT2_wat, &
+!$acc                    PSIQ2_ice, PSIT2_ice, &
+!$acc                    QSFC, QSFC_lnd, QSFC_wat, QSFC_ice, &
+!$acc                    QFLX, QFLX_lnd, QFLX_wat, QFLX_ice, &
+!$acc                    HFLX, HFLX_lnd, HFLX_wat, HFLX_ice, &
+!$acc                    QSFCMR_lnd, QSFCMR_wat, QSFCMR_ice, &
+!$acc                    QV1D, WSPD, WSPDI, CPM, TH1D, &
+!$acc                    THSK_lnd, THSK_wat, THSK_ice, &
+!$acc                    ch_lnd, ch_wat, ch_ice, &
+!$acc                    cm_lnd, cm_wat, cm_ice)
  DO I=its,ite
   if( flag_iter(i) ) then
 
@@ -2040,6 +2256,18 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
 ENDDO ! end i-loop
 
 IF (compute_diag) then
+   !$acc parallel loop present(flag_iter, dry, wet, icy, &
+   !$acc                    ZA, ZA2, T2, TH2, TH1D, Q2, QV1D, PSFCPA, &
+   !$acc                    THSK_lnd, THSK_wat, THSK_ice, &
+   !$acc                    QSFC_lnd, QSFC_wat, QSFC_ice, &
+   !$acc                    U10, V10, U1D, V1D, U1D2, V1D2, &
+   !$acc                    ZNTstoch_lnd, ZNTstoch_lnd, ZNTstoch_ice, &
+   !$acc                    PSIX_lnd, PSIX_wat, PSIX_ice, &
+   !$acc                    PSIX10_lnd, PSIX10_wat, PSIX10_ice, &
+   !$acc                    PSIT2_lnd, PSIT2_wat, PSIT2_ice, &
+   !$acc                    PSIT_lnd, PSIT_wat, PSIT_ice, &
+   !$acc                    PSIQ2_lnd, PSIQ2_wat, PSIQ2_ice, &
+   !$acc                    PSIQ_lnd, PSIQ_wat, PSIQ_ice)
    DO I=its,ite
      if( flag_iter(i) ) then
       !-----------------------------------------------------
@@ -2153,6 +2381,16 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
 !-----------------------------------------------------
 ! DEBUG - SUSPICIOUS VALUES
 !-----------------------------------------------------
+!$acc serial present(dry, wet, icy, CPM, MAVAIL, &
+!$acc       HFX, LH, wstar, RHO1D, PBLH, ZOL, ZA, MOL, &
+!$acc       PSIM, PSIH, WSTAR, T1D, TH1D, THV1D, QVSH, &
+!$acc       UST_wat, UST_lnd, UST_ice, &
+!$acc       THSK_wat, THSK_lnd, THSK_ice, &
+!$acc       THVSK_wat, THVSK_lnd, THVSK_ice, &
+!$acc       ZNTstoch_wat, ZNTstoch_lnd, ZNTstoch_ice, &
+!$acc       ZT_wat, ZT_lnd, ZT_ice, &
+!$acc       QSFC_wat, QSFC_lnd, QSFC_ice, &
+!$acc       PSIX_wat, PSIX_lnd, PSIX_ice)
 IF ( debug_code == 2) THEN
    DO I=its,ite
       yesno = 0
@@ -2257,6 +2495,62 @@ SUBROUTINE SFCLAY1D_mynn(flag_iter,                            &
       ENDIF
     ENDDO ! end i-loop
  ENDIF ! end debug option
+!$acc end serial
+
+!$acc exit data copyout(CPM, FLHC, FLQC, CHS, CH, CHS2, CQS2,&
+!$acc                   USTM, wstar, qstar, ZOL, MOL, RMOL,  &
+!$acc                   HFX, QFX, LH, QSFC, QFLX, HFLX,      &
+!$acc                   T2, TH2, Q2, WSPD, U10, V10,         &
+!$acc                   QGH, psim, psih,                     &
+!$acc                   stress_wat, stress_lnd, stress_ice,  &
+!$acc                   rb_wat,     rb_lnd,     rb_ice,      &
+!$acc                   UST_wat,    UST_lnd,    UST_ice,     &
+!$acc                   ZNT_wat,    ZNT_lnd,    ZNT_ice,     &
+!$acc                   QSFC_lnd,   QSFC_wat,   QSFC_ice,    &
+!$acc                   QFLX_lnd,   QFLX_wat,   QFLX_ice,    &
+!$acc                   HFLX_lnd,   HFLX_wat,   HFLX_ice,    &
+!$acc                   PSIX_wat,   PSIX_lnd,   PSIX_ice,    &
+!$acc                   PSIX10_wat, PSIX10_lnd, PSIX10_ice,  &
+!$acc                   PSIT2_lnd,  PSIT2_wat,  PSIT2_ice,   &
+!$acc                   PSIT_lnd,   PSIT_wat,   PSIT_ice,    &
+!$acc                   ch_lnd,     ch_wat,     ch_ice,      &
+!$acc                   cm_lnd,     cm_wat,     cm_ice,      &
+!$acc                   device_errmsg, device_errflg)
+
+! Final sync of device and host error flags and messages
+IF (device_errflg /= 0) THEN
+    errflg = device_errflg
+    errmsg = device_errmsg
+ENDIF
+
+!$acc exit data delete( flag_iter, dry, wet, icy, dx,         &
+!$acc                   MAVAIL, PBLH, PSFCPA, z0pert, ztpert, &
+!$acc                   QV1D, U1D, V1D, U1D2, V1D2, T1D, P1D, &
+!$acc                   rstoch1D, sigmaf, shdmax, vegtype,    &
+!$acc                   dz2w1d, dz8w1d,                       &
+!$acc                   snowh_wat, snowh_lnd, snowh_ice,      &
+!$acc                   tskin_wat, tskin_lnd, tskin_ice,      &
+!$acc                   tsurf_wat, tsurf_lnd, tsurf_ice)
+
+!$acc exit data delete(  ZA,     ZA2,    THV1D,  TH1D,   TC1D,   TV1D,  &
+!$acc                    RHO1D,  QVSH,   PSIH2,  PSIM10, PSIH10, WSPDI, &
+!$acc                    GOVRTH, PSFC,   THCON,                         &
+!$acc                    zratio_lnd,   zratio_ice,   zratio_wat,        &
+!$acc                    TSK_lnd,      TSK_ice,      TSK_wat,           &
+!$acc                    THSK_lnd,     THSK_ice,     THSK_wat,          &
+!$acc                    THVSK_lnd,    THVSK_ice,    THVSK_wat,         &
+!$acc                    GZ1OZ0_lnd,   GZ1OZ0_ice,   GZ1OZ0_wat,        &
+!$acc                    GZ1OZt_lnd,   GZ1OZt_ice,   GZ1OZt_wat,        &
+!$acc                    GZ2OZ0_lnd,   GZ2OZ0_ice,   GZ2OZ0_wat,        &
+!$acc                    GZ2OZt_lnd,   GZ2OZt_ice,   GZ2OZt_wat,        &
+!$acc                    GZ10OZ0_lnd,  GZ10OZ0_ice,  GZ10OZ0_wat,       &
+!$acc                    GZ10OZt_lnd,  GZ10OZt_ice,  GZ10OZt_wat,       &
+!$acc                    ZNTstoch_lnd, ZNTstoch_ice, ZNTstoch_wat,      &
+!$acc                    ZT_lnd,       ZT_ice,       ZT_wat,            &
+!$acc                    ZQ_lnd,       ZQ_ice,       ZQ_wat,            &
+!$acc                    PSIQ_lnd,     PSIQ_ice,     PSIQ_wat,          &
+!$acc                    PSIQ2_lnd,    PSIQ2_ice,    PSIQ2_wat,         &
+!$acc                    QSFCMR_lnd,   QSFCMR_ice,   QSFCMR_wat )
 
 END SUBROUTINE SFCLAY1D_mynn
 !-------------------------------------------------------------------
@@ -2272,6 +2566,7 @@ END SUBROUTINE SFCLAY1D_mynn
   SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,&
         & landsea,IZ0TLND2,spp_sfc,rstoch)
 
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)       :: Z_0,restar,ustar,KARMAN,landsea
        INTEGER,  OPTIONAL,   INTENT(IN)  :: IZ0TLND2
@@ -2341,6 +2636,7 @@ SUBROUTINE davis_etal_2008(Z_0,ustar)
     !This is an update version from Davis et al. 2008, which
     !corrects a small-bias in Z_0 (AHW real-time 2012).
 
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: ustar
        REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2368,7 +2664,7 @@ END SUBROUTINE davis_etal_2008
 !>This formulation for roughness length was designed account for.
 !!wave steepness.
    SUBROUTINE Taylor_Yelland_2001(Z_0,ustar,wsp10)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: ustar,wsp10
        REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2396,7 +2692,7 @@ END SUBROUTINE Taylor_Yelland_2001
 !! The Charnock parameter CZC is varied from .011 to .018.
 !! between 10-m wsp = 10 and 18..
    SUBROUTINE charnock_1955(Z_0,ustar,wsp10,visc,zu)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: ustar, visc, wsp10, zu
        REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2421,7 +2717,7 @@ END SUBROUTINE charnock_1955
 !!The Charnock parameter CZC is varied from about .005 to .028
 !!between 10-m wind speeds of 6 and 19 m/s.
    SUBROUTINE edson_etal_2013(Z_0,ustar,wsp10,visc,zu)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: ustar, visc, wsp10, zu
        REAL(kind_phys), INTENT(OUT) :: Z_0
@@ -2450,7 +2746,7 @@ END SUBROUTINE edson_etal_2013
 !!data. The formula for land uses a constant ratio (Z_0/7.4) taken
 !!from Garratt (1992).
    SUBROUTINE garratt_1992(Zt,Zq,Z_0,Ren,landsea)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: Ren, Z_0,landsea
        REAL(kind_phys), INTENT(OUT) :: Zt,Zq
@@ -2486,7 +2782,7 @@ END SUBROUTINE garratt_1992
 !!
 !!This is for use over water only.
     SUBROUTINE fairall_etal_2003(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)   :: Ren,ustar,visc,rstoch
        INTEGER, INTENT(IN)           :: spp_sfc
@@ -2530,7 +2826,7 @@ END SUBROUTINE fairall_etal_2003
 !! The actual reference is unknown. This was passed along by Jim Edson (personal communication).
 !! This is for use over water only, preferably open ocean.
     SUBROUTINE fairall_etal_2014(Zt,Zq,Ren,ustar,visc,rstoch,spp_sfc)
-
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: Ren,ustar,visc,rstoch
        INTEGER, INTENT(IN)          :: spp_sfc
@@ -2578,6 +2874,7 @@ END SUBROUTINE fairall_etal_2014
 !!This should only be used over land!
        SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,qst,Ren,visc)
 
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: Z_0, Ren, ustar, tstar, qst, visc
        REAL(kind_phys) ::      ht,     &! roughness height at critical Reynolds number
@@ -2613,6 +2910,7 @@ END SUBROUTINE Yang_2008
 !>\ingroup mynn_sfc
     SUBROUTINE GFS_z0_lnd(z0max,shdmax,z1,vegtype,ivegsrc,z0pert)
 
+        !$acc routine seq
         REAL(kind_phys), INTENT(OUT)  :: z0max
         REAL(kind_phys), INTENT(IN)   :: shdmax,z1,z0pert
         INTEGER, INTENT(IN)           :: vegtype,ivegsrc
@@ -2673,6 +2971,7 @@ END SUBROUTINE GFS_z0_lnd
 !>\ingroup mynn_sfc
     SUBROUTINE GFS_zt_lnd(ztmax,z0max,sigmaf,ztpert,ustar_lnd)
 
+        !$acc routine seq
         REAL(kind_phys), INTENT(OUT)  :: ztmax
         REAL(kind_phys), INTENT(IN)   :: z0max,sigmaf,ztpert,ustar_lnd
         REAL(kind_phys)               :: czilc, tem1, tem2
@@ -2701,6 +3000,7 @@ END SUBROUTINE GFS_zt_lnd
 !>\ingroup mynn_sfc
     SUBROUTINE GFS_z0_wat(z0rl_wat,ustar_wat,WSPD,z1,sfc_z0_type,redrag)
 
+        !$acc routine seq
         REAL(kind_phys), INTENT(OUT)  :: z0rl_wat
         REAL(kind_phys), INTENT(INOUT):: ustar_wat
         REAL(kind_phys), INTENT(IN)   :: wspd,z1
@@ -2752,19 +3052,27 @@ SUBROUTINE GFS_z0_wat(z0rl_wat,ustar_wat,WSPD,z1,sfc_z0_type,redrag)
     END SUBROUTINE GFS_z0_wat
 !--------------------------------------------------------------------
 !>\ingroup mynn_sfc
-    SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,errmsg,errflg)
-
+    SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,device_errmsg,device_errflg)
+        !$acc routine seq
         real(kind_phys), INTENT(OUT)  :: ztmax
         real(kind_phys), INTENT(IN)   :: wspd,z1,z0rl_wat,restar
         INTEGER, INTENT(IN)           :: sfc_z0_type
-        character(len=*), intent(out) :: errmsg
-        integer,          intent(out) :: errflg
+
+! Using device_errmsg and device_errflg rather than the CCPP errmsg and errflg
+! so that this subroutine can be run on an accelerator device with OpenACC.
+!        character(len=*), intent(out) :: errmsg
+!        integer,          intent(out) :: errflg
+        character(len=512), intent(out) :: device_errmsg
+        integer,            intent(out) :: device_errflg
+
         real(kind_phys)               :: z0,z0max,wind10m,rat,ustar_wat
         real(kind_phys), PARAMETER    :: charnock = 0.014, z0s_max=.317e-2
 
         ! Initialize error-handling
-        errflg = 0
-        errmsg = ''
+!        errflg = 0
+!        errmsg = ''
+        device_errflg = 0
+        device_errmsg = ''
 
 !            z0           = 0.01 * z0rl_wat
 !Already converted to meters in the wrapper
@@ -2795,9 +3103,12 @@ SUBROUTINE GFS_zt_wat(ztmax,z0rl_wat,restar,WSPD,z1,sfc_z0_type,errmsg,errflg)
               call znot_t_v7(wind10m, ztmax)   ! 10-m wind,m/s, ztmax(m)
             else if (sfc_z0_type > 0) then
               write(0,*)'no option for sfc_z0_type=',sfc_z0_type
-              errflg = 1
-              errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
+!              errflg = 1
+!              errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
+              device_errflg = 1
+              device_errmsg = 'ERROR(GFS_zt_wat): sfc_z0_type not valid.'
               return
+
             endif
 
     END SUBROUTINE GFS_zt_wat
@@ -2807,6 +3118,7 @@ END SUBROUTINE GFS_zt_wat
 !! Weiguo Wang, 2019-0425
 
       SUBROUTINE znot_m_v6(uref, znotm)
+      !$acc routine seq
       use machine , only : kind_phys
       IMPLICIT NONE
 ! Calculate areodynamical roughness over water with input 10-m wind
@@ -2856,6 +3168,7 @@ END SUBROUTINE znot_m_v6
 !!
       SUBROUTINE znot_t_v6(uref, znott)
 
+      !$acc routine seq
       IMPLICIT NONE
 !> Calculate scalar roughness over water with input 10-m wind
 !! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
@@ -2922,6 +3235,7 @@ END SUBROUTINE znot_t_v6
 !!
       SUBROUTINE znot_m_v7(uref, znotm)
 
+      !$acc routine seq
       IMPLICIT NONE
 !> Calculate areodynamical roughness over water with input 10-m wind
 !! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
@@ -2971,6 +3285,7 @@ END SUBROUTINE znot_m_v7
 !!
       SUBROUTINE znot_t_v7(uref, znott)
 
+      !$acc routine seq
       IMPLICIT NONE
 !> Calculate scalar roughness over water with input 10-m wind
 !! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
@@ -3040,6 +3355,7 @@ END SUBROUTINE znot_t_v7
 !! This should only be used over snow/ice!
     SUBROUTINE Andreas_2002(Z_0,bvisc,ustar,Zt,Zq)
 
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(IN)  :: Z_0, bvisc, ustar
        REAL(kind_phys), INTENT(OUT) :: Zt, Zq
@@ -3313,6 +3629,7 @@ END SUBROUTINE PSI_CB2005
 !! and Holtslag (1991) for stable conditions.
     SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt)
 
+       !$acc routine seq
        IMPLICIT NONE
        REAL(kind_phys), INTENT(OUT)  :: zL
        REAL(kind_phys), INTENT(IN) :: Rib, zaz0, z0zt
@@ -3471,6 +3788,7 @@ REAL(kind_phys) function zolri2(zol2,ri2,za,z0,zt,psi_opt)
 
       REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
 
+      !$acc routine seq
       ! This iterative algorithm to compute z/L from bulk-Ri
 
       IMPLICIT NONE
@@ -3480,7 +3798,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
       REAL(kind_phys) :: zol20,zol3,zolt,zolold
       INTEGER :: n
       INTEGER, PARAMETER :: nmax = 20
-      REAL(kind_phys), DIMENSION(nmax):: zLhux
+      !REAL(kind_phys), DIMENSION(nmax):: zLhux
       REAL(kind_phys) :: psit2,psix2
 
       !print*,"+++++++INCOMING: z/L=",zol1," ri=",ri
@@ -3522,7 +3840,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
         endif
         !print*,"n=",n," psit2=",psit2," psix2=",psix2
         zolrib=ri*psix2**2/psit2
-        zLhux(n)=zolrib
+        !zLhux(n)=zolrib
         n=n+1
       enddo
 
@@ -3530,7 +3848,7 @@ REAL(kind_phys) function zolrib(ri,za,z0,zt,logz0,logzt,zol1,psi_opt)
          !print*,"iter FAIL, n=",n," Ri=",ri," z/L=",zolri
          !if convergence fails, use approximate values:
          CALL Li_etal_2010(zolrib, ri, za/z0, z0/zt)
-         zLhux(n)=zolrib
+         !zLhux(n)=zolrib
          !print*,"FAILED, n=",n," Ri=",ri," z0=",z0
          !print*,"z/L=",zLhux(1:nmax)
       else
@@ -3595,6 +3913,7 @@ END SUBROUTINE psi_init
 !
 !>\ingroup mynn_sfc
    real(kind_phys) function psim_stable_full(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf
 
         !psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
@@ -3605,6 +3924,7 @@ real(kind_phys) function psim_stable_full(zolf)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psih_stable_full(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf
 
         !psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
@@ -3615,6 +3935,7 @@ real(kind_phys) function psih_stable_full(zolf)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psim_unstable_full(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf,x,ym,psimc,psimk
 
         x=(1.-16.*zolf)**.25
@@ -3633,6 +3954,7 @@ real(kind_phys) function psim_unstable_full(zolf)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psih_unstable_full(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf,y,yh,psihc,psihk
 
         y=(1.-16.*zolf)**.5
@@ -3654,6 +3976,7 @@ real(kind_phys) function psih_unstable_full(zolf)
 !>\ingroup mynn_sfc
 !!
    REAL(kind_phys) function psim_stable_full_gfs(zolf)
+        !$acc routine seq
         REAL(kind_phys) :: zolf
         REAL(kind_phys), PARAMETER :: alpha4 = 20.
         REAL(kind_phys) :: aa
@@ -3667,6 +3990,7 @@ REAL(kind_phys) function psim_stable_full_gfs(zolf)
 !>\ingroup mynn_sfc
 !!
    real(kind_phys) function psih_stable_full_gfs(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf
         real(kind_phys), PARAMETER :: alpha4 = 20.
         real(kind_phys) :: bb
@@ -3680,6 +4004,7 @@ real(kind_phys) function psih_stable_full_gfs(zolf)
 !>\ingroup mynn_sfc
 !!
    real(kind_phys) function psim_unstable_full_gfs(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf
         real(kind_phys) :: hl1,tem1
         real(kind_phys), PARAMETER :: a0=-3.975,  a1=12.32,  &
@@ -3700,6 +4025,7 @@ real(kind_phys) function psim_unstable_full_gfs(zolf)
 !>\ingroup mynn_sfc
 !!
    real(kind_phys) function psih_unstable_full_gfs(zolf)
+        !$acc routine seq
         real(kind_phys) :: zolf
         real(kind_phys) :: hl1,tem1
         real(kind_phys), PARAMETER :: a0p=-7.941, a1p=24.75, &
@@ -3720,6 +4046,7 @@ real(kind_phys) function psih_unstable_full_gfs(zolf)
 !>\ingroup mynn_sfc
 !! look-up table functions - or, if beyond -10 < z/L < 10, recalculate
    real(kind_phys) function psim_stable(zolf,psi_opt)
+        !$acc routine seq
         integer :: nzol,psi_opt
         real(kind_phys) :: rzol,zolf
 
@@ -3740,6 +4067,7 @@ real(kind_phys) function psim_stable(zolf,psi_opt)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psih_stable(zolf,psi_opt)
+        !$acc routine seq
         integer :: nzol,psi_opt
         real(kind_phys) :: rzol,zolf
 
@@ -3760,6 +4088,7 @@ real(kind_phys) function psih_stable(zolf,psi_opt)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psim_unstable(zolf,psi_opt)
+        !$acc routine seq
         integer :: nzol,psi_opt
         real(kind_phys) :: rzol,zolf
 
@@ -3780,6 +4109,7 @@ real(kind_phys) function psim_unstable(zolf,psi_opt)
 
 !>\ingroup mynn_sfc
    real(kind_phys) function psih_unstable(zolf,psi_opt)
+        !$acc routine seq
         integer :: nzol,psi_opt
         real(kind_phys) :: rzol,zolf
 
diff --git a/physics/mynnsfc_wrapper.F90 b/physics/mynnsfc_wrapper.F90
index 1a970c9f4..3c033e65e 100644
--- a/physics/mynnsfc_wrapper.F90
+++ b/physics/mynnsfc_wrapper.F90
@@ -191,6 +191,16 @@ SUBROUTINE mynnsfc_wrapper_run(            &
      &            IMS,IME,JMS,JME,KMS,KME,                  &
      &            ITS,ITE,JTS,JTE,KTS,KTE
 
+!$acc enter data create(hfx, znt, psim, psih, chs,          &
+!$acc                   mavail, xland, GZ1OZ0, cpm, qgh,    &
+!$acc                   qfx, snowh_wat)
+
+!$acc enter data create(dz, th, qv)
+
+!$acc enter data copyin(rmol, phii, t3d, exner, qvsh, slmsk, xland)
+
+!$acc enter data copyin(dry, wet, icy, znt_lnd, znt_wat, znt_ice, qsfc_lnd, qsfc_ice, qsfc_lnd_ruc, qsfc_ice_ruc)
+
       ! Initialize CCPP error handling variables
       errmsg = ''
       errflg = 0
@@ -203,6 +213,7 @@ SUBROUTINE mynnsfc_wrapper_run(            &
 !         write(0,*)"iter=",iter
 !      endif
 
+!$acc kernels
       ! prep MYNN-only variables
       dz(:,:) = 0
       th(:,:) = 0
@@ -210,6 +221,9 @@ SUBROUTINE mynnsfc_wrapper_run(            &
       hfx(:)  = 0
       qfx(:)  = 0
       rmol(:) = 0
+!$acc end kernels
+
+!$acc parallel loop collapse(2) present(dz, phii, th, t3d, exner, qv, qvsh)
       do k=1,2 !levs
         do i=1,im
            dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv
@@ -219,6 +233,7 @@ SUBROUTINE mynnsfc_wrapper_run(            &
         enddo
       enddo
 
+!$acc parallel loop present(slmsk, xland, qgh, mavail, cpm, snowh_wat)
       do i=1,im
           if (slmsk(i)==1. .or. slmsk(i)==2.)then !sea/land/ice mask (=0/1/2) in FV3
             xland(i)=1.0                          !but land/water = (1/2) in SFCLAY_mynn
@@ -235,6 +250,7 @@ SUBROUTINE mynnsfc_wrapper_run(            &
           snowh_wat(i) = 0.0
       enddo
 
+!$acc kernels
       ! cm -> m
       where (dry) znt_lnd=znt_lnd*0.01
       where (wet) znt_wat=znt_wat*0.01
@@ -245,6 +261,7 @@ SUBROUTINE mynnsfc_wrapper_run(            &
         where (dry) qsfc_lnd = qsfc_lnd_ruc/(1.+qsfc_lnd_ruc) ! spec. hum
         where (icy) qsfc_ice = qsfc_ice_ruc/(1.+qsfc_ice_ruc) ! spec. hum.
       end if
+!$acc end kernels
 
 !      if (lprnt) then
 !          write(0,*)"CALLING SFCLAY_mynn; input:"
@@ -274,6 +291,8 @@ SUBROUTINE mynnsfc_wrapper_run(            &
 !          write(0,*)"PBLH=",pblh(1)," xland=",xland(1)
 !       endif
 
+!$acc exit data delete(qsfc_lnd_ruc, qsfc_ice_ruc)
+!$acc exit data delete(phii, qvsh, slmsk)
 
         CALL SFCLAY_mynn(                                                     &
              u3d=u,v3d=v,t3d=t3d,qv3d=qv,p3d=prsl,dz8w=dz,                    &
@@ -318,6 +337,13 @@ SUBROUTINE mynnsfc_wrapper_run(            &
              errmsg=errmsg, errflg=errflg                                     )
         if (errflg/=0) return
 
+!$acc exit data delete(hfx, znt, psim, psih, chs,          &
+!$acc                   mavail, xland, GZ1OZ0, cpm, qgh,    &
+!$acc                   qfx, snowh_wat, t3d, exner)
+!$acc exit data delete(dz, th, qv)
+!$acc exit data copyout(rmol)
+!$acc exit data copyout(qsfc_lnd, qsfc_ice)
+
         !! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK:
         !do i = 1, im
         !   !* Taken from sfc_nst.f
@@ -336,10 +362,15 @@ SUBROUTINE mynnsfc_wrapper_run(            &
         !   znt_ice(i)=znt_ice(i)*100.
         !enddo
 
+!$acc kernels
         ! m -> cm
         where (dry) znt_lnd=znt_lnd*100.
         where (wet) znt_wat=znt_wat*100.
         where (icy) znt_ice=znt_ice*100.
+!$acc end kernels
+
+!$acc exit data delete(dry, wet, icy)
+!$acc exit data copyout(znt_lnd, znt_wat, znt_ice)
 
 !      if (lprnt) then
 !         write(0,*)