Merge pull request #454 from shansun6/iceflx_ttice_20200520

initialize ice flux and add "tiice" array

physics/GFS_PBL_generic.F90

@@ -330,8 +330,8 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
         dqdt, dusfc_cpl, dvsfc_cpl, dtsfc_cpl,                                                                                 &
         dqsfc_cpl, dusfci_cpl, dvsfci_cpl, dtsfci_cpl, dqsfci_cpl, dusfc_diag, dvsfc_diag, dtsfc_diag, dqsfc_diag,             &
         dusfci_diag, dvsfci_diag, dtsfci_diag, dqsfci_diag, dt3dt, du3dt_PBL, du3dt_OGWD, dv3dt_PBL, dv3dt_OGWD, dq3dt,        &
-        dq3dt_ozone, rd, cp,fvirt, hvap, t1, q1, prsl, hflx, ushfsfci, oceanfrac, fice, dusfc_cice, dvsfc_cice, dtsfc_cice,    &
-        dqsfc_cice, wet, dry, icy, wind, stress_ocn, hflx_ocn, evap_ocn, ugrs1, vgrs1, dkt_cpl, dkt, hffac, hefac,             &
+        dq3dt_ozone, rd, cp, fvirt, hvap, t1, q1, prsl, hflx, ushfsfci, oceanfrac, flag_cice, dusfc_cice, dvsfc_cice,          & 
+        dtsfc_cice, dqsfc_cice, wet, dry, icy, wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1, dkt_cpl, dkt, hffac, hefac, &
         errmsg, errflg)
 
       use machine,                only : kind_phys
@@ -346,13 +346,14 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       integer, intent(in) :: imp_physics_zhao_carr, imp_physics_mg, imp_physics_fer_hires
       logical, intent(in) :: ltaerosol, cplflx, cplchm, lssav, ldiag3d, lsidea
       logical, intent(in) :: hybedmf, do_shoc, satmedmf, shinhong, do_ysu
+      logical, dimension(:), intent(in) :: flag_cice
 
       real(kind=kind_phys), intent(in) :: dtf
       real(kind=kind_phys), intent(in) :: rd, cp, fvirt, hvap
-      real(kind=kind_phys), dimension(:), intent(in) :: t1, q1, hflx, oceanfrac, fice
+      real(kind=kind_phys), dimension(:), intent(in) :: t1, q1, hflx, oceanfrac
       real(kind=kind_phys), dimension(:,:), intent(in) :: prsl
       real(kind=kind_phys), dimension(:), intent(in) :: dusfc_cice, dvsfc_cice, dtsfc_cice, dqsfc_cice, &
-          wind, stress_ocn, hflx_ocn, evap_ocn, ugrs1, vgrs1
+          wind, stress_wat, hflx_wat, evap_wat, ugrs1, vgrs1
       real(kind=kind_phys), dimension(im, levs, nvdiff), intent(in) :: dvdftra
       real(kind=kind_phys), dimension(im), intent(in) :: dusfc1, dvsfc1, dtsfc1, dqsfc1, xmu
       real(kind=kind_phys), dimension(im, levs), intent(in) :: dudt, dvdt, dtdt, htrsw, htrlw
@@ -382,7 +383,6 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       real(kind=kind_phys), parameter :: zero  = 0.0d0
       real(kind=kind_phys), parameter :: one   = 1.0d0
       real(kind=kind_phys), parameter :: huge  = 9.9692099683868690E36 ! NetCDF float FillValue, same as in GFS_typedefs.F90
-      real(kind=kind_phys), parameter :: epsln = 1.0d-10 ! same as in GFS_physics_driver.F90
       integer :: i, k, kk, k1, n
       real(kind=kind_phys) :: tem, tem1, rho
 
@@ -550,24 +550,31 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac,
       if (cplflx) then
         do i=1,im
           if (oceanfrac(i) > zero) then ! Ocean only, NO LAKES
