      subroutine bulk_flux (tile)
      implicit none
      integer*4 tile, trd, omp_get_thread_num
      integer*4  LLm,Lm,MMm,Mm,N, LLm0,MMm0
      parameter (LLm0=510,   MMm0=509,   N=32)
      parameter (LLm=LLm0,  MMm=MMm0)
      integer*4 N_sl
      parameter (N_sl=0)
      integer*4 Lmmpi,Mmmpi,iminmpi,imaxmpi,jminmpi,jmaxmpi
      common /comm_setup_mpi1/ Lmmpi,Mmmpi
      common /comm_setup_mpi2/ iminmpi,imaxmpi,jminmpi,jmaxmpi
      integer*4 NSUB_X, NSUB_E, NPP
      parameter (NPP=1)
      parameter (NSUB_X=1, NSUB_E=NPP)
      integer*4,parameter :: my_acc_device = 0
      logical,parameter :: compute_on_device = .true.
      integer*4 NWEIGHT
      parameter (NWEIGHT=1000)
      integer*4 Ntides
      parameter (Ntides=8)
      integer*4 stdout, Np, NpHz, padd_X,padd_E
      parameter (stdout=6)
      parameter (Np=N+1)
      parameter (NpHz=(N+1+N_sl))
      parameter (Lm=LLm, Mm=MMm)
      parameter (padd_X=(Lm+2)/2-(Lm+1)/2)
      parameter (padd_E=(Mm+2)/2-(Mm+1)/2)
      integer*4 NSA, N2d,N3d,N3dHz, size_XI,size_ETA
      integer*4 se,sse, sz,ssz
      parameter (NSA=28)
      parameter (size_XI=7+(Lm+NSUB_X-1)/NSUB_X)
      parameter (size_ETA=7+(Mm+NSUB_E-1)/NSUB_E)
      parameter (sse=size_ETA/Np, ssz=Np/size_ETA)
      parameter (se=sse/(sse+ssz), sz=1-se)
      parameter (N2d=size_XI*(se*size_ETA+sz*Np))
      parameter (N3d=size_XI*size_ETA*Np)
      parameter (N3dHz=size_XI*size_ETA*NpHz)
      real Vtransform
      parameter (Vtransform=2)
      integer*4   NT, NTA, itemp, NTot
      integer*4   ntrc_temp, ntrc_salt, ntrc_pas, ntrc_bio, ntrc_sed
      integer*4   ntrc_subs, ntrc_substot
      integer*4   ntrc_mld
      parameter (itemp=1)
      parameter (ntrc_temp=1)
      parameter (ntrc_salt=1)
      parameter (ntrc_mld=0)
      parameter (ntrc_pas=0)
      parameter (ntrc_bio=0)
      parameter (ntrc_subs=0, ntrc_substot=0)
      parameter (ntrc_sed=0)
      parameter (NTA=itemp+ntrc_salt)
      parameter (NT=itemp+ntrc_salt+ntrc_pas+ntrc_bio+ntrc_sed+ntrc_mld)
      parameter (NTot=NT)
      integer*4   ntrc_diats, ntrc_diauv, ntrc_diabio
      integer*4   ntrc_diavrt, ntrc_diaek, ntrc_diapv
      integer*4   ntrc_diaeddy, ntrc_surf
     &          , isalt
      parameter (isalt=itemp+1)
      parameter (ntrc_diabio=0)
      parameter (ntrc_diats=0)
      parameter (ntrc_diauv=0)
      parameter (ntrc_diavrt=0)
      parameter (ntrc_diaek=0)
      parameter (ntrc_diapv=0)
      parameter (ntrc_diaeddy=0)
      parameter (ntrc_surf=0)
      real A2d(N2d,NSA,0:NPP-1), A3d(N3d,9,0:NPP-1)
     &    ,A3dHz(N3dHz,4,0:NPP-1)
      integer*4 B2d(N2d,0:NPP-1)
      common/private_scratch/ A2d,A3d,A3dHz
      common/private_scratch_bis/ B2d
      integer*4 chunk_size_X,margin_X,chunk_size_E,margin_E
      integer*4 Istr,Iend,Jstr,Jend, i_X,j_E
      chunk_size_X=(Lm+NSUB_X-1)/NSUB_X
      margin_X=(NSUB_X*chunk_size_X-Lm)/2
      chunk_size_E=(Mm+NSUB_E-1)/NSUB_E
      margin_E=(NSUB_E*chunk_size_E-Mm)/2
      j_E=tile/NSUB_X
      i_X=tile-j_E*NSUB_X
      Istr=1+i_X*chunk_size_X-margin_X
      Iend=Istr+chunk_size_X-1
      Istr=max(Istr,1)
      Iend=min(Iend,Lm)
      Jstr=1+j_E*chunk_size_E-margin_E
      Jend=Jstr+chunk_size_E-1
      Jstr=max(Jstr,1)
      Jend=min(Jend,Mm)
