Added wave forcing fields.

cesm/mod_cesm.F90

@@ -28,8 +28,8 @@ module mod_cesm
    use mod_time, only: nstep
    use mod_xc
    use mod_forcing, only: trxday, srxday, swa, nsf, lip, sop, eva, rnf, rfi, &
-                          fmltfz, sfl, ztx, mty, ustarw, slp, abswnd, atmco2,&
-                          atmbrf
+                          fmltfz, sfl, ztx, mty, ustarw, slp, abswnd, &
+                          lamult, lasl, ustokes, vstokes, atmco2, atmbrf
    use mod_ben02, only: initai, rdcsic, rdctsf, fnlzai
    use mod_seaice, only: ficem
    use mod_checksum, only: csdiag, chksummsk
@@ -70,9 +70,13 @@ module mod_cesm
       ustarw_da, &       ! Friction velocity for open water [m s-1].
       slp_da, &          ! Sea-level pressure [kg m-1 s-2].
       abswnd_da, &       ! Wind speed at measurement height (zu) [m s-1].
+      ficem_da, &        ! Ice concentration [].
+      lamult_da, &       ! Langmuir enhancement factor [].
+      lasl_da, &         ! Surface layer averaged Langmuir number [].
+      ustokes_da, &      ! u-component of surface Stokes drift [m s-1].
+      vstokes_da, &      ! v-component of surface Stokes drift [m s-1].
       atmco2_da, &       ! Atmospheric CO2 concentration [ppm].
-      atmbrf_da, &       ! Atmospheric bromoform concentration [ppt].
-      ficem_da           ! Ice concentration [].
+      atmbrf_da          ! Atmospheric bromoform concentration [ppt].
 
    logical :: &
       smtfrc             ! If true, time smooth CESM forcing fields.
@@ -83,7 +87,8 @@ module mod_cesm
    public :: runid_cesm, runtyp_cesm, ocn_cpl_dt_cesm, nstep_in_cpl, hmlt, &
              frzpot, mltpot, swa_da, nsf_da, hmlt_da, lip_da, sop_da, eva_da, &
              rnf_da, rfi_da, fmltfz_da, sfl_da, ztx_da, mty_da, ustarw_da, &
-             slp_da, abswnd_da, atmco2_da, atmbrf_da, ficem_da, smtfrc, l1ci, l2ci, &
+             slp_da, abswnd_da, ficem_da, lamult_da, lasl_da, &
+             ustokes_da, vstokes_da, atmco2_da, atmbrf_da, smtfrc, l1ci, l2ci, &
              inicon_cesm, inifrc_cesm, getfrc_cesm
 
 contains
@@ -146,6 +151,7 @@ subroutine getfrc_cesm
 #undef DIAG
 #ifdef DIAG
       use mod_nctools
+      use mod_dia, only : iotype
 #endif
 
       integer :: i, j, l
@@ -163,22 +169,26 @@ subroutine getfrc_cesm
       do j = 1, jj
         do l = 1, isp(j)
         do i = max(1, ifp(j, l)), min(ii, ilp(j, l))