-            if (fice(i) > one - epsln) then ! no open water, use results from CICE
-              dusfci_cpl(i) = dusfc_cice(i)
-              dvsfci_cpl(i) = dvsfc_cice(i)
-              dtsfci_cpl(i) = dtsfc_cice(i)
-              dqsfci_cpl(i) = dqsfc_cice(i)
+            if ( .not. wet(i)) then ! no open water
+              if (flag_cice(i)) then !use results from CICE
+                dusfci_cpl(i) = dusfc_cice(i)
+                dvsfci_cpl(i) = dvsfc_cice(i)
+                dtsfci_cpl(i) = dtsfc_cice(i)
+                dqsfci_cpl(i) = dqsfc_cice(i)
+              else !use PBL fluxes when CICE fluxes is unavailable
+                dusfci_cpl(i) = dusfc1(i)
+                dvsfci_cpl(i) = dvsfc1(i)
+                dtsfci_cpl(i) = dtsfc1(i)
+                dqsfci_cpl(i) = dqsfc1(i)
+              end if
             elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point
               tem1 = max(q1(i), 1.e-8)
               rho = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1))
               if (wind(i) > zero) then
-                tem = - rho * stress_ocn(i) / wind(i)
+                tem = - rho * stress_wat(i) / wind(i)
                 dusfci_cpl(i) = tem * ugrs1(i)   ! U-momentum flux
                 dvsfci_cpl(i) = tem * vgrs1(i)   ! V-momentum flux
               else
                 dusfci_cpl(i) = zero
                 dvsfci_cpl(i) = zero
               endif
-              dtsfci_cpl(i) = cp   * rho * hflx_ocn(i) ! sensible heat flux over open ocean
-              dqsfci_cpl(i) = hvap * rho * evap_ocn(i) ! latent heat flux over open ocean
+              dtsfci_cpl(i) = cp   * rho * hflx_wat(i) ! sensible heat flux over open ocean
+              dqsfci_cpl(i) = hvap * rho * evap_wat(i) ! latent heat flux over open ocean
             else                                       ! use results from PBL scheme for 100% open ocean
               dusfci_cpl(i) = dusfc1(i)
               dvsfci_cpl(i) = dvsfc1(i)

physics/GFS_PBL_generic.meta

@@ -1186,13 +1186,12 @@
   kind = kind_phys
   intent = in
   optional = F
-[fice]
-  standard_name = sea_ice_concentration
-  long_name = ice fraction over open water
-  units = frac
+[flag_cice]
+  standard_name = flag_for_cice
+  long_name = flag for cice
+  units = flag
   dimensions = (horizontal_dimension)
-  type = real
-  kind = kind_phys
+  type = logical
   intent = in
   optional = F
 [dusfc_cice]
@@ -1264,7 +1263,7 @@
   kind = kind_phys
   intent = in
   optional = F
-[stress_ocn]
+[stress_wat]
   standard_name = surface_wind_stress_over_ocean
   long_name = surface wind stress over ocean
   units = m2 s-2
@@ -1273,7 +1272,7 @@
   kind = kind_phys
   intent = in
   optional = F
-[hflx_ocn]
+[hflx_wat]
   standard_name = kinematic_surface_upward_sensible_heat_flux_over_ocean
   long_name = kinematic surface upward sensible heat flux over ocean
   units = K m s-1
@@ -1282,7 +1281,7 @@
   kind = kind_phys
   intent = in
   optional = F
-[evap_ocn]
+[evap_wat]
   standard_name = kinematic_surface_upward_latent_heat_flux_over_ocean
   long_name = kinematic surface upward latent heat flux over ocean
   units = kg kg-1 m s-1

physics/GFS_suite_interstitial.F90

@@ -160,7 +160,7 @@ end subroutine GFS_suite_interstitial_2_finalize
     subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplflx, flag_cice, shal_cnv, old_monin, mstrat,       &
       do_shoc, frac_grid, imfshalcnv, dtf, xcosz, adjsfcdsw, adjsfcdlw, cice, pgr, ulwsfc_cice, lwhd, htrsw, htrlw, xmu, ctei_rm,    &
       work1, work2, prsi, tgrs, prsl, qgrs_water_vapor, qgrs_cloud_water, cp, hvap, prslk, suntim, adjsfculw, adjsfculw_lnd,         &
