Merge pull request mom-ocean#28 from gustavo-marques/fix_latent_heat_…

…flux

Pass latent heat flux via ice/ocean boundary type

config_src/coupled_driver/MOM_surface_forcing.F90

@@ -152,6 +152,7 @@ module MOM_surface_forcing
   real, pointer, dimension(:,:) :: sw_flux_vis_dif      =>NULL() ! diffuse visible sw radiation (W/m2)
   real, pointer, dimension(:,:) :: sw_flux_nir_dir      =>NULL() ! direct Near InfraRed sw radiation (W/m2)
   real, pointer, dimension(:,:) :: sw_flux_nir_dif      =>NULL() ! diffuse Near InfraRed sw radiation (W/m2)
+  real, pointer, dimension(:,:) :: latent               =>NULL() ! latent heat flux (W/m2)
   real, pointer, dimension(:,:) :: lprec                =>NULL() ! mass flux of liquid precip (kg/m2/s)
   real, pointer, dimension(:,:) :: fprec                =>NULL() ! mass flux of frozen precip (kg/m2/s)
   real, pointer, dimension(:,:) :: runoff               =>NULL() ! mass flux of liquid runoff (kg/m2/s)
@@ -480,17 +481,21 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state,
       fluxes%sens(i,j) = - IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j)
 
     fluxes%latent(i,j) = 0.0
-    if (ASSOCIATED(IOB%fprec)) then
-      fluxes%latent(i,j)            = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
-      fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
-    endif
-    if (ASSOCIATED(IOB%calving)) then
-      fluxes%latent(i,j)              = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion
-      fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion
-    endif
-    if (ASSOCIATED(IOB%q_flux)) then
-      fluxes%latent(i,j)           = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
-      fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+    if (ASSOCIATED(IOB%latent)) then
+      fluxes%latent(i,j) = IOB%latent(i-i0,j-j0)
+    else
+      if (ASSOCIATED(IOB%fprec)) then
+        fluxes%latent(i,j)            = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+        fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion
+      endif
+      if (ASSOCIATED(IOB%calving)) then
+        fluxes%latent(i,j)              = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion
+        fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion
+      endif
+      if (ASSOCIATED(IOB%q_flux)) then
+        fluxes%latent(i,j)           = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+        fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor
+      endif
     endif
 
     fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j)

config_src/coupled_driver/ocean_model_MOM.F90

@@ -188,7 +188,7 @@ module ocean_model_mod
 
 !> ocean_model_init initializes the ocean model, including registering fields
 !! for restarts and reading restart files if appropriate.
-subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file)
+subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
   type(ocean_public_type), target, &
                        intent(inout) :: Ocean_sfc !< A structure containing various
                                 !! publicly visible ocean surface properties after initialization,
@@ -204,7 +204,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
                                               !! in the calculation of additional gas or other
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
-  character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read
 
 !   This subroutine initializes both the ocean state and the ocean surface type.
 ! Because of the way that indicies and domains are handled, Ocean_sfc must have
@@ -241,7 +240,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   OS%Time = Time_in
   call initialize_MOM(OS%Time, param_file, OS%dirs, OS%MOM_CSp, Time_in, &
-      offline_tracer_mode=offline_tracer_mode, input_restart_file=input_restart_file)
+      offline_tracer_mode=offline_tracer_mode)
   OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV
   OS%C_p = OS%MOM_CSp%tv%C_p
   OS%fluxes%C_p = OS%MOM_CSp%tv%C_p

config_src/mct_driver/ocn_comp_mct.F90

@@ -596,6 +596,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )
   allocate(glb%ice_ocean_boundary%v_flux(isc:iec,jsc:jec));          glb%ice_ocean_boundary%v_flux(:,:) = 0.0
   allocate(glb%ice_ocean_boundary%t_flux(isc:iec,jsc:jec));          glb%ice_ocean_boundary%t_flux(:,:) = 0.0
   allocate(glb%ice_ocean_boundary%q_flux(isc:iec,jsc:jec));          glb%ice_ocean_boundary%q_flux(:,:) = 0.0
