Merge pull request #77 from gustavo-marques/add_melt_potential_mct

Add melt/freeze potential via MCT cap

config_src/mct_driver/MOM_ocean_model.F90

@@ -139,7 +139,8 @@ module MOM_ocean_model
                          !! i.e. dzt(1) + eta_t + patm/rho0/grav (m)
     frazil =>NULL(), &   !< Accumulated heating (in Joules/m^2) from frazil
                          !! formation in the ocean.
-    area => NULL()       !< cell area of the ocean surface, in m2.
+    melt_potential => NULL(), & !< Accumulated heat used to melt sea ice (in W/m^2)
+    area => NULL()              !< cell area of the ocean surface, in m2.
   type(coupler_2d_bc_type) :: fields    !< A structure that may contain an
                                         !! array of named tracer-related fields.
   integer                  :: avg_kount !< Used for accumulating averages of this type.
@@ -337,8 +338,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i
 
   !   Consider using a run-time flag to determine whether to do the diagnostic
   ! vertical integrals, since the related 3-d sums are not negligible in cost.
-  call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, &
-                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn)
+  call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, do_integrals=.true., &
+                              gas_fields_ocn=gas_fields_ocn, use_meltpot=.true.)
 
   call surface_forcing_init(Time_in, OS%grid, param_file, OS%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)
@@ -706,20 +707,22 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap,
   endif
   call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec)
 
-  allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), &
-             Ocean_sfc%s_surf (isc:iec,jsc:jec), &
-             Ocean_sfc%u_surf (isc:iec,jsc:jec), &
-             Ocean_sfc%v_surf (isc:iec,jsc:jec), &
-             Ocean_sfc%sea_lev(isc:iec,jsc:jec), &
-             Ocean_sfc%area   (isc:iec,jsc:jec), &
-             Ocean_sfc%frazil (isc:iec,jsc:jec))
+  allocate (Ocean_sfc%t_surf (isc:iec,jsc:jec), &
+            Ocean_sfc%s_surf (isc:iec,jsc:jec), &
+            Ocean_sfc%u_surf (isc:iec,jsc:jec), &
+            Ocean_sfc%v_surf (isc:iec,jsc:jec), &
+            Ocean_sfc%sea_lev(isc:iec,jsc:jec), &
+            Ocean_sfc%area   (isc:iec,jsc:jec), &
+            Ocean_sfc%melt_potential(isc:iec,jsc:jec), &
+            Ocean_sfc%frazil (isc:iec,jsc:jec))
 
   Ocean_sfc%t_surf  = 0.0  ! time averaged sst (Kelvin) passed to atmosphere/ice model
   Ocean_sfc%s_surf  = 0.0  ! time averaged sss (psu) passed to atmosphere/ice models
   Ocean_sfc%u_surf  = 0.0  ! time averaged u-current (m/sec) passed to atmosphere/ice models
   Ocean_sfc%v_surf  = 0.0  ! time averaged v-current (m/sec)  passed to atmosphere/ice models
   Ocean_sfc%sea_lev = 0.0  ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav
   Ocean_sfc%frazil  = 0.0  ! time accumulated frazil (J/m^2) passed to ice model
+  Ocean_sfc%melt_potential  = 0.0  ! time accumulated melt potential (J/m^2) passed to ice model
   Ocean_sfc%area    = 0.0
   Ocean_sfc%axes    = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics
 
@@ -783,11 +786,13 @@ subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, patm, press_to_z)
 
   do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd
     Ocean_sfc%sea_lev(i,j) = state%sea_lev(i+i0,j+j0)
+    Ocean_sfc%area(i,j)   =  G%areaT(i+i0,j+j0)
     if (present(patm)) &
       Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + patm(i,j) * press_to_z
-      if (associated(state%frazil)) &
+    if (associated(state%frazil)) &
       Ocean_sfc%frazil(i,j) = state%frazil(i+i0,j+j0)
-    Ocean_sfc%area(i,j)   =  G%areaT(i+i0,j+j0)
+    if (allocated(state%melt_potential)) &
+      Ocean_sfc%melt_potential(i,j) = state%melt_potential(i+i0,j+j0)
   enddo ; enddo
 
   if (Ocean_sfc%stagger == AGRID) then
@@ -1012,6 +1017,7 @@ subroutine ocean_public_type_chksum(id, timestep, ocn)
   write(outunit,100) 'ocean%v_surf   ',mpp_chksum(ocn%v_surf )
   write(outunit,100) 'ocean%sea_lev  ',mpp_chksum(ocn%sea_lev)
   write(outunit,100) 'ocean%frazil   ',mpp_chksum(ocn%frazil )
+  write(outunit,100) 'ocean%melt_potential  ',mpp_chksum(ocn%melt_potential)
 
   call coupler_type_write_chksums(ocn%fields, outunit, 'ocean%')
 100 FORMAT("   CHECKSUM::",A20," = ",Z20)

config_src/mct_driver/ocn_cap_methods.F90

@@ -141,18 +141,26 @@ end subroutine ocn_import
 !=======================================================================
 
   !> Maps outgoing ocean data to MCT attribute vector real array
-  subroutine ocn_export(ind, ocn_public, grid, o2x)
+  subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
     type(cpl_indices_type),  intent(inout) :: ind        !< Structure with coupler indices and vectors
     type(ocean_public_type), intent(in)    :: ocn_public !< Ocean surface state
     type(ocean_grid_type),   intent(in)    :: grid       !< Ocean model grid
     real(kind=8),            intent(inout) :: o2x(:,:)   !< MCT outgoing bugger
+    real(kind=8), intent(in)               :: dt_int     !< Amount of time over which to advance the
+                                                         !! ocean (ocean_coupling_time_step), in sec
+    integer, intent(in)                    :: ncouple_per_day !< Number of ocean coupling calls per day
 
     ! Local variables
     real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo
     integer :: i, j, n, ig, jg  !< Grid indices
     real    :: slp_L, slp_R, slp_C, slope, u_min, u_max
+    real :: I_time_int  !< The inverse of coupling time interval in s-1.
+
     !-----------------------------------------------------------------------
 
+    ! Use Adcroft's rule of reciprocals; it does the right thing here.
+    I_time_int = 0.0 ; if (dt_int > 0.0) I_time_int = 1.0 / dt_int
+
     ! Copy from ocn_public to o2x. ocn_public uses global indexing with no halos.
     ! The mask comes from "grid" that uses the usual MOM domain that has halos
     ! and does not use global indexing.
@@ -168,6 +176,16 @@ subroutine ocn_export(ind, ocn_public, grid, o2x)
         o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
         o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j)
         o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j)
+        ! ocean melt and freeze potential (o2x_Fioo_q), W m-2
+        if (ocn_public%frazil(ig,jg) > 0.0) then
+          ! Frazil: change from J/m^2 to W/m^2
+          o2x(ind%o2x_Fioo_q, n) = ocn_public%frazil(ig,jg) * grid%mask2dT(i,j) * I_time_int
+        else
+          ! Melt_potential: change from J/m^2 to W/m^2
+          o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * I_time_int !* ncouple_per_day
+          ! make sure Melt_potential is always <= 0
+          if (o2x(ind%o2x_Fioo_q, n) > 0.0) o2x(ind%o2x_Fioo_q, n) = 0.0
+        endif
         ! Make a copy of ssh in order to do a halo update. We use the usual MOM domain
         ! in order to update halos. i.e. does not use global indexing.
         ssh(i,j) = ocn_public%sea_lev(ig,jg)

config_src/mct_driver/ocn_comp_mct.F90

@@ -373,7 +373,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )
   ! end if
 
   if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_export"
-  call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr)
+  call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr, mom_cpl_dt, ncouple_per_day)
 
   call t_stopf('MOM_mct_init')
 
@@ -423,6 +423,10 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
   integer            :: shrlogunit ! original log file unit
   integer            :: shrloglev  ! original log level
   logical, save      :: firstCall = .true.