-      adjsfculw_ice, adjsfculw_ocn, dlwsfc, ulwsfc, psmean, dt3dt_lw, dt3dt_sw, dt3dt_pbl, dt3dt_dcnv, dt3dt_scnv, dt3dt_mp,         &
+      adjsfculw_ice, adjsfculw_wat, dlwsfc, ulwsfc, psmean, dt3dt_lw, dt3dt_sw, dt3dt_pbl, dt3dt_dcnv, dt3dt_scnv, dt3dt_mp,         &
       ctei_rml, ctei_r, kinver, dry, icy, wet, frland, huge, errmsg, errflg)
 
       implicit none
@@ -181,7 +181,7 @@ subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplfl
 
       integer,              intent(inout), dimension(im) :: kinver
       real(kind=kind_phys), intent(inout), dimension(im) :: suntim, dlwsfc, ulwsfc, psmean, ctei_rml, ctei_r
-      real(kind=kind_phys), intent(in   ), dimension(im) :: adjsfculw_lnd, adjsfculw_ice, adjsfculw_ocn
+      real(kind=kind_phys), intent(in   ), dimension(im) :: adjsfculw_lnd, adjsfculw_ice, adjsfculw_wat
       real(kind=kind_phys), intent(  out), dimension(im) :: adjsfculw
 
       ! These arrays are only allocated if ldiag3d is .true.
@@ -232,11 +232,11 @@ subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplfl
             if (flag_cice(i)) then
               adjsfculw(i) = adjsfculw_lnd(i) * frland(i)               &
                            + ulwsfc_cice(i)   * tem                     &
-                           + adjsfculw_ocn(i) * (one - frland(i) - tem)
+                           + adjsfculw_wat(i) * (one - frland(i) - tem)
             else
               adjsfculw(i) = adjsfculw_lnd(i) * frland(i)               &
                            + adjsfculw_ice(i) * tem                     &
-                           + adjsfculw_ocn(i) * (one - frland(i) - tem)
+                           + adjsfculw_wat(i) * (one - frland(i) - tem)
             endif
           enddo
         else
@@ -246,20 +246,20 @@ subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplfl
             elseif (icy(i)) then                 ! ice (and water)
               tem = one - cice(i)
               if (flag_cice(i)) then
-                if (wet(i) .and. adjsfculw_ocn(i) /= huge) then
-                  adjsfculw(i) = ulwsfc_cice(i)*cice(i) + adjsfculw_ocn(i)*tem
+                if (wet(i) .and. adjsfculw_wat(i) /= huge) then
+                  adjsfculw(i) = ulwsfc_cice(i)*cice(i) + adjsfculw_wat(i)*tem
                 else
                   adjsfculw(i) = ulwsfc_cice(i)
                 endif
               else
-                if (wet(i) .and. adjsfculw_ocn(i) /= huge) then
-                  adjsfculw(i) = adjsfculw_ice(i)*cice(i) + adjsfculw_ocn(i)*tem
+                if (wet(i) .and. adjsfculw_wat(i) /= huge) then
+                  adjsfculw(i) = adjsfculw_ice(i)*cice(i) + adjsfculw_wat(i)*tem
                 else
                   adjsfculw(i) = adjsfculw_ice(i)
                 endif
               endif
             else                                 ! all water
-              adjsfculw(i) = adjsfculw_ocn(i)
+              adjsfculw(i) = adjsfculw_wat(i)
             endif
           enddo
         endif

physics/GFS_suite_interstitial.meta

@@ -604,7 +604,7 @@
   kind = kind_phys
   intent = in
   optional = F
-[adjsfculw_ocn]
+[adjsfculw_wat]
   standard_name = surface_upwelling_longwave_flux_over_ocean_interstitial
   long_name = surface upwelling longwave flux at current time over ocean (temporary use as interstitial)
   units = W m-2