+  allocate(glb%ice_ocean_boundary%latent(isc:iec,jsc:jec));          glb%ice_ocean_boundary%latent(:,:) = 0.0
   allocate(glb%ice_ocean_boundary%salt_flux(isc:iec,jsc:jec));       glb%ice_ocean_boundary%salt_flux(:,:) = 0.0
   allocate(glb%ice_ocean_boundary%lw_flux(isc:iec,jsc:jec));         glb%ice_ocean_boundary%lw_flux(:,:) = 0.0
   allocate(glb%ice_ocean_boundary%sw_flux_vis_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dir(:,:) = 0.0
@@ -1130,54 +1131,55 @@ subroutine ocn_export(ind, ocn_public, grid, o2x)
   ! Update halo of ssh so we can calculate gradients
   call pass_var(ssh, grid%domain)
 
-
   ! d/dx ssh
   n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
     n = n+1
     ! This is a simple second-order difference
-    ! o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j)
+    o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
-    slp_L = (ssh(i,j) - ssh(i-1,j)) * grid%mask2dCu(I-1,j)
+!    slp_L = (ssh(I,j) - ssh(I-1,j)) * grid%mask2dCu(I-1,j)
     !if (grid%mask2dCu(I-1,j)==0.) slp_L = 0.
-    slp_R = (ssh(i+1,j) - ssh(i,j)) * grid%mask2dCu(I,j)
+!    slp_R = (ssh(I+1,j) - ssh(I,j)) * grid%mask2dCu(I,j)
     !if (grid%mask2dCu(I,j)==0.) slp_R = 0.
-    slp_C = 0.5 * (slp_L + slp_R)
-    if ( (slp_L * slp_R) > 0.0 ) then
+!    slp_C = 0.5 * (slp_L + slp_R)
+!    if ( (slp_L * slp_R) > 0.0 ) then
       ! This limits the slope so that the edge values are bounded by the
       ! two cell averages spanning the edge.
-      u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
-      u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
-      slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
-    else
+!      u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
+!      u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) )
+!      slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
+!    else
       ! Extrema in the mean values require a PCM reconstruction avoid generating
       ! larger extreme values.
-      slope = 0.0
-    end if
-    o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
+!      slope = 0.0
+!    end if
+!    o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
   end do; end do
 
   ! d/dy ssh
+  n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
+    n = n+1
     ! This is a simple second-order difference
-    ! o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j)
+    o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
-    slp_L = ssh(i,j) - ssh(i,j-1)
-    slp_R = ssh(i,j+1) - ssh(i,j)
-    slp_L=0.; slp_R=0.
-    slp_C = 0.5 * (slp_L + slp_R)
-    if ((slp_L * slp_R) > 0.0) then
+!    slp_L = ssh(i,J) - ssh(i,J-1) * grid%mask2dCv(i,J-1)
+!    slp_R = ssh(i,J+1) - ssh(i,J) * grid%mask2dCv(i,J)
+!    slp_C = 0.5 * (slp_L + slp_R)
+!    write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R
+!    if ((slp_L * slp_R) > 0.0) then
       ! This limits the slope so that the edge values are bounded by the
       ! two cell averages spanning the edge.
-      u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
-      u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
-      slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
-    else
+!      u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
+!      u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) )
+!      slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C )
+!    else
       ! Extrema in the mean values require a PCM reconstruction avoid generating
       ! larger extreme values.
-      slope = 0.0
-    end if
-    o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
+!      slope = 0.0
+!    end if
+!    o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
   end do; end do
 
 end subroutine ocn_export
@@ -1258,6 +1260,8 @@ subroutine fill_ice_ocean_bnd(ice_ocean_boundary, grid, x2o_o, ind, &
       ice_ocean_boundary%v_flux(ig,jg) = x2o_o(ind%x2o_Foxx_tauy,k)
       ! sensible heat flux (W/m2)
       ice_ocean_boundary%t_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_sen,k)
+      ! latent heat flux (W/m^2)
+      ice_ocean_boundary%latent(ig,jg) = x2o_o(ind%x2o_Foxx_lat,k)
       ! salt flux
       ice_ocean_boundary%salt_flux(ig,jg) = -x2o_o(ind%x2o_Fioi_salt,k)
       ! heat content from frozen runoff
@@ -1309,6 +1313,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
   type(ESMF_timeInterval) :: ocn_cpl_interval !< The length of one ocean coupling interval
   integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc
   logical :: write_restart_at_eod      !< Controls if restart files must be written
+  logical :: debug=.false.
   type(time_type) :: time_start        !< Start of coupled time interval to pass to MOM6
   type(time_type) :: coupling_timestep !< Coupled time interval to pass to MOM6
   character(len=128) :: err_msg        !< Error message
@@ -1367,6 +1372,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
   call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, &
                           time_start, coupling_timestep)
 
+  ! return export state to driver
+  call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr)
+
   !--- write out intermediate restart file when needed.
   ! Check alarms for flag to write restart at end of day
   write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock)