From a840dc0f45810d363713568b90f17606839b970f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Jul 2018 07:12:23 -0600 Subject: [PATCH 01/12] Pass melt potential (o2x_Fioo_q) via mct driver --- config_src/mct_driver/MOM_ocean_model.F90 | 30 ++++++++++++++--------- config_src/mct_driver/ocn_cap_methods.F90 | 18 +++++++++++++- config_src/mct_driver/ocn_comp_mct.F90 | 11 +++++++-- config_src/mct_driver/ocn_cpl_indices.F90 | 2 +- 4 files changed, 45 insertions(+), 16 deletions(-) diff --git a/config_src/mct_driver/MOM_ocean_model.F90 b/config_src/mct_driver/MOM_ocean_model.F90 index c894f42270..f7092865c6 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 (W/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..fdc1a619b4 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,14 @@ 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 already is in W/m^2 (ncouple_per_day is unitless) + o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * ncouple_per_day + 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..f924f50078 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -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 From 3cdd97bfab98e62ccb6b1dbd2eb21334975d9d2d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Jul 2018 07:22:49 -0600 Subject: [PATCH 02/12] Add option to calculate melt potential --- src/core/MOM.F90 | 47 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index bdd1f159cf..f3abb6bc69 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -72,7 +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_EOS, only : EOS_init, calculate_density +use MOM_EOS, only : EOS_init, calculate_density, calculate_TFreeze use MOM_debugging, only : check_redundant use MOM_grid, only : ocean_grid_type, set_first_direction use MOM_grid, only : MOM_grid_init, MOM_grid_end @@ -807,7 +807,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & endif if (showCallTree) call callTree_waypoint("calling extract_surface_state (step_MOM)") - call extract_surface_state(CS, sfc_state) + call extract_surface_state(CS, sfc_state, dt) ! Do diagnostics that only occur at the end of a complete forcing step. if (cycle_end) then @@ -2631,28 +2631,29 @@ end subroutine adjust_ssh_for_p_atm !> This subroutine sets the surface (return) properties of the ocean !! model by setting the appropriate fields in sfc_state. Unused fields !! are set to NULL or are unallocated. -subroutine extract_surface_state(CS, sfc_state) +subroutine extract_surface_state(CS, sfc_state, dt) type(MOM_control_struct), pointer :: CS !< Master MOM control structure type(surface), intent(inout) :: sfc_state !< transparent ocean surface state !! structure shared with the calling routine !! data in this structure is intent out. + real, optional, intent(in) :: dt !< Thermodynamic time step, in s. ! 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 @@ -2810,6 +2811,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 + sfc_state%melt_potential(i,j) = 0.0 + ! calculate freezing temp. + call calculate_TFreeze(sfc_state%SSS(i,j), CS%tv%P_Ref, T_freeze, CS%tv%eqn_of_state) + if (present(dt)) then + ! melt_potential, in W/m^2 + sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%GV%Rho0 * (sfc_state%SST(i,j) - T_freeze) * sfc_state%Hml(i,j)/dt + else + sfc_state%melt_potential(i,j) = 0.0 + endif + 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 From b52db1ff5dfc0e06b879e53ba7062c06c9649950 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Jul 2018 07:40:16 -0600 Subject: [PATCH 03/12] Allocate/deallocate melt_potential --- src/core/MOM_variables.F90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 02b0b622a3..32a1f75c0c 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -38,6 +38,8 @@ module MOM_variables 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, & !< Instantaneous amount of heat that can be used to melt sea ice, + !! in J m-2. This is computed w.r.t. the 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. @@ -280,7 +282,7 @@ module MOM_variables !> This subroutine 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) + 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. @@ -292,8 +294,10 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! 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 @@ -303,6 +307,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 @@ -317,6 +322,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. @@ -342,6 +351,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) From d981745f8d87b7d3ad7381df7409a38203c631b3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Jul 2018 09:11:37 -0600 Subject: [PATCH 04/12] Fix a call to ocn_export --- config_src/mct_driver/ocn_comp_mct.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index f924f50078..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') From 60eb7ec30667a000f23c9ab5af751b058cce1906 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 30 Jul 2018 09:21:43 -0600 Subject: [PATCH 05/12] Add comment about icefrq used in Hycom --- config_src/mct_driver/ocn_cap_methods.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index fdc1a619b4..36a4faef09 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/config_src/mct_driver/ocn_cap_methods.F90 @@ -182,6 +182,7 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day) o2x(ind%o2x_Fioo_q, n) = ocn_public%frazil(ig,jg) * grid%mask2dT(i,j) * I_time_int else ! Melt_potential already is in W/m^2 (ncouple_per_day is unitless) + ! GMM - Hycom cap uses icefrq rather than oceanfrq (ncouple_per_day) o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * ncouple_per_day endif ! Make a copy of ssh in order to do a halo update. We use the usual MOM domain From 4e5bf74943a363dc84710b0ffee0cb2d39187209 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 1 Aug 2018 13:09:49 -0600 Subject: [PATCH 06/12] Change calculation of melt_potential sfc_state%melt_potential is now accumulated over time_in_thermo_cycle. --- config_src/mct_driver/MOM_ocean_model.F90 | 2 +- config_src/mct_driver/ocn_cap_methods.F90 | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/config_src/mct_driver/MOM_ocean_model.F90 b/config_src/mct_driver/MOM_ocean_model.F90 index f7092865c6..9971747afd 100644 --- a/config_src/mct_driver/MOM_ocean_model.F90 +++ b/config_src/mct_driver/MOM_ocean_model.F90 @@ -722,7 +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 (W/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 diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index 36a4faef09..5d2aac1317 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/config_src/mct_driver/ocn_cap_methods.F90 @@ -181,9 +181,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day) ! 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 already is in W/m^2 (ncouple_per_day is unitless) - ! GMM - Hycom cap uses icefrq rather than oceanfrq (ncouple_per_day) - o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * ncouple_per_day + ! 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 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. From 7347ad973cc518173e02f04e28bb27c17d56887c Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 1 Aug 2018 13:12:22 -0600 Subject: [PATCH 07/12] Update calculation of melt_potential --- src/core/MOM.F90 | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f3abb6bc69..73b09820b6 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -195,7 +195,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 @@ -526,6 +525,8 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (therm_reset) then CS%time_in_thermo_cycle = 0.0 + ! GMM + 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 @@ -807,7 +808,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & endif if (showCallTree) call callTree_waypoint("calling extract_surface_state (step_MOM)") - call extract_surface_state(CS, sfc_state, dt) + call extract_surface_state(CS, sfc_state, dt_therm) ! Do diagnostics that only occur at the end of a complete forcing step. if (cycle_end) then @@ -2815,15 +2816,17 @@ subroutine extract_surface_state(CS, sfc_state, dt) !$OMP parallel do default(shared) do j=js,je ; do i=is,ie ! set melt_potential to zero to avoid passing values set previously - sfc_state%melt_potential(i,j) = 0.0 - ! calculate freezing temp. - call calculate_TFreeze(sfc_state%SSS(i,j), CS%tv%P_Ref, T_freeze, CS%tv%eqn_of_state) + 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) if (present(dt)) then - ! melt_potential, in W/m^2 - sfc_state%melt_potential(i,j) = CS%tv%C_p * CS%GV%Rho0 * (sfc_state%SST(i,j) - T_freeze) * sfc_state%Hml(i,j)/dt + ! 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 + endif! G%mask2dT enddo ; enddo endif From 05d6c805abc50fbb6dc61d62ed128400c1f0c108 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 1 Aug 2018 14:26:33 -0600 Subject: [PATCH 08/12] Update melt_potential again, but values are still ~ 3 x smaller --- config_src/mct_driver/ocn_cap_methods.F90 | 2 +- src/core/MOM.F90 | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index 5d2aac1317..5b7c341424 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/config_src/mct_driver/ocn_cap_methods.F90 @@ -182,7 +182,7 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day) 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 + o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) !* I_time_int * ncouple_per_day 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. diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 73b09820b6..2c7201028d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -525,7 +525,6 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (therm_reset) then CS%time_in_thermo_cycle = 0.0 - ! GMM 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 @@ -2820,9 +2819,9 @@ subroutine extract_surface_state(CS, sfc_state, dt) ! calculate freezing pot. temp. @ surface call calculate_TFreeze(sfc_state%SSS(i,j), 0.0, T_freeze, CS%tv%eqn_of_state) if (present(dt)) then - ! time accumulated melt_potential, in J/m^2 + ! time accumulated melt_potential, in W/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) + (sfc_state%SST(i,j) - T_freeze) * CS%Hmix)/dt else sfc_state%melt_potential(i,j) = 0.0 endif From be0e7a844aada4541cc9f8dec341b5406bc4bccd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 3 Aug 2018 16:48:24 -0600 Subject: [PATCH 09/12] Remove dt from extract_surface_state --- src/core/MOM.F90 | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2c7201028d..a6206c2fec 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -807,7 +807,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & endif if (showCallTree) call callTree_waypoint("calling extract_surface_state (step_MOM)") - call extract_surface_state(CS, sfc_state, dt_therm) + call extract_surface_state(CS, sfc_state) ! Do diagnostics that only occur at the end of a complete forcing step. if (cycle_end) then @@ -2631,12 +2631,11 @@ end subroutine adjust_ssh_for_p_atm !> This subroutine sets the surface (return) properties of the ocean !! model by setting the appropriate fields in sfc_state. Unused fields !! are set to NULL or are unallocated. -subroutine extract_surface_state(CS, sfc_state, dt) +subroutine extract_surface_state(CS, sfc_state) type(MOM_control_struct), pointer :: CS !< Master MOM control structure type(surface), intent(inout) :: sfc_state !< transparent ocean surface state !! structure shared with the calling routine !! data in this structure is intent out. - real, optional, intent(in) :: dt !< Thermodynamic time step, in s. ! local real :: hu, hv @@ -2818,13 +2817,11 @@ subroutine extract_surface_state(CS, sfc_state, dt) 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) - if (present(dt)) then - ! time accumulated melt_potential, in W/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)/dt - else + ! time accumulated melt_potential, in W/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 endif! G%mask2dT enddo ; enddo endif From eb2c69f1dde9c955257f98fc183d62f606cf3e64 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 3 Aug 2018 16:53:25 -0600 Subject: [PATCH 10/12] Replace units (W/m^2 to J/m^2) in the comments --- src/core/MOM.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a6206c2fec..5a648c4c9c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2817,7 +2817,7 @@ subroutine extract_surface_state(CS, sfc_state) 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 W/m^2 + ! 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 From c6b253400839a91678d4278519ad9c62131e47f9 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 3 Aug 2018 17:02:21 -0600 Subject: [PATCH 11/12] Pass melt/freeze potential to coupler; reset melt_potential --- config_src/mct_driver/ocn_cap_methods.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index 5b7c341424..6e3317f376 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/config_src/mct_driver/ocn_cap_methods.F90 @@ -182,7 +182,9 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day) 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 + o2x(ind%o2x_Fioo_q, n) = -ocn_public%melt_potential(ig,jg) * grid%mask2dT(i,j) * I_time_int !* ncouple_per_day + ! reset melt_potential + ocn_public%melt_potential(ig,jg) = 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. From d0bc3bc2fcbfec4a6646bd3f21475e5c87d5f230 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 6 Aug 2018 14:25:49 -0600 Subject: [PATCH 12/12] Fix bug and make sure melt potential is always <= 0 --- config_src/mct_driver/ocn_cap_methods.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config_src/mct_driver/ocn_cap_methods.F90 b/config_src/mct_driver/ocn_cap_methods.F90 index 6e3317f376..17894dc966 100644 --- a/config_src/mct_driver/ocn_cap_methods.F90 +++ b/config_src/mct_driver/ocn_cap_methods.F90 @@ -183,8 +183,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day) 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 - ! reset melt_potential - ocn_public%melt_potential(ig,jg) = 0.0 + ! 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.