+  real (kind=8), parameter  ::  seconds_in_day = 86400.0 !< number of seconds in one day
+  integer                   :: ocn_cpl_dt   !< one ocn coupling interval in seconds. (to be received from cesm)
+  real (kind=8)             :: mom_cpl_dt   !< one ocn coupling interval in seconds. (internal)
+  integer                   :: ncouple_per_day !< number of ocean coupled call in one day (non-dim)
 
   ! reset shr logging to ocn log file:
   if (is_root_pe()) then
@@ -441,6 +445,10 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
   call ESMF_TimeIntervalGet(ocn_cpl_interval, yy=year, mm=month, d=day, s=seconds, sn=seconds_n, sd=seconds_d, rc=rc)
   coupling_timestep = set_time(seconds, days=day, err_msg=err_msg)
 
+  call seq_timemgr_EClockGetData(EClock, dtime=ocn_cpl_dt)
+  ncouple_per_day = seconds_in_day / ocn_cpl_dt
+  mom_cpl_dt = seconds_in_day / ncouple_per_day
+
   ! The following if-block is to correct monthly mean outputs:
   ! With this change, MOM6 starts at the same date as the other components, and runs for the same
   ! duration as other components, unlike POP, which would have one missing interval due to ocean
@@ -502,7 +510,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)
   call update_ocean_model(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)
+  call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr, mom_cpl_dt, ncouple_per_day)
 
   !--- write out intermediate restart file when needed.
   ! Check alarms for flag to write restart at end of day
@@ -806,6 +814,5 @@ end subroutine ocean_model_init_sfc
 !! Boundary layer depth
 !! CO2
 !! DMS
-!! o2x_Fioo_q       !< Heat flux?
 
 end module ocn_comp_mct

config_src/mct_driver/ocn_cpl_indices.F90

@@ -16,7 +16,7 @@ module ocn_cpl_indices
     integer :: o2x_So_dhdx       !< Zonal slope in the sea surface height
     integer :: o2x_So_dhdy       !< Meridional lope in the sea surface height
     integer :: o2x_So_bldepth    !< Boundary layer depth (m)
-    integer :: o2x_Fioo_q        !< Heat flux?
+    integer :: o2x_Fioo_q        !< Ocean melt and freeze potential (W/m2)
     integer :: o2x_Faoo_fco2_ocn !< CO2 flux
     integer :: o2x_Faoo_fdms_ocn !< DMS flux
 

src/core/MOM.F90

@@ -5,6 +5,7 @@ module MOM
 
 ! Infrastructure modules
 use MOM_debugging,            only : MOM_debugging_init, hchksum, uvchksum
+use MOM_debugging,            only : check_redundant
 use MOM_checksum_packages,    only : MOM_thermo_chksum, MOM_state_chksum
 use MOM_checksum_packages,    only : MOM_accel_chksum, MOM_surface_chksum
 use MOM_cpu_clock,            only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
@@ -71,8 +72,7 @@ module MOM
 use MOM_dynamics_unsplit_RK2,  only : initialize_dyn_unsplit_RK2, end_dyn_unsplit_RK2
 use MOM_dynamics_unsplit_RK2,  only : MOM_dyn_unsplit_RK2_CS
 use MOM_dyn_horgrid,           only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
-use MOM_debugging,             only : check_redundant
-use MOM_EOS,                   only : EOS_init, calculate_density
+use MOM_EOS,                   only : EOS_init, calculate_density, calculate_TFreeze
 use MOM_fixed_initialization,  only : MOM_initialize_fixed
 use MOM_grid,                  only : ocean_grid_type, set_first_direction
 use MOM_grid,                  only : MOM_grid_init, MOM_grid_end
@@ -193,7 +193,6 @@ module MOM
                     !! bottom drag viscosities, and related fields
   type(MEKE_type), pointer :: MEKE => NULL() !<  structure containing fields
                     !! related to the Mesoscale Eddy Kinetic Energy
-
   logical :: adiabatic !< If true, there are no diapycnal mass fluxes, and no calls
                     !! to routines to calculate or apply diapycnal fluxes.
   logical :: use_legacy_diabatic_driver!< If true (default), use the a legacy version of the diabatic