!$acc set device_num(tile)
      trd=omp_get_thread_num()
      call bulk_flux_tile (Istr,Iend,Jstr,Jend,
     &                               A2d(1,1,trd),A2d(1,2,trd))
      return
      end
      subroutine bulk_flux_tile (Istr,Iend,Jstr,Jend, aer,cer)
      implicit none
      integer*4  LLm,Lm,MMm,Mm,N, LLm0,MMm0
      parameter (LLm0=510,   MMm0=509,   N=32)
      parameter (LLm=LLm0,  MMm=MMm0)
      integer*4 N_sl
      parameter (N_sl=0)
      integer*4 Lmmpi,Mmmpi,iminmpi,imaxmpi,jminmpi,jmaxmpi
      common /comm_setup_mpi1/ Lmmpi,Mmmpi
      common /comm_setup_mpi2/ iminmpi,imaxmpi,jminmpi,jmaxmpi
      integer*4 NSUB_X, NSUB_E, NPP
      parameter (NPP=1)
      parameter (NSUB_X=1, NSUB_E=NPP)
      integer*4,parameter :: my_acc_device = 0
      logical,parameter :: compute_on_device = .true.
      integer*4 NWEIGHT
      parameter (NWEIGHT=1000)
      integer*4 Ntides
      parameter (Ntides=8)
      integer*4 stdout, Np, NpHz, padd_X,padd_E
      parameter (stdout=6)
      parameter (Np=N+1)
      parameter (NpHz=(N+1+N_sl))
      parameter (Lm=LLm, Mm=MMm)
      parameter (padd_X=(Lm+2)/2-(Lm+1)/2)
      parameter (padd_E=(Mm+2)/2-(Mm+1)/2)
      integer*4 NSA, N2d,N3d,N3dHz, size_XI,size_ETA
      integer*4 se,sse, sz,ssz
      parameter (NSA=28)
      parameter (size_XI=7+(Lm+NSUB_X-1)/NSUB_X)
      parameter (size_ETA=7+(Mm+NSUB_E-1)/NSUB_E)
      parameter (sse=size_ETA/Np, ssz=Np/size_ETA)
      parameter (se=sse/(sse+ssz), sz=1-se)
      parameter (N2d=size_XI*(se*size_ETA+sz*Np))
      parameter (N3d=size_XI*size_ETA*Np)
      parameter (N3dHz=size_XI*size_ETA*NpHz)
      real Vtransform
      parameter (Vtransform=2)
      integer*4   NT, NTA, itemp, NTot
      integer*4   ntrc_temp, ntrc_salt, ntrc_pas, ntrc_bio, ntrc_sed
      integer*4   ntrc_subs, ntrc_substot
      integer*4   ntrc_mld
      parameter (itemp=1)
      parameter (ntrc_temp=1)
      parameter (ntrc_salt=1)
      parameter (ntrc_mld=0)
      parameter (ntrc_pas=0)
      parameter (ntrc_bio=0)
      parameter (ntrc_subs=0, ntrc_substot=0)
      parameter (ntrc_sed=0)
      parameter (NTA=itemp+ntrc_salt)
      parameter (NT=itemp+ntrc_salt+ntrc_pas+ntrc_bio+ntrc_sed+ntrc_mld)
      parameter (NTot=NT)
      integer*4   ntrc_diats, ntrc_diauv, ntrc_diabio
      integer*4   ntrc_diavrt, ntrc_diaek, ntrc_diapv
      integer*4   ntrc_diaeddy, ntrc_surf
     &          , isalt
      parameter (isalt=itemp+1)
      parameter (ntrc_diabio=0)
      parameter (ntrc_diats=0)
      parameter (ntrc_diauv=0)
      parameter (ntrc_diavrt=0)
      parameter (ntrc_diaek=0)
      parameter (ntrc_diapv=0)
      parameter (ntrc_diaeddy=0)
      parameter (ntrc_surf=0)
      real h(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real hinv(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real f(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real fomn(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /grid_h/h /grid_hinv/hinv /grid_f/f /grid_fomn/fomn
      real angler(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /grid_angler/angler
      real latr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real lonr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real latu(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real lonu(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real latv(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real lonv(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /grid_latr/latr /grid_lonr/lonr
      common /grid_latu/latu /grid_lonu/lonu
      common /grid_latv/latv /grid_lonv/lonv
      real pm(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pn(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real om_r(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real on_r(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real om_u(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real on_u(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real om_v(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real on_v(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real om_p(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real on_p(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pn_u(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pm_v(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pm_u(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pn_v(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /metrics_pm/pm    /metrics_pn/pn
      common /metrics_omr/om_r /metrics_on_r/on_r
      common /metrics_omu/om_u /metrics_on_u/on_u
      common /metrics_omv/om_v /metrics_on_v/on_v
      common /metrics_omp/om_p /metrics_on_p/on_p
      common /metrics_pnu/pn_u /metrics_pmv/pm_v
      common /metrics_pmu/pm_u /metrics_pnv/pn_v
      real dmde(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real dndx(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /metrics_dmde/dmde    /metrics_dndx/dndx
      real pmon_p(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pmon_r(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pmon_u(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pnom_p(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pnom_r(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pnom_v(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real grdscl(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /metrics_pmon_p/pmon_p /metrics_pnom_p/pnom_p
      common /metrics_pmon_r/pmon_r /metrics_pnom_r/pnom_r
      common /metrics_pmon_u/pmon_u /metrics_pnom_v/pnom_v
      common /metrics_grdscl/grdscl
      real rmask(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pmask(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real umask(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real vmask(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real pmask2(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /mask_r/rmask
      common /mask_p/pmask
      common /mask_u/umask
      common /mask_v/vmask
      common /mask_p2/pmask2
      real zob(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /Z0B_VAR/zob
      real u(0:Lm+1+padd_X,0:Mm+1+padd_E,N,3)
      real v(0:Lm+1+padd_X,0:Mm+1+padd_E,N,3)
      real t(0:Lm+1+padd_X,0:Mm+1+padd_E,N,3,NT)
      common /ocean_u/u /ocean_v/v /ocean_t/t
      real Hz(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      real Hz_bak(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      real z_r(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      real z_w(0:Lm+1+padd_X,0:Mm+1+padd_E,0:N)
      real Huon(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      real Hvom(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      common /grid_Hz_bak/Hz_bak /grid_zw/z_w /grid_Huon/Huon
      common /grid_Hvom/Hvom
      real We(0:Lm+1+padd_X,0:Mm+1+padd_E,0:N)
      common /grid_Hz/Hz /grid_zr/z_r /grid_We/We
      real rho1(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      real rho(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      common /ocean_rho1/rho1 /ocean_rho/rho
      real qp1(0:Lm+1+padd_X,0:Mm+1+padd_E,N)
      common /ocean_qp1/qp1
      real qp2
      parameter (qp2=0.0000172D0)
      real sustr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real svstr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /forces_sustr/sustr /forces_svstr/svstr
      real sustrg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real svstrg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      common /smsdat_sustrg/sustrg /smsdat_svstrg/svstrg
      real    sustrp(2), svstrp(2), sms_time(2)
      real    sms_cycle, sms_scale
      integer*4 itsms, sms_ncycle, sms_rec, lsusgrd
      integer*4 lsvsgrd,sms_tid, susid, svsid
      real    sms_origin_date_in_sec
      common /smsdat1/ sustrp, svstrp, sms_time
      common /smsdat2/ sms_origin_date_in_sec
      common /smsdat3/ sms_cycle, sms_scale
      common /smsdat4/ itsms, sms_ncycle, sms_rec, lsusgrd
      common /smsdat5/ lsvsgrd,sms_tid, susid, svsid
      real wspd_cfb(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /smsdat_wspd_cfb/ wspd_cfb
      integer*4 lwgrd, wid
      common /smsdat5/ lwgrd, wid
      real bustr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real bvstr(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /forces_bustr/bustr /forces_bvstr/bvstr
      real bustrg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real bvstrg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      common /bmsdat_bustrg/bustrg /bmsdat_bvstrg/bvstrg
      real bms_tintrp(2), bustrp(2),    bvstrp(2), tbms(2)
      real bmsclen, bms_tstart, bms_tend,  tsbms, sclbms
      integer*4 itbms,      bmstid,busid, bvsid,     tbmsindx
      logical bmscycle,   bms_onerec,   lbusgrd,   lbvsgrd
      common /bmsdat1/bms_tintrp, bustrp,       bvstrp,    tbms
      common /bmsdat2/bmsclen, bms_tstart, bms_tend, tsbms, sclbms
      common /bmsdat3/itbms,      bmstid,busid, bvsid,     tbmsindx
      common /bmsdat4/bmscycle,   bms_onerec,   lbusgrd,   lbvsgrd
      real stflx(0:Lm+1+padd_X,0:Mm+1+padd_E,NT)
      common /forces_stflx/stflx
      real shflx_rsw(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /frc_shflx_rsw/shflx_rsw
      real shflx_rlw(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /frc_shflx_rlw/shflx_rlw
      real shflx_lat(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /frc_shflx_lat/shflx_lat
      real shflx_sen(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /frc_shflx_sen/shflx_sen
      real stflxg(0:Lm+1+padd_X,0:Mm+1+padd_E,2,NT)
      common /stfdat_stflxg/stflxg
      real stflxp(2,NT), stf_time(2,NT)
      real stf_cycle(NT), stf_scale(NT)
      integer*4 itstf(NT), stf_ncycle(NT), stf_rec(NT)
      integer*4 lstfgrd(NT), stf_tid(NT), stf_id(NT)
      REAL(kind=8) :: stf_origin_date_in_sec
      common /stfdat1/ stflxp,  stf_time, stf_cycle, stf_scale
      common /stfdat2/ stf_origin_date_in_sec
      common /stfdat3/ itstf, stf_ncycle, stf_rec, lstfgrd
      common /stfdat4/  stf_tid, stf_id
      real btflx(0:Lm+1+padd_X,0:Mm+1+padd_E,NT)
      common /forces_btflx/btflx
      real tair(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real rhum(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real prate(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real radlw(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real radsw(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real wspd(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real uwnd(0:Lm+1+padd_X,0:Mm+1+padd_E)
      real vwnd(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /bulk_tair/ tair
      common /bulk_rhum/ rhum
      common /bulk_prate/ prate
      common /bulk_radlw/ radlw
      common /bulk_radsw/ radsw
      common /bulk_wspd/ wspd
      common /bulk_uwnd/ uwnd
      common /bulk_vwnd/ vwnd
      real tairg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real rhumg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real prateg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real radlwg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real radswg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real uwndg_norot(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real radswg_down(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real uwndg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real vwndg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      real wspdg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      common /bulkdat_tairg/tairg
      common /bulkdat_rhumg/rhumg
      common /bulkdat_prateg/prateg
      common /bulkdat_radlwg/radlwg
      common /bulkdat_radswg/radswg
      common /bulk_uwndg_norot/uwndg_norot
      common /bulk_radswg_down/radswg_down
      common /bulk_uwndg/uwndg
      common /bulk_vwndg/vwndg
      common /bulkdat_wspdg/wspdg
      real    tairp(2),rhump(2),pratep(2),radlwp(2),radswp(2)
      real    uwndp(2),vwndp(2)
      real    bulk_time(2), bulk_cycle
      integer*4 tair_id,rhum_id,prate_id,radlw_id,radsw_id
      integer*4 ltairgrd,lrhumgrd,lprategrd,lradlwgrd,lradswgrd
      REAL(kind=8) :: blk_origin_date_in_sec
      integer*4 uwnd_id,vwnd_id,luwndgrd,lvwndgrd
      integer*4 itbulk,bulk_ncycle,bulk_rec,bulk_tid
      integer*4 bulkunused
      common /bulkdat1_for/ tair_id,rhum_id,prate_id,radlw_id,radsw_id
      common /bulkdat1_grd
     &                 / ltairgrd,lrhumgrd,lprategrd,lradlwgrd,lradswgrd
      common /bulkdat1_tim/ itbulk, bulk_ncycle, bulk_rec, bulk_tid
      common /bulkdat1_uns/ bulkunused
      common /bulkdat1_wnd/ uwnd_id,vwnd_id,luwndgrd,lvwndgrd
      common /bulkdat2_for/ tairp,rhump,pratep,radlwp,radswp
      common /bulkdat2_tim/ bulk_time, bulk_cycle, 
     &                                            blk_origin_date_in_sec
      common /bulkdat2_wnd/ uwndp,vwndp
      real srflx(0:Lm+1+padd_X,0:Mm+1+padd_E)
      common /forces_srflx/srflx
      real srflxg(0:Lm+1+padd_X,0:Mm+1+padd_E,2)
      common /srfdat_srflxg/srflxg
      real srflxp(2),srf_time(2)
      real srf_cycle, srf_scale
      integer*4 itsrf, srf_ncycle, srf_rec
      integer*4 lsrfgrd, srf_tid, srf_id
      REAL(kind=8) :: srf_origin_date_in_sec
      common /srfdat1/ srflxp, srf_time, srf_cycle, srf_scale
      common /srfdat2/ srf_origin_date_in_sec
      common /srfdat3/ itsrf,srf_ncycle,srf_rec,lsrfgrd,srf_tid,srf_id
      real dt, dtfast, time, time2, time_start, tdays, start_time
      integer*4 ndtfast, iic, kstp, krhs, knew, next_kstp
     &      , iif, nstp, nrhs, nnew, nbstep3d
      logical PREDICTOR_2D_STEP
      common /time_indices/  dt,dtfast, time, time2,time_start, tdays,
     &     ndtfast, iic, kstp, krhs, knew, next_kstp,
     &     start_time,
     &                       iif, nstp, nrhs, nnew, nbstep3d,
     &                       PREDICTOR_2D_STEP
      real time_avg, time2_avg, rho0
     &               , rdrg, rdrg2, Cdb_min, Cdb_max, Zobt
     &               , xl, el, visc2, visc4, gamma2
      real  theta_s,   theta_b,   Tcline,  hc
      real  sc_w(0:N), Cs_w(0:N), sc_r(N), Cs_r(N)
      real  rx0, rx1
      real  tnu2(NT),tnu4(NT)
      real weight(6,0:NWEIGHT)
      real  x_sponge,   v_sponge
       real  tauT_in, tauT_out, tauM_in, tauM_out
      integer*4 numthreads,     ntstart,   ntimes,  ninfo
     &      , nfast,  nrrec,     nrst,    nwrt
     &                                 , ntsavg,  navg
      logical ldefhis
      logical got_tini(NT)
      common /scalars_main/
     &             time_avg, time2_avg,  rho0,      rdrg,    rdrg2
     &           , Zobt,       Cdb_min,   Cdb_max
     &           , xl, el,    visc2,     visc4,   gamma2
     &           , theta_s,   theta_b,   Tcline,  hc
     &           , sc_w,      Cs_w,      sc_r,    Cs_r
     &           , rx0,       rx1
     &           ,       tnu2,    tnu4
     &                      , weight
     &                      , x_sponge,   v_sponge
     &                      , tauT_in, tauT_out, tauM_in, tauM_out
     &      , numthreads,     ntstart,   ntimes,  ninfo
     &      , nfast,  nrrec,     nrst,    nwrt
     &                                 , ntsavg,  navg
     &                      , got_tini
     &                      , ldefhis
      logical synchro_flag
      common /sync_flag/ synchro_flag
      integer*4 may_day_flag
      integer*4 tile_count, first_time, bc_count
      common /communicators_i/
     &        may_day_flag, tile_count, first_time, bc_count
      real hmin, hmax, grdmin, grdmax, Cu_min, Cu_max
      common /communicators_r/
     &     hmin, hmax, grdmin, grdmax, Cu_min, Cu_max
      real lonmin, lonmax, latmin, latmax
      common /communicators_lonlat/
     &     lonmin, lonmax, latmin, latmax
      real*8 Cu_Adv3d,  Cu_W, Cu_Nbq_X, Cu_Nbq_Y, Cu_Nbq_Z
      integer*4 i_cx_max, j_cx_max, k_cx_max
      common /diag_vars/ Cu_Adv3d,  Cu_W,
     &        i_cx_max, j_cx_max, k_cx_max
      real*8 volume, avgke, avgpe, avgkp, bc_crss
      common /communicators_rq/
     &          volume, avgke, avgpe, avgkp, bc_crss
      real*4 CPU_time(0:31,0:NPP)
      integer*4 proc(0:31,0:NPP),trd_count
      common /timers_roms/CPU_time,proc,trd_count
      real pi, deg2rad, rad2deg
      parameter (pi=3.14159265358979323846D0, deg2rad=pi/180.D0,
     &                                      rad2deg=180.D0/pi)
      real Eradius, Erotation, g, day2sec,sec2day, jul_off,
     &     year2day,day2year
      parameter (Eradius=6371315.0D0,  Erotation=7.292115090D-5,
     &           day2sec=86400.D0, sec2day=1.D0/86400.D0,
     &           year2day=365.25D0, day2year=1.D0/365.25D0,
     &           jul_off=2440000.D0)
      parameter (g=9.81D0)
      real Cp
      parameter (Cp=3985.0D0)
      real vonKar
      parameter (vonKar=0.41D0)
      real spval
      parameter (spval=-999.0D0)
      logical mask_val
      parameter (mask_val = .true.)
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      INTEGER*4              :: i, j, m, iter
      INTEGER*4              :: Istr, Iend, Jstr, Jend
      INTEGER*4              :: imin, imax, jmin, jmax
      INTEGER*4, PARAMETER   :: IterFl = 3
      INTEGER*4, PARAMETER :: mb = 2
      REAL           :: cff  , cff2
      REAL           :: rho0i, cpi  , RH   , Hlv     , wspd0
      REAL           :: TseaC, TseaK, Qsea , spec_hum, Tvstar
      REAL           :: TairC, TairK, Q    , rhoAir  , Bf
      REAL           :: delW(2) , delT , delQ , ddelW(2) , ddelT , ddelQ
      REAL           :: Cd   , Ch   , Ce   , qsat    , iexns , iexna
      REAL           :: Wstar(2), Tstar, Qstar, Wstar0  , Tstar0, Qstar0
      REAL           :: psi_u(2), ZoLu(2) , psi_t, ZoLt    , patm
      REAL           :: hfsen, hflat, hflw , upvel   , evap
!$acc declare create( cff  , cff2
!$acc&               ,rho0i, cpi  , RH   , Hlv     , wspd0
!$acc&               ,TseaC, TseaK, Qsea , Tvstar
!$acc&               ,TairC, TairK, Q    , rhoAir  , Bf
!$acc&               ,delW , delT , delQ , ddelW , ddelT , ddelQ
!$acc&               ,Cd   , Ch   , Ce   , iexns , iexna
!$acc&               ,Wstar, Tstar, Qstar, Wstar0  , Tstar0, Qstar0
!$acc&               ,psi_u, ZoLu, psi_t, ZoLt    , patm
!$acc&               ,hfsen, hflat, hflw , upvel   , evap
!$acc&               )
      REAL           :: wspd0_cfb
      REAL            :: blk_ZW, blk_ZT, blk_ZToZW
      PARAMETER ( blk_ZW    =  10.0D0 )
      PARAMETER ( blk_ZT    =   2.0D0 )
      PARAMETER ( blk_ZToZW = blk_ZT/blk_ZW )
      REAL            ::  logus10,logts10
      REAL           :: Log10oLogZw
      PARAMETER ( Log10oLogZw = LOG(  10.0D0*10000.D0)/
     &                          LOG(blk_ZW*10000.D0) )
      REAL           :: bulk_psiu_coare, bulk_psit_coare, air_visc
      REAL           :: iZoW , iZoT  , CC, Ch10, Ribcu , charn
      REAL           :: iZo10, iZoT10, Ri, Rr  , VisAir
      REAL           ::  aer(Istr-2:Iend+2,Jstr-2:Jend+2)
      REAL           ::  cer(Istr-2:Iend+2,Jstr-2:Jend+2)
      REAL           ::  stau(Istr-2:Iend+2,Jstr-2:Jend+2)
!$acc declare create( aer, cer
!$acc&                ,stau
!$acc&               )
      if (istr.eq.1) then
        imin=Istr-1
      else
        imin=Istr-2
      endif
      if (iend.eq.Lm) then
        imax=Iend+1
      else
        imax=Iend+2
      endif
      if (jstr.eq.1) then
        jmin=Jstr-1
      else
        jmin=Jstr-2
      endif
      if (jend.eq.Mm) then
        jmax=Jend+1
      else
        jmax=Jend+2
      endif
      rho0i=1.0D0/rho0
      cpi=1.0D0/cp
      iexns  = (psurf*ip00)**(-rdocpd)
!$acc update device( iexns, iexna, rho0i, hfsen, cpi )
!$acc kernels if(compute_on_device) default(present)
!$acc loop seq
      DO j=jmin,jmax
!$acc loop seq
        DO i=imin,imax
          wspd0     = wspd(i,j)
          wspd0_cfb = wspd_cfb(i,j)
          wspd0_cfb = MAX ( wspd0_cfb , 0.1D0 * MIN(10.D0, blk_ZW) )
          TairC    = tair(i,j)
          TairK    = TairC + CtoK
          RH       = rhum(i,j)
          Q        = spec_hum (RH,psurf,TairC)
          wspd0     = MAX ( wspd0 , 0.1D0 * MIN(10.D0, blk_ZW) )
          TseaC = t(i,j,N,nrhs,itemp)
          TseaK = TseaC + CtoK
          CALL exner_patm_from_tairabs(iexna,patm,Q,TairK,blk_ZT,psurf)
          rhoAir = patm*(1.D0+Q) / ( blk_Rgas*TairK*(1.D0+MvoMa*Q) )
          Qsea   = qsat(TseaK,psurf,0.98D0)
          delW(1) = SQRT(wspd0*wspd0+0.25D0)
          delW(2) = SQRT(wspd0_cfb*wspd0_cfb+0.25D0)
          delQ = Q-Qsea
          cff  = CtoK*(iexna-iexns)
          delT = TairC*iexna - TseaC*iexns + cff
          Wstar(1:mb)  = 0.035D0*delW(1:mb)*Log10oLogZw
          VisAir = air_visc ( TairC )
          charn  = 0.011D0
          Ch10   = 0.00115D0
          Ribcu  = -    blk_ZW / ( blk_Zabl * 0.004D0 * blk_beta**3 )
!$acc loop seq
          DO m=1,mb
            iZo10  = g*Wstar(m) /
     &               (charn*Wstar(m)*Wstar(m)*Wstar(m)+0.11D0*g*VisAir)
            iZoT10 = 0.1D0 * exp(vonKar*vonKar /
     &                     ( Ch10*LOG( 10.0D0*iZo10 ) ) )
            CC     = LOG( blk_ZW*iZo10 )*LOG( blk_ZW*iZo10 )/
     &               LOG( blk_ZT*iZoT10 )
            Ri     =  g * blk_ZW * ( delT+cpvir*TairK*delQ )/
     &                             ( TairK*delW(m)*delW(m) )
            IF ( Ri < 0.0D0 ) THEN
              ZoLu(m)=CC*Ri/(1.0D0+Ri/Ribcu)
            ELSE
              ZoLu(m)=CC*Ri/(1.0D0+3.0D0*Ri/CC)
            ENDIF
            psi_u(m)     = bulk_psiu_coare ( ZoLu(m) )
            logus10   = LOG(blk_ZW* iZo10)
            Wstar(m)  = delW(m)*vonKar/(logus10-psi_u(m))
          ENDDO
          ZoLt      = ZoLu(mb)*blk_ZToZW
          psi_t     = bulk_psit_coare ( ZoLt    )
          logts10   = LOG(blk_ZT*iZoT10)
          cff       =      vonKar/(logts10-psi_t)
          Tstar     = delT*cff
          Qstar     = delQ*cff
          charn = 0.011D0
          IF     ( delW(1) > 18.0D0 ) then
            charn = 0.018D0
          ELSEIF ( delW(1) > 10.0D0 ) then
            charn = 0.011D0+0.125D0*(0.018D0-0.011D0)*(delW(1)-10.D0)
          ENDIF
!$acc loop seq
          DO Iter=1,IterFl
!$acc loop seq
            DO m=1,mb
              iZoW    = g*Wstar(m) / ( charn*Wstar(m)*Wstar(m)*Wstar(m)
     &                                 +0.11D0*g*VisAir )
              Rr      = Wstar(m)/(iZow*VisAir)
              iZoT    = MAX(8695.65D0,18181.8D0*(Rr**0.6D0))
              ZoLu(m)    = vonKar*g*blk_ZW*
     &             (Tstar*(1.0D0+cpvir*Q)+cpvir*TairK*Qstar)/
     &             (TairK*Wstar(m)*Wstar(m)*(1.0D0+cpvir*Q)+eps)
              psi_u(m)     = bulk_psiu_coare ( ZoLu(m) )
              logus10 = LOG(blk_ZW*iZoW)
              Wstar(m)   = delW(m)*vonKar/(logus10-psi_u(m))
            ENDDO
            ZoLt      = ZoLu(mb)*blk_ZToZW
            psi_t     = bulk_psit_coare ( ZoLt )
            logts10   = LOG(blk_ZT*iZoT)
            cff       = vonKar/(logts10-psi_t)
            Tstar     = delT*cff
            Qstar     = delQ*cff
              Bf=-g/TairK*Wstar(1)*(Tstar+cpvir*TairK*Qstar)
              if (Bf.gt.0.0D0) then
                cff=blk_beta*(Bf*blk_Zabl)**r3
              else
                cff=0.2D0
              endif
              delW(1)  = SQRT(wspd0*wspd0+cff*cff)
              Bf=-g/TairK*Wstar(2)*(Tstar+cpvir*TairK*Qstar)
              if (Bf.gt.0.0D0) then
                cff=blk_beta*(Bf*blk_Zabl)**r3
              else
                cff=0.2D0
              endif
              delW(2)  = SQRT(wspd0_cfb*wspd0_cfb+cff*cff)
            ddelW(1:mb) = delW(1:mb)
          ENDDO
          hfsen = - blk_Cpa*rhoAir*Wstar(mb)*Tstar
          Hlv   = (2.5008D0 - 0.0023719D0*TseaC)*1.0D+6
          hflat = - Hlv*rhoAir*Wstar(mb)*Qstar
          hflw   = radlw(i,j)
     &           - emiss_lw*rho0i*cpi*SigmaSB*TseaK*TseaK*TseaK*TseaK
          upvel=-1.61D0*Wstar(mb)*Qstar-(1.0D0+1.61D0*Q)*Wstar(mb)
     &                                            *Tstar/TairK
          hflat=hflat+rhoAir*Hlv*upvel*Q
          hflat=-hflat*rho0i*cpi
          hfsen=-hfsen*rho0i*cpi
          stflx(i,j,itemp)=srflx(i,j)+hflw+hflat+hfsen
          evap=-cp*hflat/Hlv
          stflx(i,j,isalt)=(evap-prate(i,j))*t(i,j,N,nrhs,isalt)
          stflx(i,j,itemp)=stflx(i,j,itemp)*rmask(i,j)
          stflx(i,j,isalt)=stflx(i,j,isalt)*rmask(i,j)
          aer(i,j)  = rhoAir*delW(1)
          Cd        = (Wstar(1)/ddelW(1))**2
          cer(i,j)  = Cd
          stau(i,j) = cfb_slope * delW(1) + cfb_offset
          shflx_rsw(i,j)=srflx(i,j)
          shflx_lat(i,j)=hflat
          shflx_sen(i,j)=hfsen
          shflx_rlw(i,j)=hflw
        ENDDO
      ENDDO
!$acc loop seq
      do j=jmin,jmax
!$acc loop seq
        do i=imin+1,imax
          cff =0.5D0*(cer(i-1,j)+cer(i,j))
          cff2=0.5D0*(aer(i-1,j)+aer(i,j))
          sustr(i,j)= ( cff*cff2*uwnd(i,j)
     &     + 0.5D0*(stau(i,j)+stau(i-1,j))*u(i,j,N,nrhs)
     &                             )*rho0i
          sustr(i,j)=sustr(i,j)*umask(i,j)
        enddo
      enddo
!$acc end kernels
!$acc kernels if(compute_on_device) default(present)
!$acc loop seq
      do j=jmin+1,jmax
!$acc loop seq
        do i=imin,imax
          cff =0.5D0*(cer(i,j-1)+cer(i,j))
          cff2=0.5D0*(aer(i,j-1)+aer(i,j))
          svstr(i,j)=( cff*cff2*vwnd(i,j)
     &   + 0.5D0*(stau(i,j)+stau(i,j-1))*v(i,j,N,nrhs)
     &                                        )*rho0i
          svstr(i,j)=svstr(i,j)*vmask(i,j)
        enddo
      enddo
!$acc end kernels
      return
      end
      FUNCTION spec_hum (RH,psfc,TairC)
!$acc routine(spec_hum) seq
      IMPLICIT NONE
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      REAL       ::  RH   , spec_hum , cff
      REAL       ::  psfc, TairC
!$acc declare create(  RH,psfc,TairC )
      cff=(1.0007D0+3.46D-6*0.01D0*psfc)*6.1121D0*
     &        exp(17.502D0*TairC/(240.97D0+TairC))
      IF (RH.lt.2.0D0) then
         cff=cff*RH
         spec_hum=MvoMa*(cff/(psfc*0.01D0-0.378D0*cff))
      ELSE
         spec_hum=0.001D0*RH
      ENDIF
      END FUNCTION spec_hum
      SUBROUTINE exner_patm_from_tairabs (iexn,pair,q,tairabs,z,psfc)
!$acc routine(exner_patm_from_tairabs) seq
      IMPLICIT NONE
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      REAL,INTENT(  out)   :: iexn,pair
      REAL,INTENT(in   )   :: q, tairabs, z, psfc
      REAL                 :: xm,q_sat,qsat
      REAL, PARAMETER      ::  g    = 9.80665D0
      INTEGER*4              :: iter
      INTEGER*4, PARAMETER   :: Niter = 3
      pair = psfc
!$acc loop seq
      DO Iter = 1, Niter
        q_sat = qsat(tairabs, pair, 1.D0)
        xm    =  mm_dryair + (q/q_sat) * ( mm_water - mm_dryair )
        pair  = psfc * EXP( -g * xm * z / ( r_gas * tairabs ) )
      ENDDO
      iexn =  (pair*ip00)**(-rdocpd)
      return
      END SUBROUTINE exner_patm_from_tairabs
      FUNCTION qsat (TairK, patm, coeff)
!$acc routine(qsat) seq
      IMPLICIT NONE
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      REAL                 ::  qsat
      REAL                 ::  TairK, patm, coeff
      REAL                 ::  psat
      REAL, PARAMETER      ::  alpw    = 60.2227554D0
      REAL, PARAMETER      ::  betaw   = 6822.40088D0
      REAL, PARAMETER      ::  gamw    = 5.13926744D0
      REAL, PARAMETER      ::  alpi    = 32.62117980819471D0
      REAL, PARAMETER      ::  betai   = 6295.421338904806D0
      REAL, PARAMETER      ::  gami    = 0.5631331575423155D0
      IF (TairK .LE. CtoK) then
        psat = EXP( alpi - betai/TairK - gami*LOG(TairK) )
      ELSE
        psat = EXP( alpw - betaw/TairK - gamw*LOG(TairK) )
      ENDIF
      psat = coeff * psat
      qsat = (MvoMa*psat)/(patm+(MvoMa-1.0D0)*psat)
      return
      END FUNCTION qsat
      FUNCTION air_visc(TairC)
!$acc routine(air_visc) seq
      REAL                 :: air_visc,cff
      REAL, PARAMETER      :: c0 = 1.326D-5
      REAL, PARAMETER      :: c1 = 6.542D-3
      REAL, PARAMETER      :: c2 = 8.301D-6
      REAL, PARAMETER      :: c3 = 4.84D-9
      cff      = TairC*TairC
      air_visc = c0*(1.D0+c1*TairC+c2*cff-c3*cff*TairC)
      return
      END FUNCTION air_visc
      FUNCTION bulk_psiu_coare (ZoL)
!$acc routine(bulk_psiu_coare) seq
      IMPLICIT NONE
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      REAL                 ::  bulk_psiu_coare
      REAL                 ::  ZoL
      REAL                 ::  chik, psik
      REAL                 ::  chic, psic
      IF (ZoL <= 0.0D0) then
        chik = (1.0D0-15.0D0*ZoL)**0.25D0
        psik = 2.0D0*LOG(0.5D0*(1.0D0+chik))+LOG(0.5D0*(1.0D0+chik**2))
     &                                  -2.0D0*ATAN(chik)+pis2
        chic = (1.0D0-10.15D0*ZoL)**r3
        psic  = 1.5D0*LOG(r3*(chic**2+chic+1.0D0))
     &          - sqr3*ATAN((2.0D0*chic+1.0D0)/sqr3)+2.D0*pis2osqr3
        bulk_psiu_coare=psic+(psik-psic)/(1.0D0+ZoL**2)
      ELSE
        chic=-MIN(50.0D0,0.35D0*ZoL)
        bulk_psiu_coare=-((1.0D0+ZoL)+0.6667D0*(ZoL-14.28D0)*EXP(chic)+
     &                                                          8.525D0)
      ENDIF
      return
      END FUNCTION bulk_psiu_coare
      FUNCTION bulk_psit_coare (ZoL)
!$acc routine(bulk_psit_coare) seq
      IMPLICIT NONE
      REAL, PARAMETER ::  CtoK      =  273.16D0
      REAL, PARAMETER ::  blk_Rgas  =  287.0596736665907D0
      REAL, PARAMETER ::  blk_Rvap  =  461.5249933083879D0
      REAL, PARAMETER ::  blk_Cpa   = 1004.708857833067D0
      REAL, PARAMETER ::  ip00      =    1.D-5
      REAL, PARAMETER ::  p00       =    1.D+5
      REAL, PARAMETER ::  rdocpd    = blk_Rgas/blk_Cpa
      REAL, PARAMETER ::  cpdord    = blk_Cpa/blk_Rgas
      REAL, PARAMETER ::  r_gas     = 8.314510D0
      REAL, PARAMETER ::  mm_dryair = 28.9644D-3
      REAL, PARAMETER ::  mm_water  = 18.0153D-3
      REAL, PARAMETER ::  cpvir     = blk_Rvap/blk_Rgas - 1.D0
      REAL, PARAMETER ::  MvoMa     = mm_water/mm_dryair
      REAL, PARAMETER ::  eps       = 1.D-8
      REAL, PARAMETER ::  r3        = 1.D0/3.D0
      REAL, PARAMETER ::  pis2      = 2.D0*ATAN(1.D0)
      REAL, PARAMETER ::  sqr3      = SQRT(3.D0)
      REAL, PARAMETER ::  pis2osqr3 = pis2/sqr3
      REAL, PARAMETER ::  LMOmin    = -200.D0
      REAL, PARAMETER ::  LMOmax    = 0.25D0
      REAL, PARAMETER ::  blk_beta  = 1.2D0
      REAL, PARAMETER ::  blk_Zabl  = 600.D0
      REAL, PARAMETER ::  dWstar0   = 1.D-6
      REAL, PARAMETER ::  dTstar0   = 1.D-6
      REAL, PARAMETER ::  dQstar0   = 1.D-9
      REAL, PARAMETER ::  emiss_lw  =  0.985D0
      REAL, PARAMETER ::  SigmaSB   =  5.6697D-8
      REAL            :: psurf
      PARAMETER (psurf=100000.0D0)
      real swparam
      parameter (swparam=0.3D0)
      real cfb_slope, cfb_offset
      parameter (cfb_slope=-0.0029D0)
      parameter (cfb_offset=0.008D0)
      REAL                 ::  bulk_psit_coare
      REAL                 ::  ZoL
      REAL                 ::  chik, psik
      REAL                 ::  chic, psic
      IF (ZoL < 0.0D0) THEN
        chik  = (1.0D0-15.0D0*ZoL)**0.25D0
        psik  = 2.0D0*LOG(0.5D0*(1.0D0+chik**2))
        chic  = (1.0D0-34.15D0*ZoL)**r3
        psic  = 1.5D0*LOG((chic**2+chic+1.0D0)*r3)
     &               -sqr3*ATAN((2.0D0*chic+1.0D0)/sqr3)
     &               +2.D0*pis2osqr3
        bulk_psit_coare = psic+(psik-psic)/(1.0D0+ZoL**2)
      ELSE
        chic=-MIN(50.0D0,0.35D0*ZoL)
        bulk_psit_coare = -((1.0D0+2.0D0*ZoL/3.0D0)**1.5D0+
     &            0.6667D0*(ZoL-14.28D0)*EXP(chic)+8.525D0)
      ENDIF
      return
      END FUNCTION bulk_psit_coare