-           ustarw(i, j) = w1*ustarw_da(i, j, l1ci) + w2*ustarw_da(i, j, l2ci)
-           lip(i, j)    = w1*lip_da(i, j, l1ci)    + w2*lip_da(i, j, l2ci)
-           sop(i, j)    = w1*sop_da(i, j, l1ci)    + w2*sop_da(i, j, l2ci)
-           eva(i, j)    = w1*eva_da(i, j, l1ci)    + w2*eva_da(i, j, l2ci)
-           rnf(i, j)    = w1*rnf_da(i, j, l1ci)    + w2*rnf_da(i, j, l2ci)
-           rfi(i, j)    = w1*rfi_da(i, j, l1ci)    + w2*rfi_da(i, j, l2ci)
-           fmltfz(i, j) = w1*fmltfz_da(i, j, l1ci) + w2*fmltfz_da(i, j, l2ci)
-           sfl(i, j)    = w1*sfl_da(i, j, l1ci)    + w2*sfl_da(i, j, l2ci)
-           swa(i, j)    = w1*swa_da(i, j, l1ci)    + w2*swa_da(i, j, l2ci)
-           nsf(i, j)    = w1*nsf_da(i, j, l1ci)    + w2*nsf_da(i, j, l2ci)
-           hmlt(i, j)   = w1*hmlt_da(i, j, l1ci)   + w2*hmlt_da(i, j, l2ci)
-           slp(i, j)    = w1*slp_da(i, j, l1ci)    + w2*slp_da(i, j, l2ci)
-           ficem(i, j)  = w1*ficem_da(i, j, l1ci)  + w2*ficem_da(i, j, l2ci)
-           abswnd(i, j) = w1*abswnd_da(i, j, l1ci) + w2*abswnd_da(i, j, l2ci)
-           atmco2(i, j) = w1*atmco2_da(i, j, l1ci) + w2*atmco2_da(i, j, l2ci)
-           atmbrf(i, j) = w1*atmbrf_da(i, j, l1ci) + w2*atmbrf_da(i, j, l2ci)
+           ustarw(i, j)  = w1*ustarw_da(i, j, l1ci)  + w2*ustarw_da(i, j, l2ci)
+           lip(i, j)     = w1*lip_da(i, j, l1ci)     + w2*lip_da(i, j, l2ci)
+           sop(i, j)     = w1*sop_da(i, j, l1ci)     + w2*sop_da(i, j, l2ci)
+           eva(i, j)     = w1*eva_da(i, j, l1ci)     + w2*eva_da(i, j, l2ci)
+           rnf(i, j)     = w1*rnf_da(i, j, l1ci)     + w2*rnf_da(i, j, l2ci)
+           rfi(i, j)     = w1*rfi_da(i, j, l1ci)     + w2*rfi_da(i, j, l2ci)
+           fmltfz(i, j)  = w1*fmltfz_da(i, j, l1ci)  + w2*fmltfz_da(i, j, l2ci)
+           sfl(i, j)     = w1*sfl_da(i, j, l1ci)     + w2*sfl_da(i, j, l2ci)
+           swa(i, j)     = w1*swa_da(i, j, l1ci)     + w2*swa_da(i, j, l2ci)
+           nsf(i, j)     = w1*nsf_da(i, j, l1ci)     + w2*nsf_da(i, j, l2ci)
+           hmlt(i, j)    = w1*hmlt_da(i, j, l1ci)    + w2*hmlt_da(i, j, l2ci)
+           slp(i, j)     = w1*slp_da(i, j, l1ci)     + w2*slp_da(i, j, l2ci)
+           abswnd(i, j)  = w1*abswnd_da(i, j, l1ci)  + w2*abswnd_da(i, j, l2ci)
+           ficem(i, j)   = w1*ficem_da(i, j, l1ci)   + w2*ficem_da(i, j, l2ci)
+           lamult(i, j)  = w1*lamult_da(i, j, l1ci)  + w2*lamult_da(i, j, l2ci)
+           lasl(i, j)    = w1*lasl_da(i, j, l1ci)    + w2*lasl_da(i, j, l2ci)
+           ustokes(i, j) = w1*ustokes_da(i, j, l1ci) + w2*ustokes_da(i, j, l2ci)
+           vstokes(i, j) = w1*vstokes_da(i, j, l1ci) + w2*vstokes_da(i, j, l2ci)
+           atmco2(i, j)  = w1*atmco2_da(i, j, l1ci)  + w2*atmco2_da(i, j, l2ci)
+           atmbrf(i, j)  = w1*atmbrf_da(i, j, l1ci)  + w2*atmbrf_da(i, j, l2ci)
         enddo
         enddo
         do l = 1, isu(j)
@@ -210,8 +220,12 @@ subroutine getfrc_cesm
       call ncdefvar('nsf_da', 'x y', ndouble, 8)
       call ncdefvar('hmlt_da', 'x y', ndouble, 8)
       call ncdefvar('slp_da', 'x y', ndouble, 8)
-      call ncdefvar('ficem_da', 'x y', ndouble, 8)
       call ncdefvar('abswnd_da', 'x y', ndouble, 8)