@@ -545,6 +544,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, &
 
   if (therm_reset) then
     CS%time_in_thermo_cycle = 0.0
+    if (allocated(sfc_state%melt_potential)) sfc_state%melt_potential(:,:)  = 0.0
     if (associated(CS%tv%frazil))        CS%tv%frazil(:,:)        = 0.0
     if (associated(CS%tv%salt_deficit))  CS%tv%salt_deficit(:,:)  = 0.0
     if (associated(CS%tv%TempxPmE))      CS%tv%TempxPmE(:,:)      = 0.0
@@ -2660,20 +2660,20 @@ subroutine extract_surface_state(CS, sfc_state)
 
   ! local
   real :: hu, hv
-  type(ocean_grid_type), pointer :: G => NULL() ! pointer to a structure containing
-                                      ! metrics and related information
+  type(ocean_grid_type), pointer :: G => NULL() !< pointer to a structure containing
+                                      !! metrics and related information
   type(verticalGrid_type), pointer :: GV => NULL()
   real, dimension(:,:,:), pointer :: &
-    u => NULL(), & ! u : zonal velocity component (m/s)
-    v => NULL(), & ! v : meridional velocity component (m/s)
-    h => NULL()    ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss))
-  real :: depth(SZI_(CS%G))           ! distance from the surface (meter)
-  real :: depth_ml                    ! depth over which to average to
-                                      ! determine mixed layer properties (meter)
-  real :: dh                          ! thickness of a layer within mixed layer (meter)
-  real :: mass                        ! mass per unit area of a layer (kg/m2)
-
-  logical :: use_temperature   ! If true, temp and saln used as state variables.
+    u => NULL(), & !< u : zonal velocity component (m/s)
+    v => NULL(), & !< v : meridional velocity component (m/s)
+    h => NULL()    !< h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss))
+  real :: depth(SZI_(CS%G))           !< distance from the surface (meter)
+  real :: depth_ml                    !< depth over which to average to
+                                      !< determine mixed layer properties (meter)
+  real :: dh                          !< thickness of a layer within mixed layer (meter)
+  real :: mass                        !< mass per unit area of a layer (kg/m2)
+  real :: T_freeze                    !< freezing temperature (oC)
+  logical :: use_temperature   !< If true, temp and saln used as state variables.
   integer :: i, j, k, is, ie, js, je, nz, numberOfErrors
   integer :: isd, ied, jsd, jed
   integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB
@@ -2831,6 +2831,22 @@ subroutine extract_surface_state(CS, sfc_state)
     endif
   endif  ! (CS%Hmix >= 0.0)
 
+  if (allocated(sfc_state%melt_potential)) then
+    !$OMP parallel do default(shared)
+    do j=js,je ; do i=is,ie
+      ! set melt_potential to zero to avoid passing values set previously
+      if (G%mask2dT(i,j)>0.) then
+        ! calculate freezing pot. temp. @ surface
+        call calculate_TFreeze(sfc_state%SSS(i,j), 0.0, T_freeze, CS%tv%eqn_of_state)
+        ! time accumulated melt_potential, in J/m^2
+        sfc_state%melt_potential(i,j) = sfc_state%melt_potential(i,j) +  (CS%tv%C_p * CS%GV%Rho0 * &
+                                        (sfc_state%SST(i,j) - T_freeze) * CS%Hmix)
+      else
+          sfc_state%melt_potential(i,j) = 0.0
+      endif! G%mask2dT
+    enddo ; enddo
+  endif
+
   if (allocated(sfc_state%salt_deficit) .and. associated(CS%tv%salt_deficit)) then
     !$OMP parallel do default(shared)
     do j=js,je ; do i=is,ie

src/core/MOM_variables.F90

@@ -35,10 +35,13 @@ module MOM_variables
     SST, &         !< The sea surface temperature in C.
     SSS, &         !< The sea surface salinity in psu.
     sfc_density, & !< The mixed layer density in kg m-3.
-    Hml, &         !< The mixed layer depth in m.
-    u, &           !< The mixed layer zonal velocity in m s-1.
-    v, &           !< The mixed layer meridional velocity in m s-1.
-    sea_lev, &     !< The sea level in m. If a reduced surface gravity is used, it is compensated in sea_lev.
+    Hml, &      !< The mixed layer depth in m.
+    u, &        !< The mixed layer zonal velocity in m s-1.
+    v, &        !< The mixed layer meridional velocity in m s-1.
+    sea_lev, &  !< The sea level in m.  If a reduced surface gravity is
+                !! used, that is compensated for in sea_lev.
+    melt_potential, & !< Amount of heat that can be used to melt sea ice, in J m-2.
+                      !! This is computed w.r.t. surface freezing temperature.
     ocean_mass, &  !< The total mass of the ocean in kg m-2.
     ocean_heat, &  !< The total heat content of the ocean in C kg m-2.
     ocean_salt, &  !< The total salt content of the ocean in kgSalt m-2.
@@ -289,20 +292,22 @@ module MOM_variables
 !> Allocates the fields for the surface (return) properties of
 !! the ocean model. Unused fields are unallocated.
 subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
-                                  gas_fields_ocn)
-  type(ocean_grid_type), intent(in)    :: G               !< ocean grid structure
-  type(surface),         intent(inout) :: sfc_state       !< ocean surface state type to be allocated.
-  logical,     optional, intent(in)    :: use_temperature !< If true, allocate the space for thermodynamic variables.
-  logical,     optional, intent(in)    :: do_integrals    !< If true, allocate the space for vertically
-                                                          !! integrated fields.
+                                  gas_fields_ocn, use_meltpot)
+  type(ocean_grid_type), intent(in)    :: G                !< ocean grid structure
+  type(surface),         intent(inout) :: sfc_state        !< ocean surface state type to be allocated.
+  logical,     optional, intent(in)    :: use_temperature  !< If true, allocate the space for thermodynamic variables.
+  logical,     optional, intent(in)    :: do_integrals     !< If true, allocate the space for vertically
+                                                           !! integrated fields.
   type(coupler_1d_bc_type), &
                optional, intent(in)    :: gas_fields_ocn  !< If present, this type describes the ocean
                                               !! ocean and surface-ice fields that will participate
                                               !! in the calculation of additional gas or other
                                               !! tracer fluxes, and can be used to spawn related
                                               !! internal variables in the ice model.
+  logical,     optional, intent(in)    :: use_meltpot      !< If true, allocate the space for melt potential
 
-  logical :: use_temp, alloc_integ
+  ! local variables
+  logical :: use_temp, alloc_integ, use_melt_potential
   integer :: is, ie, js, je, isd, ied, jsd, jed
   integer :: isdB, iedB, jsdB, jedB
 
@@ -312,6 +317,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
 
   use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature
   alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals
+  use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot
 
   if (sfc_state%arrays_allocated) return
 
@@ -326,6 +332,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, &
   allocate(sfc_state%u(IsdB:IedB,jsd:jed)) ; sfc_state%u(:,:) = 0.0
   allocate(sfc_state%v(isd:ied,JsdB:JedB)) ; sfc_state%v(:,:) = 0.0
 
+  if (use_melt_potential) then
+    allocate(sfc_state%melt_potential(isd:ied,jsd:jed)) ; sfc_state%melt_potential(:,:) = 0.0
+  endif
+
   if (alloc_integ) then
     ! Allocate structures for the vertically integrated ocean_mass, ocean_heat, and ocean_salt.
     allocate(sfc_state%ocean_mass(isd:ied,jsd:jed)) ; sfc_state%ocean_mass(:,:) = 0.0
@@ -350,6 +360,7 @@ subroutine deallocate_surface_state(sfc_state)
 
   if (.not.sfc_state%arrays_allocated) return
 
+  if (allocated(sfc_state%melt_potential)) deallocate(sfc_state%melt_potential)
   if (allocated(sfc_state%SST)) deallocate(sfc_state%SST)
   if (allocated(sfc_state%SSS)) deallocate(sfc_state%SSS)
   if (allocated(sfc_state%sfc_density)) deallocate(sfc_state%sfc_density)