diff --git a/config_src/mct_driver/MOM_ocean_model.F90 b/config_src/mct_driver/MOM_ocean_model.F90 index c894f42270..9971747afd 100644 --- a/config_src/mct_driver/MOM_ocean_model.F90 +++ b/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,13 +707,14 @@ 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 @@ -720,6 +722,7 @@ subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, 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) diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index cc214306f0..17894dc966 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/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) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 63a24b153d..7608ef4579 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/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 diff --git a/config_src/mct_driver/ocn_cpl_indices.F90 b/config_src/mct_driver/ocn_cpl_indices.F90 index 4bd9c1f383..52f94f6106 100644 --- a/config_src/mct_driver/ocn_cpl_indices.F90 +++ b/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 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c554e4f92e..dbe20bd861 100644 --- a/src/core/MOM.F90 +++ b/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 diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 4165fb0e11..1b9dae9a6d 100644 --- a/src/core/MOM_variables.F90 +++ b/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)