+      call ncdefvar('ficem_da', 'x y', ndouble, 8)
+      call ncdefvar('lamult_da', 'x y', ndouble, 8)
+      call ncdefvar('lasl_da', 'x y', ndouble, 8)
+      call ncdefvar('ustokes_da', 'x y', ndouble, 8)
+      call ncdefvar('vstokes_da', 'x y', ndouble, 8)
       call ncdefvar('atmco2_da', 'x y', ndouble, 8)
       call ncdefvar('atmbrf_da', 'x y', ndouble, 8)
       call ncdefvar('ztx_da', 'x y', ndouble, 8)
@@ -242,14 +256,22 @@ subroutine getfrc_cesm
                   ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('slp_da', 'x y', slp_da(1 - nbdy, 1 - nbdy, l2ci), &
                   ip, 1, 1._r8, 0._r8, 8)
+      call ncwrtr('abswnd_da', 'x y', abswnd_da(1 - nbdy, 1 - nbdy, l2ci), &
+                  ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('ficem_da', 'x y', ficem_da(1 - nbdy, 1 - nbdy, l2ci), &
                   ip, 1, 1._r8, 0._r8, 8)
-      call ncwrtr('abswnd_da', 'x y', abswnd_da(1 - nbdy, 1 - nbdy, l2ci), &
+      call ncwrtr('lamult_da', 'x y', lamult_da(1 - nbdy, 1 - nbdy, l2ci), &
+                  ip, 1, 1._r8, 0._r8, 8)
+      call ncwrtr('lasl_da', 'x y', lasl_da(1 - nbdy, 1 - nbdy, l2ci), &
+                  ip, 1, 1._r8, 0._r8, 8)
+      call ncwrtr('ustokes_da', 'x y', ustokes_da(1 - nbdy, 1 - nbdy, l2ci), &
+                  ip, 1, 1._r8, 0._r8, 8)
+      call ncwrtr('vstokes_da', 'x y', vstokes_da(1 - nbdy, 1 - nbdy, l2ci), &
                   ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('atmco2_da', 'x y', atmco2_da(1 - nbdy, 1 - nbdy, l2ci), &
                   ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('atmbrf_da', 'x y', atmbrf_da(1 - nbdy, 1 - nbdy, l2ci), &
-     .            ip, 1, 1._r8, 0._r8, 8)
+                  ip, 1, 1._r8, 0._r8, 8)
       call ncwrtr('ztx_da', 'x y', ztx_da(1 - nbdy, 1 - nbdy, l2ci), &
                   iu, 1, 1._r8, 0._r8, 8)
       call ncwrtr('mty_da', 'x y', mty_da(1 - nbdy, 1 - nbdy, l2ci), &
@@ -277,8 +299,12 @@ subroutine getfrc_cesm
          call chksummsk(nsf, ip, 1, 'nsf')
          call chksummsk(hmlt, ip, 1, 'hmlt')
          call chksummsk(slp, ip, 1, 'slp')
-         call chksummsk(ficem, ip, 1, 'ficem')
          call chksummsk(abswnd, ip, 1, 'abswnd')
+         call chksummsk(ficem, ip, 1, 'ficem')
+         call chksummsk(lamult, ip, 1, 'lamult')
+         call chksummsk(lasl, ip, 1, 'lasl')
+         call chksummsk(ustokes, ip, 1, 'ustokes')
+         call chksummsk(vstokes, ip, 1, 'vstokes')
          call chksummsk(atmco2, ip, 1, 'atmco2')
          call chksummsk(atmbrf, ip, 1, 'atmbrf')
       endif

drivers/nuopc/mod_nuopc_methods.F90

@@ -36,8 +36,9 @@ module mod_nuopc_methods
    use mod_cesm, only: frzpot, mltpot, &
                        swa_da, nsf_da, hmlt_da, lip_da, sop_da, eva_da, &
                        rnf_da, rfi_da, fmltfz_da, sfl_da, ztx_da, mty_da, &
-                       ustarw_da, slp_da, abswnd_da, atmco2_da, atmbrf_da, &
-                       ficem_da, l1ci, l2ci
+                       ustarw_da, slp_da, abswnd_da, ficem_da, lamult_da, &
+                       lasl_da, ustokes_da, vstokes_da, atmco2_da, atmbrf_da, &
+                       l1ci, l2ci
    use mod_utility, only: util1, util2
    use mod_checksum, only: csdiag, chksummsk
    use shr_const_mod, only: SHR_CONST_RHOSW, SHR_CONST_LATICE, SHR_CONST_TKFRZ
@@ -492,10 +493,10 @@ subroutine blom_importflds(fldlist_num, fldlist)
          index_Sw_ustokes  = - 1, &
          index_Sw_vstokes  = - 1, &
          index_Sw_hstokes  = - 1, &
-         index_Sa_pslv     = - 1, &
          index_Faxa_lwdn   = - 1, &
          index_Faxa_snow   = - 1, &
          index_Faxa_rain   = - 1, &
+         index_Sa_pslv     = - 1, &
          index_Sa_co2diag  = - 1, &
          index_Sa_co2prog  = - 1, &
          index_Sa_brfprog  = - 1
@@ -512,7 +513,7 @@ subroutine blom_importflds(fldlist_num, fldlist)
       call getfldindex(fldlist_num, fldlist, 'Foxx_taux', index_Foxx_taux)
       call getfldindex(fldlist_num, fldlist, 'Foxx_tauy', index_Foxx_tauy)
 
-   !$omp parallel do private(i, n)
+   !$omp parallel do private(i, n, afac, utmp, vtmp)
       do j = 1, jjcpl
          do i = 1, ii
             if     (ip(i,j) == 0) then
@@ -526,6 +527,7 @@ subroutine blom_importflds(fldlist_num, fldlist)
             else
                n = (j - 1)*ii + i
                afac = med2mod_areacor(n)
+
                utmp = fldlist(index_Foxx_taux)%dataptr(n)*afac
                vtmp = fldlist(index_Foxx_tauy)%dataptr(n)*afac
                util1(i,j) =   utmp*cosang(i,j) + vtmp*sinang(i,j)
@@ -577,10 +579,10 @@ subroutine blom_importflds(fldlist_num, fldlist)
       call getfldindex(fldlist_num, fldlist, 'Faxa_lwdn', index_Faxa_lwdn)
       call getfldindex(fldlist_num, fldlist, 'Fioi_melth', index_Fioi_melth)
       call getfldindex(fldlist_num, fldlist, 'Sa_pslv', index_Sa_pslv)
-      call getfldindex(fldlist_num, fldlist, 'Si_ifrac', index_Si_ifrac)
       call getfldindex(fldlist_num, fldlist, 'So_duu10n', index_So_duu10n)
+      call getfldindex(fldlist_num, fldlist, 'Si_ifrac', index_Si_ifrac)
 
-   !$omp parallel do private(i, n)
+   !$omp parallel do private(i, n, afac)
       do j = 1, jjcpl
          do i = 1, ii
 
@@ -596,8 +598,8 @@ subroutine blom_importflds(fldlist_num, fldlist)
                nsf_da(i,j,l2ci) = mval
                hmlt_da(i,j,l2ci) = mval
                slp_da(i,j,l2ci) = mval
-               ficem_da(i,j,l2ci) = mval
                abswnd_da(i,j,l2ci) = mval
+               ficem_da(i,j,l2ci) = mval
             elseif (cplmsk(i,j) == 0) then
                lip_da(i,j,l2ci) = 0._r8
                sop_da(i,j,l2ci) = 0._r8
@@ -610,8 +612,8 @@ subroutine blom_importflds(fldlist_num, fldlist)
                nsf_da(i,j,l2ci) = 0._r8
                hmlt_da(i,j,l2ci) = 0._r8
                slp_da(i,j,l2ci) = fval
-               ficem_da(i,j,l2ci) = fval
                abswnd_da(i,j,l2ci) = fval
+               ficem_da(i,j,l2ci) = fval
             else
                n = (j - 1)*ii + i
                afac = med2mod_areacor(n)
@@ -656,12 +658,12 @@ subroutine blom_importflds(fldlist_num, fldlist)
                ! Sea level pressure [kg m-1 s-2].
                slp_da(i,j,l2ci) = fldlist(index_Sa_pslv)%dataptr(n)
 
-               ! Ice fraction [].
-               ficem_da(i,j,l2ci) = fldlist(index_Si_ifrac)%dataptr(n)
-
                ! 10m wind speed [m s-1].
                abswnd_da(i,j,l2ci) = sqrt(fldlist(index_So_duu10n)%dataptr(n))
 
+               ! Ice fraction [].
+               ficem_da(i,j,l2ci) = fldlist(index_Si_ifrac)%dataptr(n)
+
             endif
 
          enddo
@@ -682,8 +684,70 @@ subroutine blom_importflds(fldlist_num, fldlist)
       endif
 
       call fill_global(mval, fval, halo_ps, slp_da(1-nbdy,1-nbdy,l2ci))
-      call fill_global(mval, fval, halo_ps, ficem_da(1-nbdy,1-nbdy,l2ci))
       call fill_global(mval, fval, halo_ps, abswnd_da(1-nbdy,1-nbdy,l2ci))
+      call fill_global(mval, fval, halo_ps, ficem_da(1-nbdy,1-nbdy,l2ci))
+
+      call getfldindex(fldlist_num, fldlist, 'Sw_lamult', index_Sw_lamult)
+      call getfldindex(fldlist_num, fldlist, 'Sw_ustokes', index_Sw_ustokes)
+      call getfldindex(fldlist_num, fldlist, 'Sw_vstokes', index_Sw_vstokes)
+      call getfldindex(fldlist_num, fldlist, 'Sw_hstokes', index_Sw_hstokes)
+
+   !$omp parallel do private(i, n, utmp, vtmp)
+      do j = 1, jjcpl
+         do i = 1, ii
+            if     (ip(i,j) == 0) then
+               util1(i,j) = mval
+               util2(i,j) = mval
+               lamult_da(i,j,l2ci) = mval
+               lasl_da(i,j,l2ci) = mval
+            elseif (cplmsk(i,j) == 0) then
+               util1(i,j) = fval
+               util2(i,j) = fval
+               lamult_da(i,j,l2ci) = fval
+               lasl_da(i,j,l2ci) = fval
+            else
+               n = (j - 1)*ii + i
+
+               utmp = fldlist(index_Sw_ustokes)%dataptr(n)
+               vtmp = fldlist(index_Sw_vstokes)%dataptr(n)
+               util1(i,j) =   utmp*cosang(i,j) + vtmp*sinang(i,j)
+               util2(i,j) = - utmp*sinang(i,j) + vtmp*cosang(i,j)
+
+               ! Langmuir enhancement factor [].
+               lamult_da(i,j,l2ci) = fldlist(index_Sw_lamult)%dataptr(n)
+
+               ! Surface layer averaged Langmuir number [].
+               lasl_da(i,j,l2ci) = fldlist(index_Sw_hstokes)%dataptr(n)
+
+            endif
+         enddo
+      enddo
+   !$omp end parallel do
+
+      call fill_global(mval, fval, halo_pv, util1)
+      call fill_global(mval, fval, halo_pv, util2)
+      call fill_global(mval, fval, halo_ps, lamult_da(1-nbdy,1-nbdy,l2ci))
+      call fill_global(mval, fval, halo_ps, lasl_da(1-nbdy,1-nbdy,l2ci))
+
+      call xctilr(util1, 1,1, 1,1, halo_pv)
+      call xctilr(util2, 1,1, 1,1, halo_pv)
+
+   !$omp parallel do private(l, i)
+      do j = 1, jj
+         do l = 1, isu(j)
+         do i = max(1,ifu(j,l)), min(ii,ilu(j,l))
+            ! x-component of surface Stokes drift [m s-1].
+            ustokes_da(i,j,l2ci) = .5_r8*(util1(i-1,j) + util1(i,j))
+         enddo
+         enddo
+         do l = 1,isv(j)
+         do i = max(1,ifv(j,l)), min(ii,ilv(j,l))
+            ! y-component of surface Stokes drift [m s-1].
+            vstokes_da(i,j,l2ci) = .5_r8*(util2(i,j-1) + util2(i,j))
+         enddo
+         enddo
+      enddo
+   !$omp end parallel do
 
 #ifdef PROGCO2
       call getfldindex(fldlist_num, fldlist, 'Sa_co2prog', index_Sa_co2prog)
@@ -831,8 +895,8 @@ subroutine blom_importflds(fldlist_num, fldlist)
          call chksummsk(nsf_da(1-nbdy,1-nbdy,l2ci),ip,1,'nsf')
          call chksummsk(hmlt_da(1-nbdy,1-nbdy,l2ci),ip,1,'hmlt')
          call chksummsk(slp_da(1-nbdy,1-nbdy,l2ci),ip,1,'slp')
-         call chksummsk(ficem_da(1-nbdy,1-nbdy,l2ci),ip,1,'ficem')
          call chksummsk(abswnd_da(1-nbdy,1-nbdy,l2ci),ip,1,'abswnd')
+         call chksummsk(ficem_da(1-nbdy,1-nbdy,l2ci),ip,1,'ficem')
          call chksummsk(atmco2_da(1-nbdy,1-nbdy,l2ci),ip,1,'atmco2')
          call chksummsk(atmbrf_da(1-nbdy,1-nbdy,l2ci),ip,1,'atmbrf')
       endif

phy/mod_difest.F

@@ -214,6 +214,8 @@ subroutine init_difest
 c --- ------------------------------------------------------------------
 c --- Initialize CVmix variables.
 c --- ------------------------------------------------------------------
+c
+      integer :: i,j,l
 c
 c -- ------- Background diapycnal mixing.
 c         The Bryan-Lewis parameterization is based on the following:
@@ -283,6 +285,15 @@ subroutine init_difest
 c    .          lnoDGat1=.true.                 ,
 c    .          CVMix_kpp_params_user=KPP_params )
 c
+c$OMP PARALLEL DO PRIVATE(l,i)
+      do j=1,jj
+        do l=1,isp(j)
+        do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
+          OBLdepth(i,j)=10.
+        enddo
+        enddo
+      enddo
+c$OMP END PARALLEL DO
 c
       end subroutine init_difest
 c

phy/mod_forcing.F90

@@ -106,6 +106,10 @@ module mod_forcing
       ustarw, &       ! Friction velocity for open water [m s-1].
       slp, &          ! Sea-level pressure [kg m-1 s-2].
       abswnd, &       ! Wind speed at measurement height (zu) [m s-1].
+      lamult, &       ! Langmuir enhancement factor [].
+      lasl, &         ! Surface layer averaged Langmuir number [].
+      ustokes, &      ! u-component of surface Stokes drift [m s-1].
+      vstokes, &      ! v-component of surface Stokes drift [m s-1].
       atmco2, &       ! Atmospheric CO2 concentration [ppm].
       flxco2, &       ! Air-sea CO2 flux [kg m-2 s-1].
       flxdms, &       ! Sea-air DMS flux [kg m-2 s-1].
@@ -138,7 +142,8 @@ module mod_forcing
              sref, tflxap, sflxap, tflxdi, sflxdi, nflxdi, &
              sstclm, ricclm, sssclm, prfac, eiacc, pracc, &
              swa, nsf, hmltfz, lip, sop, eva, rnf, rfi, fmltfz, sfl, ztx, mty, &
-             ustarw, slp, abswnd, atmco2, flxco2, flxdms, flxbrf, atmbrf, &
+             ustarw, slp, abswnd, lamult, lasl, ustokes, vstokes, &
+             atmco2, flxco2, flxdms, flxbrf, atmbrf, &
              surflx, surrlx, sswflx, salflx, brnflx, salrlx, taux, tauy, &
              ustar, ustarb, ustar3, buoyfl, t_sw_nonloc, &
              inivar_forcing, fwbbal
@@ -172,6 +177,10 @@ subroutine inivar_forcing
             ustarw(i, j) = spval
             slp(i, j) = spval
             abswnd(i, j) = spval
+            lamult(i, j) = spval
+            lasl(i, j) = spval
+            ustokes(i, j) = spval
+            vstokes(i, j) = spval
             atmco2(i, j) = spval
             flxco2(i, j) = spval
             flxdms(i, j) = spval