From 14e233438ff175813e4593474e2de1206f62b648 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Sat, 22 May 2021 08:53:35 -0600 Subject: [PATCH 01/28] Deallocate Kd_bkgnd and Kv_bkgnd These arrays were not being deallocated which was causing a memory leak problem. --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 0cb39b2f15..d3ced12160 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -694,6 +694,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%drho_rat)) deallocate(dd%drho_rat) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) + if (associated(dd%Kd_bkgnd)) deallocate(dd%Kd_bkgnd) + if (associated(dd%Kv_bkgnd)) deallocate(dd%Kv_bkgnd) if (showCallTree) call callTree_leave("set_diffusivity()") From 0d591560846c96bb3b667e2f59bc7324bfe106ed Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Sun, 13 Jun 2021 22:56:47 -0600 Subject: [PATCH 02/28] introduce a new wave coupling method: VR12-MA (not fully implemented yet) --- config_src/drivers/nuopc_cap/mom_cap.F90 | 49 ++++++++++--------- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 16 ++++-- src/user/MOM_wave_interface.F90 | 14 +++++- 3 files changed, 51 insertions(+), 28 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..8d383c3d3f 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -696,17 +696,21 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%frunoff = 0.0 if (ocean_state%use_waves) then - Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands - allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & - Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & - Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & - Ice_ocean_boundary% vstkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & - Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) - Ice_ocean_boundary%ustk0 = 0.0 - Ice_ocean_boundary%vstk0 = 0.0 - Ice_ocean_boundary%stk_wavenumbers = ocean_state%Waves%WaveNum_Cen - Ice_ocean_boundary%ustkb = 0.0 - Ice_ocean_boundary%vstkb = 0.0 + if (ocean_state%wave_method == "VR12-MA") then + continue + else + Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands + allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & + Ice_ocean_boundary% vstk0 (isc:iec,jsc:jec), & + Ice_ocean_boundary% ustkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & + Ice_ocean_boundary% vstkb (isc:iec,jsc:jec,Ice_ocean_boundary%num_stk_bands), & + Ice_ocean_boundary%stk_wavenumbers (Ice_ocean_boundary%num_stk_bands)) + Ice_ocean_boundary%ustk0 = 0.0 + Ice_ocean_boundary%vstk0 = 0.0 + Ice_ocean_boundary%stk_wavenumbers = ocean_state%Waves%WaveNum_Cen + Ice_ocean_boundary%ustkb = 0.0 + Ice_ocean_boundary%vstkb = 0.0 + endif endif ocean_internalstate%ptr%ocean_state_type_ptr => ocean_state @@ -723,9 +727,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_hstokes" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_melth" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_meltw" , "will provide") - !call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_fswpen" , "will provide") else !call fld_list_add(fldsToOcn_num, fldsToOcn, "mass_of_overlying_sea_ice" , "will provide") !call fld_list_add(fldsFrOcn_num, fldsFrOcn, "sea_lev" , "will provide") @@ -753,15 +754,19 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") if (ocean_state%use_waves) then - if (Ice_ocean_boundary%num_stk_bands > 3) then - call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this") + if (ocean_state%wave_method == "VR12-MA") then + continue + else + if (Ice_ocean_boundary%num_stk_bands > 3) then + call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this") + endif + call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_1" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_1", "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_2" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_2", "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_3" , "will provide") + call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_3", "will provide") endif - call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_1" , "will provide") - call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_1", "will provide") - call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_2" , "will provide") - call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_2", "will provide") - call fld_list_add(fldsToOcn_num, fldsToOcn, "eastward_partitioned_stokes_drift_3" , "will provide") - call fld_list_add(fldsToOcn_num, fldsToOcn, "northward_partitioned_stokes_drift_3", "will provide") endif !--------- export fields ------------- diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 493762f4bc..7e7716d386 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -145,6 +145,7 @@ module MOM_ocean_model_nuopc integer :: nstep = 0 !< The number of calls to update_ocean. logical :: use_ice_shelf !< If true, the ice shelf model is enabled. logical,public :: use_waves !< If true use wave coupling. + character(len=40),public :: wave_method !< Wave coupling method. logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the !! ocean dynamics and forcing fluxes. @@ -387,10 +388,15 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & "If true, enables surface wave modules.", default=.false.) if (OS%use_waves) then + call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.) call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) - call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS",OS%Waves%WaveNum_Cen, & - "Central wavenumbers for surface Stokes drift bands.",units='rad/m', & - default=0.12566) + if (OS%wave_method == "VR12-MA") then + continue + else + call get_param(param_file,mdl,"SURFBAND_WAVENUMBERS",OS%Waves%WaveNum_Cen, & + "Central wavenumbers for surface Stokes drift bands.",units='rad/m', & + default=0.12566) + endif else call MOM_wave_interface_init_lite(param_file) endif @@ -575,7 +581,9 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid, OS%US) if (OS%use_waves) then - call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) + if (OS%wave_method /= "VR12-MA") then + call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) + endif endif if (OS%nstep==0) then diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index a8e3e207a6..1ec1e9314c 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -141,6 +141,7 @@ module MOM_wave_interface !! 1 - Surface Stokes Drift Bands !! 2 - DHH85 !! 3 - LF17 + !! 4 - VR12MA !! -99 - No waves computed, but empirical Langmuir number used. !! \todo Module variable! Move into a control structure. @@ -178,7 +179,7 @@ module MOM_wave_interface !! into a control structure. ! Switches needed in import_stokes_drift integer, parameter :: TESTPROF = 0, SURFBANDS = 1, & - DHH85 = 2, LF17 = 3, NULL_WaveMethod=-99, & + DHH85 = 2, LF17 = 3, VR12MA = 4, NULL_WaveMethod=-99, & DATAOVR = 1, COUPLER = 2, INPUT = 3 ! Options For Test Prof @@ -208,6 +209,7 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" character*(5), parameter :: DHH85_STRING = "DHH85" character*(4), parameter :: LF17_STRING = "LF17" + character*(7), parameter :: VR12MA_STRING = "VR12-MA" character*(12), parameter :: DATAOVR_STRING = "DATAOVERRIDE" character*(7), parameter :: COUPLER_STRING = "COUPLER" character*(5), parameter :: INPUT_STRING = "INPUT" @@ -268,7 +270,11 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) " DHH85 - Uses Donelan et al. 1985 empirical \n"// & " wave spectrum with prescribed values. \n"// & " LF17 - Infers Stokes drift profile from wind \n"// & - " speed following Li and Fox-Kemper 2017.\n", & + " speed following Li and Fox-Kemper 2017.\n"// & + " VR12-MA - Applies an enhancement factor to the KPP\n"//& + " turbulent velocity scale received from \n"// & + " WW3 and is based on the surface layer \n"// & + " and projected Langmuir number. (Li 2016)\n", & units='', default=NULL_STRING) select case (TRIM(TMPSTRING1)) case (NULL_STRING)! No Waves @@ -363,6 +369,8 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) default=.false.) case (LF17_STRING)!Li and Fox-Kemper 17 wind-sea Langmuir number WaveMethod = LF17 + case (VR12MA_STRING)!Li and Fox-Kemper 16 + WaveMethod = VR12MA case default call MOM_error(FATAL,'Check WAVE_METHOD.') end select @@ -732,6 +740,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo DHH85_is_set = .true. endif + elseif (WaveMethod==VR12MA) then + return ! todo else! Keep this else, fallback to 0 Stokes drift do kk= 1,GV%ke do jj = G%jsd,G%jed From 276954f3cf5005858a713e324c348b7fd677aac4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 14 Jun 2021 08:40:19 -0600 Subject: [PATCH 03/28] Fixed a misplaced parentheses This commit fixes misplaced parentheses in the tridiagonal solvers used in subroutines tracer_vertdiff_Eulerian and tracer_vertdiff. This bug has been reported in https://github.com/NOAA-GFDL/MOM6/issues/1415 This commit also changes the mask condition used in these subroutines. --- src/tracer/MOM_tracer_diabatic.F90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index 9be4af08dc..cbfe2d7f0a 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -152,7 +152,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) @@ -185,14 +185,14 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ea(i,j,1) b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ea(i,j,k) @@ -200,7 +200,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = eb(i,j,nz-1) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ea(i,j,nz) @@ -208,7 +208,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ea(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo @@ -267,8 +267,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & !! ensure positive definiteness [H ~> m or kg m-2]. real :: h_neglect !< A thickness that is so small it is usually lost !! in roundoff and can be neglected [H ~> m or kg m-2]. - logical :: convert_flux = .true. - + logical :: convert_flux integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -279,6 +278,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & return endif + convert_flux = .true. if (present(convert_flux_in)) convert_flux = convert_flux_in h_neglect = GV%H_subroundoff sink_dist = 0.0 @@ -351,7 +351,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) @@ -384,14 +384,14 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ent(i,j,1) b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ent(i,j,K) @@ -399,7 +399,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = ent(i,j,nz) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ent(i,j,nz) @@ -407,7 +407,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ent(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo From 3aade32dd6b8033603d3e1b21678dfc26c8d3281 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 14 Jun 2021 12:30:37 -0600 Subject: [PATCH 04/28] CFCs implementation via NUOPC cap This commit adds a new module (MOM_CFC_cap) that can used to simulate chlorofluorocarbons (CFCs) tracers. It also includes the modifications needed to 1) allocate arrays in the fluxes and surface types, and 2) activate the MOM_CFC_cap module via NUOPC cap. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 12 +- .../drivers/nuopc_cap/mom_cap_methods.F90 | 17 +- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 14 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 66 +- src/core/MOM_forcing_type.F90 | 74 +- src/core/MOM_unit_tests.F90 | 4 +- src/core/MOM_variables.F90 | 17 +- src/tracer/MOM_CFC_cap.F90 | 704 ++++++++++++++++++ src/tracer/MOM_tracer_flow_control.F90 | 33 +- 9 files changed, 922 insertions(+), 19 deletions(-) create mode 100644 src/tracer/MOM_CFC_cap.F90 diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..42cc579546 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -31,7 +31,7 @@ module MOM_cap_mod use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type -use MOM_ocean_model_nuopc, only: ocean_model_init_sfc +use MOM_ocean_model_nuopc, only: ocean_model_init_sfc, ocean_model_flux_init use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh use MOM_cap_time, only: AlarmInit @@ -649,6 +649,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ocean_public%is_ocean_pe = .true. call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles)) + ! GMM, this call is not needed for NCAR. Check with EMC. + ! If this can be deleted, perhaps we should also delete ocean_model_flux_init + call ocean_model_flux_init(ocean_state) + call ocean_model_init_sfc(ocean_state, ocean_public) call mpp_get_compute_domain(ocean_public%domain, isc, iec, jsc, jec) @@ -668,6 +672,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),& Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), & Ice_ocean_boundary% mi (isc:iec,jsc:jec), & + Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec), & + Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec), & Ice_ocean_boundary% p (isc:iec,jsc:jec), & Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec), & Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec), & @@ -689,6 +695,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%seaice_melt = 0.0 Ice_ocean_boundary%seaice_melt_heat= 0.0 Ice_ocean_boundary%mi = 0.0 + Ice_ocean_boundary%ice_fraction = 0.0 + Ice_ocean_boundary%u10_sqr = 0.0 Ice_ocean_boundary%p = 0.0 Ice_ocean_boundary%lrunoff_hflx = 0.0 Ice_ocean_boundary%frunoff_hflx = 0.0 @@ -747,6 +755,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface" , "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl" , "will provide") !-> liquid runoff call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi" , "will provide") !-> ice runoff + call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac" , "will provide") !-> ice fraction + call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n" , "will provide") !-> wind^2 at 10m call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn" , "will provide") !These are not currently used and changing requires a nuopc dictionary change diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 0a4e12c647..b7e3b13c1e 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -250,12 +250,27 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ! Note - preset values to 0, if field does not exist in importState, then will simply return ! and preset value will be used - ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mass_of_overlying_ice', & isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- + ! sea-ice fraction + !---- + ice_ocean_boundary%ice_fraction(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Si_ifrac', & + isc, iec, jsc, jec, ice_ocean_boundary%ice_fraction, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---- + ! 10m wind squared + !---- + ice_ocean_boundary%u10_sqr(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'So_duu10n', & + isc, iec, jsc, jec, ice_ocean_boundary%u10_sqr, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- ! Partitioned Stokes Drift Components !---- diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 493762f4bc..9d5f81f079 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -40,7 +40,6 @@ module MOM_ocean_model_nuopc use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real use time_interp_external_mod,only : time_interp_external_init -use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface @@ -240,6 +239,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! 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 + ! Local variables real :: Rho0 ! The Boussinesq ocean density, in kg m-3. real :: G_Earth ! The gravitational acceleration in m s-2. @@ -247,7 +247,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! The actual depth over which melt potential is computed will !! min(HFrz, OBLD), where OBLD is the boundary layer depth. !! If HFrz <= 0 (default), melt potential will not be computed. - logical :: use_melt_pot!< If true, allocate melt_potential array + logical :: use_melt_pot !< If true, allocate melt_potential array + logical :: use_CFC !< If true, allocated arrays for surface CFCs. + ! This include declares and sets the variable "version". #include "version_variable.h" @@ -366,10 +368,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i use_melt_pot=.false. endif + call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, & + default=.false., do_not_log=.true.) + ! 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, use_meltpot=use_melt_pot) + do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & + use_meltpot=use_melt_pot, use_cfcs=use_CFC) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) @@ -416,6 +422,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif + call extract_surface_state(OS%MOM_CSp, OS%sfc_state) + call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 7168823fbc..1054234e6b 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -21,6 +21,7 @@ module MOM_surface_forcing_nuopc use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state @@ -79,7 +80,7 @@ module MOM_surface_forcing_nuopc !! the correction for the atmospheric (and sea-ice) !! pressure limited by max_p_surf instead of the !! full atmospheric pressure. The default is true. - + logical :: use_CFC !< enables the MOM_CFC_cap tracer package. real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa] logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied !! from an input file. @@ -128,6 +129,7 @@ module MOM_surface_forcing_nuopc type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing character(len=200) :: inputdir !< directory where NetCDF input files are + character(len=200) :: CFC_BC_file !< filename with cfc11 and cfc12 data character(len=200) :: salt_restore_file !< filename for salt restoring data character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file logical :: mask_srestore !< if true, apply a 2-dimensional mask to the surface @@ -141,9 +143,13 @@ module MOM_surface_forcing_nuopc !! temperature restoring fluxes. The masking file should be !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' + character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file + character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_cfc11_atm = -1 !< id number for time_interp_external. + integer :: id_cfc12_atm = -1 !< id number for time_interp_external. ! Diagnostics handles type(forcing_diags), public :: handles @@ -178,6 +184,8 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: frunoff_hflx =>NULL() !< heat content of frozen runoff [W/m2] real, pointer, dimension(:,:) :: p =>NULL() !< pressure of overlying ice and atmosphere !< on ocean surface [Pa] + real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< mass of ice [nondim] + real, pointer, dimension(:,:) :: u10_sqr =>NULL() !< wind speed squared at 10m [m2/s2] real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] real, pointer, dimension(:,:) :: ice_rigidity =>NULL() !< rigidity of the sea ice, sea-ice and !! ice-shelves, expressed as a coefficient @@ -234,6 +242,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! local variables real, dimension(SZI_(G),SZJ_(G)) :: & + cfc11_atm, & !< CFC11 concentration in the atmopshere [???????] + cfc12_atm, & !< CFC11 concentration in the atmopshere [???????] data_restore, & !< The surface value toward which to restore [g/kg or degC] SST_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C] SSS_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg] @@ -293,7 +303,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! flux type has been used. if (fluxes%dt_buoy_accum < 0) then call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & - press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug) + press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & + cfc=CS%use_CFC) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -433,6 +444,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo endif + ! Check that liquid runoff has a place to go if (CS%liquid_runoff_from_data .and. .not. associated(IOB%lrunoff)) then call MOM_error(FATAL, "liquid runoff is being added via data_override but "// & @@ -496,6 +508,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%seaice_melt)) & fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0) + ! sea ice fraction [nondim] + if (associated(IOB%ice_fraction)) & + fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) + + ! 10-m wind speed squared [m2/s2] + if (associated(IOB%u10_sqr)) & + fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) + fluxes%latent(i,j) = 0.0 ! notice minus sign since fprec is positive into the ocean if (associated(IOB%fprec)) then @@ -550,6 +570,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure. endif + ! CFCs + if (CS%use_CFC) then + call CFC_cap_fluxes(fluxes, sfc_state, G, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + endif + if (associated(IOB%salt_flux)) then do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*IOB%salt_flux(i-i0,j-j0)) @@ -1161,6 +1186,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "coupler. This is used for testing and should be =1.0 for any "//& "production runs.", default=1.0) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, & + default=.false., do_not_log=.true.) + if (restore_salt) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//& @@ -1173,7 +1201,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "The name of the surface salinity variable to read from "//& "SALT_RESTORE_FILE for restoring salinity.", & default="salt") - call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & "If true, the restoring of salinity is applied as a salt "//& "flux instead of as a freshwater flux.", default=.false.) @@ -1269,7 +1296,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, enddo ; enddo endif - call time_interp_external_init + ! initialize time interpolator module + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global ! constant. @@ -1320,8 +1348,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, & "If true, makes available diagnostics of fluxes from icebergs "//& "as seen by MOM6.", default=.false.) + call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, & - use_berg_fluxes=iceberg_flux_diags) + use_berg_fluxes=iceberg_flux_diags, use_cfcs=CS%use_CFC) call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & "If true, allows flux adjustments to specified via the "//& @@ -1355,6 +1384,27 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, endif endif ; endif + ! Do not log these params here since they are logged in the CFC cap module + if (CS%use_CFC) then + call get_param(param_file, mdl, "CFC_BC_FILE", CS%CFC_BC_file, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ", do_not_log=.true.) + if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then + ! Add the directory if CFC_BC_file is not already a complete path. + CS%CFC_BC_file = trim(slasher(CS%inputdir))//trim(CS%CFC_BC_file) + call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.) + call get_param(param_file, mdl, "CFC12_VARIABLE", CS%cfc12_var_name, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) + + CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) + CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) + endif + endif + ! Set up any restart fields associated with the forcing. call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") call restart_init_end(CS%restart_CSp) @@ -1424,6 +1474,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) write(outunit,100) 'iobt%lrunoff ' , mpp_chksum( iobt%lrunoff ) write(outunit,100) 'iobt%frunoff ' , mpp_chksum( iobt%frunoff ) write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p ) + write(outunit,100) 'iobt%ice_fraction ' , mpp_chksum( iobt%ice_fraction ) + write(outunit,100) 'iobt%u10_sqr ' , mpp_chksum( iobt%u10_sqr ) if (associated(iobt%ustar_berg)) & write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 682ad03397..7d8375a8af 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -183,6 +183,13 @@ module MOM_forcing_type real :: C_p !< heat capacity of seawater [Q degC-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. + ! CFC-related arrays needed in the MOM_CFC_cap module + real, pointer, dimension(:,:) :: & + cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [mol m-2 s-1]. + cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [mol m-2 s-1]. + ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. + u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] + ! passive tracer surface fluxes type(coupler_2d_bc_type) :: tr_fluxes !< This structure contains arrays of !! of named fields used for passive tracer fluxes. @@ -348,6 +355,12 @@ module MOM_forcing_type integer :: id_TKE_tidal = -1 integer :: id_buoy = -1 + ! cfc-related diagnostics handles + integer :: id_cfc11 = -1 + integer :: id_cfc12 = -1 + integer :: id_ice_fraction = -1 + integer :: id_u10_sqr = -1 + ! iceberg diagnostic handles integer :: id_ustar_berg = -1 integer :: id_area_berg = -1 @@ -1087,6 +1100,10 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%p_surf)) & call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift , scale=US%RL2_T2_to_Pa) + if (associated(fluxes%u10_sqr)) & + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%Z_to_m**2*US%s_to_T**2) + if (associated(fluxes%ice_fraction)) & + call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & @@ -1239,13 +1256,14 @@ end subroutine forcing_SinglePointPrint !> Register members of the forcing type for diagnostics -subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_cfcs) type(time_type), intent(in) :: Time !< time type type(diag_ctrl), intent(inout) :: diag !< diagnostic control type type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< True if T/S are in use type(forcing_diags), intent(inout) :: handles !< handles for diagnostics logical, optional, intent(in) :: use_berg_fluxes !< If true, allow iceberg flux diagnostics + logical, optional, intent(in) :: use_cfcs !< If true, allow cfc related diagnostics ! Clock for forcing diagnostics handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1288,6 +1306,34 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, endif endif + ! units for cfc11_flux and cfc12_flux are mol m-2 s-1 + ! See: + ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html + ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html + if (present(use_cfcs)) then + if (use_cfcs) then + handles%id_cfc11 = register_diag_field('ocean_model', 'cfc11_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC11 into the ocean ', 'mol m-2 s-1', & + conversion= US%m_to_Z**2*US%s_to_T,& + cmor_field_name='fgcfc11', & + cmor_long_name='Surface Downward CFC11 Flux', & + cmor_standard_name='surface_downward_cfc11_flux') + + handles%id_cfc12 = register_diag_field('ocean_model', 'cfc12_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC12 into the ocean ', 'mol m-2 s-1', & + conversion= US%m_to_Z**2*US%s_to_T,& + cmor_field_name='fgcfc12', & + cmor_long_name='Surface Downward CFC12 Flux', & + cmor_standard_name='surface_downward_cfc12_flux') + + handles%id_ice_fraction = register_diag_field('ocean_model', 'ice_fraction', diag%axesT1, Time, & + 'Fraction of cell area covered by sea ice', 'm2 m-2') + + handles%id_u10_sqr = register_diag_field('ocean_model', 'u10_sqr', diag%axesT1, Time, & + 'Wind magnitude at 10m, squared', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + endif + endif + handles%id_psurf = register_diag_field('ocean_model', 'p_surf', diag%axesT1, Time, & 'Pressure at ice-ocean or atmosphere-ocean interface', & 'Pa', conversion=US%RL2_T2_to_Pa, cmor_field_name='pso', & @@ -2834,6 +2880,19 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_netFWGlobalScl > 0) & call post_data(handles%id_netFWGlobalScl, fluxes%netFWGlobalScl, diag) + ! post diagnostics related to cfcs ==================================== + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc11_flux)) & + call post_data(handles%id_cfc11, fluxes%cfc11_flux, diag) + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc12_flux)) & + call post_data(handles%id_cfc12, fluxes%cfc12_flux, diag) + + if ((handles%id_ice_fraction > 0) .and. associated(fluxes%ice_fraction)) & + call post_data(handles%id_ice_fraction, fluxes%ice_fraction, diag) + + if ((handles%id_u10_sqr > 0) .and. associated(fluxes%u10_sqr)) & + call post_data(handles%id_u10_sqr, fluxes%u10_sqr, diag) ! remaining boundary terms ================================================== @@ -2872,7 +2931,7 @@ end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & - shelf, iceberg, salt, fix_accum_bug) + shelf, iceberg, salt, fix_accum_bug, cfc) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields logical, optional, intent(in) :: water !< If present and true, allocate water fluxes @@ -2884,6 +2943,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: salt !< If present and true, allocate salt fluxes logical, optional, intent(in) :: fix_accum_bug !< If present and true, avoid using a bug in !! accumulation of ustar_gustless + logical, optional, intent(in) :: cfc !< If present and true, allocate cfc fluxes ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -2939,6 +2999,12 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%area_berg,isd,ied,jsd,jed, iceberg) call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) + !These fields should only on allocated when USE_CFC_CAP is activated. + call myAlloc(fluxes%cfc11_flux,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%cfc12_flux,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) + if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group @@ -3185,6 +3251,10 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%ustar_berg)) deallocate(fluxes%ustar_berg) if (associated(fluxes%area_berg)) deallocate(fluxes%area_berg) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) + if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) + if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) + if (associated(fluxes%cfc11_flux)) deallocate(fluxes%cfc11_flux) + if (associated(fluxes%cfc12_flux)) deallocate(fluxes%cfc12_flux) call coupler_type_destructor(fluxes%tr_fluxes) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 86493aad93..08f8dea634 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -11,7 +11,7 @@ module MOM_unit_tests use MOM_diag_vkernels, only : diag_vkernels_unit_tests use MOM_random, only : random_unit_tests use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests - +use MOM_CFC_cap, only : CFC_cap_unit_tests implicit none ; private public unit_tests @@ -41,6 +41,8 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: random_unit_tests FAILED") if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: near_boundary_unit_tests FAILED") + if (CFC_cap_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: CFC_cap_unit_tests FAILED") endif end subroutine unit_tests diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 886ee77510..2ea8943e67 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -42,6 +42,8 @@ module MOM_variables SST, & !< The sea surface temperature [degC]. SSS, & !< The sea surface salinity [ppt ~> psu or gSalt/kg]. sfc_density, & !< The mixed layer density [R ~> kg m-3]. + sfc_cfc11, & !< Sea surface concentration of CFC11 [mol kg-1]. + sfc_cfc12, & !< Sea surface concentration of CFC12 [mol kg-1]. Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. @@ -300,7 +302,8 @@ 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, use_meltpot, use_iceshelves, omit_frazil) + gas_fields_ocn, use_meltpot, use_iceshelves, & + omit_frazil, use_cfcs) 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. @@ -313,13 +316,14 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! 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, optional, intent(in) :: use_cfcs !< If true, allocate the space for cfcs logical, optional, intent(in) :: use_iceshelves !< If true, allocate the space for the stresses !! under ice shelves. logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to !! pass frazil fluxes to the coupler ! local variables - logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil + logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_cfcs integer :: is, ie, js, je, isd, ied, jsd, jed integer :: isdB, iedB, jsdB, jedB @@ -330,6 +334,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 + alloc_cfcs = .false. ; if (present(use_cfcs)) alloc_cfcs = use_cfcs alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil @@ -353,6 +358,11 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & allocate(sfc_state%melt_potential(isd:ied,jsd:jed)) ; sfc_state%melt_potential(:,:) = 0.0 endif + if (alloc_cfcs) then + allocate(sfc_state%sfc_cfc11(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc11(:,:) = 0.0 + allocate(sfc_state%sfc_cfc12(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc12(:,:) = 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 @@ -396,7 +406,8 @@ subroutine deallocate_surface_state(sfc_state) if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat) if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt) if (allocated(sfc_state%salt_deficit)) deallocate(sfc_state%salt_deficit) - + if (allocated(sfc_state%sfc_cfc11)) deallocate(sfc_state%sfc_cfc11) + if (allocated(sfc_state%sfc_cfc12)) deallocate(sfc_state%sfc_cfc12) call coupler_type_destructor(sfc_state%tr_fields) sfc_state%arrays_allocated = .false. diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 new file mode 100644 index 0000000000..58ba419c62 --- /dev/null +++ b/src/tracer/MOM_CFC_cap.F90 @@ -0,0 +1,704 @@ +!> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover +!! provided via cap (only NUOPC cap is implemented so far). +module MOM_CFC_cap + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_hor_index, only : hor_index_type +use MOM_grid, only : ocean_grid_type +use MOM_io, only : file_exists, MOM_read_data, slasher +use MOM_io, only : vardesc, var_desc, query_vardesc, stdout +use MOM_open_boundary, only : ocean_OBC_type +use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_time_manager, only : time_type +use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public register_CFC_cap, initialize_CFC_cap, CFC_cap_unit_tests +public CFC_cap_column_physics, CFC_cap_surface_state, CFC_cap_fluxes +public CFC_cap_stock, CFC_cap_end + +integer, parameter :: NTR = 2 !< the number of tracers in this module. + +!> The control structure for the CFC_cap tracer package +type, public :: CFC_cap_CS ; private + character(len=200) :: IC_file !< The file in which the CFC initial values can + !! be found, or an empty string for internal initilaization. + logical :: Z_IC_file !< If true, the IC_file is in Z-space. The default is false. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the MOM6 tracer registry + real, pointer, dimension(:,:,:) :: & + CFC11 => NULL(), & !< The CFC11 concentration [mol kg-1]. + CFC12 => NULL() !< The CFC12 concentration [mol kg-1]. + ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12. + real :: CFC11_IC_val = 0.0 !< The initial value assigned to CFC11 [mol kg-1]. + real :: CFC12_IC_val = 0.0 !< The initial value assigned to CFC12 [mol kg-1]. + real :: CFC11_land_val = -1.0 !< The value of CFC11 used where land is masked out [mol kg-1]. + real :: CFC12_land_val = -1.0 !< The value of CFC12 used where land is masked out [mol kg-1]. + logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code + !! if they are not found in the restart files. + character(len=16) :: CFC11_name !< CFC11 variable name + character(len=16) :: CFC12_name !< CFC12 variable name + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate + !! the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Model restart control structure + + ! The following vardesc types contain a package of metadata about each tracer. + type(vardesc) :: CFC11_desc !< A set of metadata for the CFC11 tracer + type(vardesc) :: CFC12_desc !< A set of metadata for the CFC12 tracer + !>@{ Diagnostic IDs + integer :: id_cfc11_cmor = -1, id_cfc12_cmor = -1 + !>@} +end type CFC_cap_CS + +contains + +!> Register the CFCs to be used with MOM and read the parameters +!! that are used with this tracer package +function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(CFC_cap_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module. + type(tracer_registry_type), & + pointer :: tr_Reg !< A pointer to the tracer registry. + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. + + ! Local variables + character(len=40) :: mdl = "MOM_CFC_cap" ! This module's name. + character(len=200) :: inputdir ! The directory where NetCDF input files are. + ! This include declares and sets the variable "version". +#include "version_variable.h" + real, dimension(:,:,:), pointer :: tr_ptr => NULL() + real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients + real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and + real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers. + character(len=48) :: flux_units ! The units for tracer fluxes. + character(len=48) :: dummy ! Dummy variable to store params that need to be logged here. + logical :: register_CFC_cap + integer :: isd, ied, jsd, jed, nz, m + + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(WARNING, "register_CFC_cap called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "CFC_IC_FILE", CS%IC_file, & + "The file in which the CFC initial values can be "//& + "found, or an empty string for internal initialization.", & + default=" ") + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + ! Add the directory if CS%IC_file is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file) + endif + call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", CS%Z_IC_file, & + "If true, CFC_IC_FILE is in depth space, not layer space", & + default=.false.) + call get_param(param_file, mdl, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + + ! the following params are not used in this module. Instead, they are used in + ! the cap but are logged here to keep all the CFC cap params together. + call get_param(param_file, mdl, "CFC_BC_FILE", dummy, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ") + if ((len_trim(dummy) > 0) .and. (scan(dummy,'/') == 0)) then + call get_param(param_file, mdl, "CFC11_VARIABLE", dummy, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11") + call get_param(param_file, mdl, "CFC12_VARIABLE", dummy, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12") + endif + + ! The following vardesc types contain a package of metadata about each tracer, + ! including, the name; units; longname; and grid information. + CS%CFC11_name = "CFC_11" ; CS%CFC12_name = "CFC_12" + CS%CFC11_desc = var_desc(CS%CFC11_name,"mol kg-1","Moles Per Unit Mass of CFC-11 in sea water", caller=mdl) + CS%CFC12_desc = var_desc(CS%CFC12_name,"mol kg-1","Moles Per Unit Mass of CFC-12 in sea water", caller=mdl) + if (GV%Boussinesq) then ; flux_units = "mol s-1" + else ; flux_units = "mol m-3 kg s-1" ; endif + + allocate(CS%CFC11(isd:ied,jsd:jed,nz)) ; CS%CFC11(:,:,:) = 0.0 + allocate(CS%CFC12(isd:ied,jsd:jed,nz)) ; CS%CFC12(:,:,:) = 0.0 + + ! This pointer assignment is needed to force the compiler not to do a copy in + ! the registration calls. Curses on the designers and implementers of F90. + tr_ptr => CS%CFC11 + ! Register CFC11 for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC11_desc, registry_diags=.true., & + flux_units=flux_units, & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + ! Do the same for CFC12 + tr_ptr => CS%CFC12 + call register_tracer(tr_ptr, Tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC12_desc, registry_diags=.true., & + flux_units=flux_units, & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + register_CFC_cap = .true. + +end function register_CFC_cap + +!> Initialize the CFC tracer fields and set up the tracer output. +subroutine initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate + !! diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type + !! specifies whether, where, and what + !! open boundary conditions are used. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: from_file = .false. + + if (.not.associated(CS)) return + + CS%Time => day + CS%diag => diag + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC11, CS%CFC11_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC11, CS%CFC11_name, CS%CFC11_land_val, & + CS%CFC11_IC_val, G, GV, US, CS) + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC12, CS%CFC12_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC12, CS%CFC12_name, CS%CFC12_land_val, & + CS%CFC12_IC_val, G, GV, US, CS) + + + ! cmor diagnostics + ! CFC11 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/42625c97b8fe75124a345962c4430982.html + CS%id_cfc11_cmor = register_diag_field('ocean_model', 'cfc11', diag%axesTL, day, & + 'Mole Concentration of CFC11 in Sea Water', 'mol m-3') + ! CFC12 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/3ab8e10027d7014f18f9391890369235.html + CS%id_cfc12_cmor = register_diag_field('ocean_model', 'cfc12', diag%axesTL, day, & + 'Mole Concentration of CFC12 in Sea Water', 'mol m-3') + + if (associated(OBC)) then + ! Steal from updated DOME in the fullness of time. + ! GMM: TODO this must be coded if we intend to use this module in regional applications + endif + +end subroutine initialize_CFC_cap + +!>This subroutine initializes a tracer array. +subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, GV, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: tr !< The tracer concentration array + character(len=*), intent(in) :: name !< The tracer name + real, intent(in) :: land_val !< A value the tracer takes over land + real, intent(in) :: IC_val !< The initial condition value for the tracer + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: OK + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file, G%Domain)) & + call MOM_error(FATAL, "initialize_CFC_cap: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr, h, CS%IC_file, name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr, h, CS%IC_file, trim(name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_CFC_cap: "//& + "Unable to read "//trim(name)//" from "//& + trim(CS%IC_file)//".") + endif + else + call MOM_read_data(CS%IC_file, trim(name), tr, G%Domain) + endif + else + do k=1,nz ; do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) < 0.5) then + tr(i,j,k) = land_val + else + tr(i,j,k) = IC_val + endif + enddo ; enddo ; enddo + endif + +end subroutine init_tracer_CFC + +!> Applies diapycnal diffusion, souces and sinks and any other column +!! tracer physics to the CFC cap tracers. CFCs are relatively simple, +!! as they are passive tracers with only a surface flux as a source. +subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, & + evap_CFL_limit, minimum_forcing_depth) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes!< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [H ~> m or kg m-2] + + ! The arguments to this subroutine are redundant in that + ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + + ! Local variables + real, pointer, dimension(:,:,:) :: CFC11 => NULL(), CFC12 => NULL() + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] + real :: scale_factor ! convert from [Conc. m s-1] to [Conc. kg m-2 T-1] + integer :: i, j, k, m, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + ! CFC_flux **unscaled** units are [mol m-2 s-1] which is the same as [CU kg m-2 s-1], + ! where CU = mol kg-1. However, CFC_flux has been scaled already. + ! CFC_flux **scaled** units are [CU L T-1 R]. We want [CU kg m-2 T-1] because + ! these units are what tracer_vertdiff needs. + ! Therefore, we need to convert [L R] to [kg m-2] using scale_factor. + scale_factor = US%L_to_Z * US%RZ_to_kg_m2 ! this will give [CU kg m-2 T-1], which is what verdiff wants + + if (.not.associated(CS)) return + + CFC11 => CS%CFC11 ; CFC12 => CS%CFC12 + + ! Use a tridiagonal solver to determine the concentrations after the + ! surface source is applied and diapycnal advection and diffusion occurs. + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC11, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC12, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + else + call tracer_vertdiff(h_old, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + call tracer_vertdiff(h_old, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + endif + + ! If needed, write out any desired diagnostics from tracer sources & sinks here. + if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, CFC11*GV%Rho0, CS%diag) + if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, CFC12*GV%Rho0, CS%diag) + +end subroutine CFC_cap_column_physics + +!> Calculates the mass-weighted integral of all tracer stocks, +!! returning the number of stocks it has calculated. If the stock_index +!! is present, only the stock corresponding to that coded index is returned. +function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of each + !! tracer, in kg times concentration units [kg conc]. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< The coded index of a specific + !! stock being sought. + integer :: CFC_cap_stock !< The number of stocks calculated here. + + ! Local variables + real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or nondim] + real :: mass ! The cell volume or mass [H L2 ~> m3 or kg] + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + CFC_cap_stock = 0 + if (.not.associated(CS)) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + call query_vardesc(CS%CFC11_desc, name=names(1), units=units(1), caller="CFC_cap_stock") + call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="CFC_cap_stock") + units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg" + + stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2 + stocks(1) = 0.0 ; stocks(2) = 0.0 + do k=1,nz ; do j=js,je ; do i=is,ie + mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k) + stocks(1) = stocks(1) + CS%CFC11(i,j,k) * mass + stocks(2) = stocks(2) + CS%CFC12(i,j,k) * mass + enddo ; enddo ; enddo + stocks(1) = stock_scale * stocks(1) + stocks(2) = stock_scale * stocks(2) + + CFC_cap_stock = 2 + +end function CFC_cap_stock + +!> Extracts the ocean surface CFC concentrations and copies them to sfc_state. +subroutine CFC_cap_surface_state(sfc_state, G, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + type(CFC_cap_CS), pointer :: CS!< The control structure returned by a previous + !! call to register_CFC_cap. + + ! Local variables + integer :: i, j, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + if (.not.associated(CS)) return + + do j=js,je ; do i=is,ie + sfc_state%sfc_cfc11(i,j) = CS%CFC11(i,j,1) + sfc_state%sfc_cfc12(i,j) = CS%CFC12(i,j,1) + enddo ; enddo + +end subroutine CFC_cap_surface_state + +!> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM +!! concentration, and calculating the solubility, Schmidt number, and gas exchange. +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc12_atm) + type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. + type(surface), intent(in ) :: sfc_state !< A structure containing fields + !! that describe the surface state of the ocean. + type(forcing), intent(inout) :: fluxes !< A structure containing pointers + !! to thermodynamic and tracer forcing fields. Unused fields + !! have NULL ptrs. + real, intent(in ) :: Rho0 !< The mean ocean density [kg m-3] + type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the + !! CFC's concentration in the atmosphere. + integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. + integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: & + CFC11_Csurf, & ! The CFC-11 surface concentrations times the Schmidt number term [mol kg-1]. + CFC12_Csurf, & ! The CFC-12 surface concentrations times the Schmidt number term [mol kg-1]. + CFC11_alpha, & ! The CFC-11 solubility [mol kg-1 atm-1]. + CFC12_alpha, & ! The CFC-12 solubility [mol kg-1 atm-1]. + kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [m s-1]. + cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) + ! [mol kg-1]. + cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] + cfc12_atm !< CFC11 concentration in the atmopshere [pico mol/mol] + real :: ta ! Absolute sea surface temperature [hectoKelvin] + real :: sal ! Surface salinity [PSU]. + real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. + real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. + real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12. + real :: sc_no_term ! A term related to the Schmidt number. + real :: kw_coeff ! A coefficient used to scale the piston velocity [L T-1 ~> m s-1] + real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. + integer :: i, j, m, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! CFC11 ATM concentration + if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then + call time_interp_external(id_cfc11_atm, Time, cfc11_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc11_atm = cfc11_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc11_atm internally" //& + "has not been implemented yet.") + endif + + ! CFC12 ATM concentration + if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then + call time_interp_external(id_cfc12_atm, Time, cfc12_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc12_atm = cfc12_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc12_atm internally" //& + "has not been implemented yet.") + endif + + do j=js,je ; do i=is,ie + ! ta in hectoKelvin + ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) + sal = sfc_state%SSS(i,j) + + ! Calculate solubilities + call get_solubility(alpha_11, alpha_12, ta, sal , G%mask2dT(i,j)) + + ! Calculate Schmidt numbers using coefficients given by + ! Wanninkhof (2014); doi:10.4319/lom.2014.12.351. + call comp_CFC_schmidt(sfc_state%SST(i,j), sc_11, sc_12) + sc_no_term = sqrt(660.0 / sc_11) + + CFC11_alpha(i,j) = alpha_11 * sc_no_term + CFC11_Csurf(i,j) = sfc_state%sfc_CFC11(i,j) * sc_no_term + sc_no_term = sqrt(660.0 / sc_12) + CFC12_alpha(i,j) = alpha_12 * sc_no_term + CFC12_Csurf(i,j) = sfc_state%sfc_CFC12(i,j) * sc_no_term + + !--------------------------------------------------------------------- + ! Gas exchange/piston velocity parameter + !--------------------------------------------------------------------- + ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 + ! 6.97e-07 is used to convert from cm/hr to m/s, kw = m/s [L/T] + kw_coeff = 6.97e-07 * G%US%m_to_L * G%US%T_to_s + + kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) + + ! air concentrations and cfcs BC's fluxes + ! CFC flux units: mol kg-1 s-1 kg m-2 + cair(i,j) = pa_to_atm * CFC11_alpha(i,j) * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC11_Csurf(i,j)) * Rho0 + cair(i,j) = pa_to_atm * CFC12_alpha(i,j) * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC12_Csurf(i,j)) * Rho0 + enddo ; enddo + +end subroutine CFC_cap_fluxes + +!> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32. +subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask) + real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1] + real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1] + real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin] + real, intent(in ) :: sal !< Surface salinity [PSU]. + real, intent(in ) :: mask !< ocean mask + + ! Local variables + real :: d11_dflt(4), d12_dflt(4) ! values of the various coefficients + real :: e11_dflt(3), e12_dflt(3) ! in the expressions for the solubility + real :: d1_11, d1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [nondim] + real :: d2_11, d2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-1] + real :: d3_11, d3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [log(hectoKelvin)-1] + real :: d4_11, d4_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-2] + real :: e1_11, e1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1] + real :: e2_11, e2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1 hectoKelvin-1] + real :: e3_11, e3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-2 hectoKelvin-2] + real :: factor ! factor to use in the solubility conversion + + !----------------------------------------------------------------------- + ! Solubility coefficients for alpha in mol/(kg atm) for CFC11 (_11) and CFC12 (_12) + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + !----------------------------------------------------------------------- + d11_dflt(:) = (/ -232.0411, 322.5546, 120.4956, -1.39165 /) + e11_dflt(:) = (/ -0.146531, 0.093621, -0.0160693 /) + d12_dflt(:) = (/ -220.2120, 301.8695, 114.8533, -1.39165 /) + e12_dflt(:) = (/ -0.147718, 0.093175, -0.0157340 /) + + d1_11 = d11_dflt(1) + d2_11 = d11_dflt(2) + d3_11 = d11_dflt(3) + d4_11 = d11_dflt(4) + + e1_11 = e11_dflt(1) + e2_11 = e11_dflt(2) + e3_11 = e11_dflt(3) + + d1_12 = d12_dflt(1) + d2_12 = d12_dflt(2) + d3_12 = d12_dflt(3) + d4_12 = d12_dflt(4) + + e1_12 = e12_dflt(1) + e2_12 = e12_dflt(2) + e3_12 = e12_dflt(3) + + ! Calculate solubilities using Warner and Weiss (1985) DSR, vol 32. + ! The following is from Eq. 9 in Warner and Weiss (1985) + ! The factor 1.0e+03 is for the conversion from mol/(l * atm) to mol/(m3 * atm) + ! The factor 1.e-09 converts from mol/(l * atm) to mol/(m3 * pptv) + factor = 1.0 + alpha_11 = exp(d1_11 + d2_11/ta + d3_11*log(ta) + d4_11*ta**2 +& + sal * ((e3_11 * ta + e2_11) * ta + e1_11)) * & + factor * mask + alpha_12 = exp(d1_12 + d2_12/ta + d3_12*log(ta) + d4_12*ta**2 +& + sal * ((e3_12 * ta + e2_12) * ta + e1_12)) * & + factor * mask + +end subroutine get_solubility + + +!> Compute Schmidt numbers of CFCs following Wanninkhof (2014); doi:10.4319/lom.2014.12.351 +!! Range of validity of fit is -2:40. +subroutine comp_CFC_schmidt(sst_in, cfc11_sc, cfc12_sc) + real, intent(in) :: sst_in !< The sea surface temperature [degC]. + real, intent(inout) :: cfc11_sc !< Schmidt number of CFC11 [nondim]. + real, intent(inout) :: cfc12_sc !< Schmidt number of CFC12 [nondim]. + + !local variables + real , parameter :: a_11 = 3579.2 + real , parameter :: b_11 = -222.63 + real , parameter :: c_11 = 7.5749 + real , parameter :: d_11 = -0.14595 + real , parameter :: e_11 = 0.0011874 + real , parameter :: a_12 = 3828.1 + real , parameter :: b_12 = -249.86 + real , parameter :: c_12 = 8.7603 + real , parameter :: d_12 = -0.1716 + real , parameter :: e_12 = 0.001408 + real :: sst + + + ! clip SST to avoid bad values + sst = MAX(-2.0, MIN(40.0, sst_in)) + cfc11_sc = a_11 + sst * (b_11 + sst * (c_11 + sst * (d_11 + sst * e_11))) + cfc12_sc = a_12 + sst * (b_12 + sst * (c_12 + sst * (d_12 + sst * e_12))) + +end subroutine comp_CFC_schmidt + +!> Deallocate any memory associated with the CFC cap tracer package +subroutine CFC_cap_end(CS) + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + ! local variables + integer :: m + + if (associated(CS)) then + if (associated(CS%CFC11)) deallocate(CS%CFC11) + if (associated(CS%CFC12)) deallocate(CS%CFC12) + + deallocate(CS) + endif +end subroutine CFC_cap_end + +!> Unit tests for the CFC cap module. +logical function CFC_cap_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, output additional + !! information for debugging unit tests + + ! Local variables + real :: dummy1, dummy2, ta, sal + character(len=120) :: test_name ! Title of the unit test + + CFC_cap_unit_tests = .false. + write(stdout,*) '==== MOM_CFC_cap =======================' + + ! test comp_CFC_schmidt, Table 1 in Wanninkhof (2014); doi:10.4319/lom.2014.12.351 + test_name = 'Schmidt number calculation' + call comp_CFC_schmidt(20.0, dummy1, dummy2) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 1179.0, 0.5) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 1188.0, 0.5) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)') "Passed "//test_name + + test_name = 'Solubility function, SST = 1.0 C, and SSS = 10 psu' + ta = max(0.01, (1.0 + 273.15) * 0.01); sal = 10. + ! cfc1 = 3.238 10-2 mol kg-1 atm-1 + ! cfc2 = 7.943 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 3.238e-2, 5.0e-6) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 7.943e-3, 5.0e-6) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + + test_name = 'Solubility function, SST = 20.0 C, and SSS = 35 psu' + ta = max(0.01, (20.0 + 273.15) * 0.01); sal = 35. + ! cfc1 = 0.881 10-2 mol kg-1 atm-1 + ! cfc2 = 2.446 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 8.8145e-3, 5.0e-8) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 2.4462e-3, 5.0e-8) + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + +end function CFC_cap_unit_tests + +!> Test that ans and calc are approximately equal by computing the difference +!! and comparing it against limit. +logical function compare_values(verbose, test_name, calc, ans, limit) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + real, intent(in) :: calc !< computed value + real, intent(in) :: ans !< correct value + real, intent(in) :: limit !< value above which test fails + + ! Local variables + real :: diff + + diff = ans - calc + + compare_values = .false. + if (diff > limit ) then + compare_values = .true. + write(stdout,*) "CFC_cap_unit_tests, UNIT TEST FAILED: ", test_name + write(stdout,10) calc, ans + elseif (verbose) then + write(stdout,10) calc, ans + endif + +10 format("calc=",f20.16," ans",f20.16) +end function compare_values + +!> \namespace mom_CFC_cap +!! +!! This module contains the code that is needed to simulate +!! CFC-11 and CFC-12 using atmospheric and sea ice variables +!! provided via cap (only NUOPC cap is implemented so far). + +end module MOM_CFC_cap diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 26ef197ae2..439e9b5396 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -42,6 +42,9 @@ module MOM_tracer_flow_control use MOM_OCMIP2_CFC, only : register_OCMIP2_CFC, initialize_OCMIP2_CFC, flux_init_OCMIP2_CFC use MOM_OCMIP2_CFC, only : OCMIP2_CFC_column_physics, OCMIP2_CFC_surface_state use MOM_OCMIP2_CFC, only : OCMIP2_CFC_stock, OCMIP2_CFC_end, OCMIP2_CFC_CS +use MOM_CFC_cap, only : register_CFC_cap, initialize_CFC_cap +use MOM_CFC_cap, only : CFC_cap_column_physics, CFC_cap_surface_state +use MOM_CFC_cap, only : CFC_cap_stock, CFC_cap_end, CFC_cap_CS use oil_tracer, only : register_oil_tracer, initialize_oil_tracer use oil_tracer, only : oil_tracer_column_physics, oil_tracer_surface_state use oil_tracer, only : oil_stock, oil_tracer_end, oil_tracer_CS @@ -80,6 +83,7 @@ module MOM_tracer_flow_control logical :: use_oil = .false. !< If true, use the oil tracer package logical :: use_advection_test_tracer = .false. !< If true, use the advection_test_tracer package logical :: use_OCMIP2_CFC = .false. !< If true, use the OCMIP2_CFC tracer package + logical :: use_CFC_cap = .false. !< If true, use the CFC_cap tracer package logical :: use_MOM_generic_tracer = .false. !< If true, use the MOM_generic_tracer packages logical :: use_pseudo_salt_tracer = .false. !< If true, use the psuedo_salt tracer package logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package @@ -94,6 +98,7 @@ module MOM_tracer_flow_control type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL() type(advection_test_tracer_CS), pointer :: advection_test_tracer_CSp => NULL() type(OCMIP2_CFC_CS), pointer :: OCMIP2_CFC_CSp => NULL() + type(CFC_cap_CS), pointer :: CFC_cap_CSp => NULL() type(MOM_generic_tracer_CS), pointer :: MOM_generic_tracer_CSp => NULL() type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL() type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL() @@ -114,7 +119,7 @@ subroutine call_tracer_flux_init(verbosity) type(param_file_type) :: param_file ! A structure to parse for run-time parameters character(len=40) :: mdl = "call_tracer_flux_init" ! This module's name. - logical :: use_OCMIP_CFCs, use_MOM_generic_tracer + logical :: use_OCMIP_CFCs, use_MOM_generic_tracer, use_CFC_caps ! Determine which tracer routines with tracer fluxes are to be called. Note ! that not every tracer package is required to have a flux_init call. @@ -193,6 +198,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_OCMIP2_CFC", CS%use_OCMIP2_CFC, & "If true, use the MOM_OCMIP2_CFC tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC_cap, & + "If true, use the MOM_CFC_cap tracer package.", & + default=.false.) call get_param(param_file, mdl, "USE_generic_tracer", CS%use_MOM_generic_tracer, & "If true and _USE_GENERIC_TRACER is defined as a "//& "preprocessor macro, use the MOM_generic_tracer packages.", & @@ -237,6 +245,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_OCMIP2_CFC) CS%use_OCMIP2_CFC = & register_OCMIP2_CFC(HI, GV, param_file, CS%OCMIP2_CFC_CSp, & tr_Reg, restart_CS) + if (CS%use_CFC_cap) CS%use_CFC_cap = & + register_CFC_cap(HI, GV, param_file, CS%CFC_cap_CSp, & + tr_Reg, restart_CS) if (CS%use_MOM_generic_tracer) CS%use_MOM_generic_tracer = & register_MOM_generic_tracer(HI, GV, param_file, CS%MOM_generic_tracer_CSp, & tr_Reg, restart_CS) @@ -317,6 +328,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag if (CS%use_OCMIP2_CFC) & call initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS%OCMIP2_CFC_CSp, & sponge_CSp) + if (CS%use_CFC_cap) & + call initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS%CFC_cap_CSp) + if (CS%use_MOM_generic_tracer) & call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp) @@ -462,6 +476,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, G, GV, US, CS%OCMIP2_CFC_CSp, & evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -516,6 +535,9 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -624,6 +646,12 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni set_pkg_name, max_ns, ns_tot, stock_names, stock_units) endif + if (CS%use_CFC_cap) then + ns = CFC_cap_stock(h, values, G, GV, CS%CFC_cap_CSp, names, units, stock_index) + call store_stocks("MOM_CFC_cap", ns, names, units, values, index, stock_values, & + set_pkg_name, max_ns, ns_tot, stock_names, stock_units) + endif + if (CS%use_advection_test_tracer) then ns = advection_test_stock( h, values, G, GV, CS%advection_test_tracer_CSp, & names, units, stock_index ) @@ -753,6 +781,8 @@ subroutine call_tracer_surface_state(sfc_state, h, G, GV, CS) call advection_test_tracer_surface_state(sfc_state, h, G, GV, CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_surface_state(sfc_state, h, G, GV, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_surface_state(sfc_state, G, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) & call MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS%MOM_generic_tracer_CSp) @@ -772,6 +802,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) call advection_test_tracer_end(CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) call OCMIP2_CFC_end(CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) call CFC_cap_end(CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) call end_MOM_generic_tracer(CS%MOM_generic_tracer_CSp) if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp) if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp) From f45659d47d48ef82fd2c26d7eb94ad49928d037f Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Mon, 14 Jun 2021 20:22:57 -0600 Subject: [PATCH 05/28] introduce the lamult import from ww3 in nuopc cap and add a diag --- config_src/drivers/nuopc_cap/mom_cap.F90 | 17 +++---- .../drivers/nuopc_cap/mom_cap_methods.F90 | 8 ++++ .../nuopc_cap/mom_ocean_model_nuopc.F90 | 7 +-- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 19 ++++++-- src/core/MOM_forcing_type.F90 | 46 +++++++++++++++---- 5 files changed, 71 insertions(+), 26 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 8d383c3d3f..e0aaff8534 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -697,7 +697,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (ocean_state%use_waves) then if (ocean_state%wave_method == "VR12-MA") then - continue + allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) ) + Ice_ocean_boundary%lamult = 0.0 else Ice_ocean_boundary%num_stk_bands=ocean_state%Waves%NumBands allocate ( Ice_ocean_boundary% ustk0 (isc:iec,jsc:jec), & @@ -722,15 +723,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsFrOcn_num, fldsFrOcn, trim(scalar_field_name), "will_provide") end if - if (cesm_coupled) then - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes" , "will provide") - !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_hstokes" , "will provide") - else - !call fld_list_add(fldsToOcn_num, fldsToOcn, "mass_of_overlying_sea_ice" , "will provide") - !call fld_list_add(fldsFrOcn_num, fldsFrOcn, "sea_lev" , "will provide") - endif !--------- import fields ------------- call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_salt_rate" , "will provide") ! from ice @@ -755,7 +747,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") if (ocean_state%use_waves) then if (ocean_state%wave_method == "VR12-MA") then - continue + call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide") + !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_ustokes" , "will provide") + !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_vstokes" , "will provide") + !call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_hstokes" , "will provide") else if (Ice_ocean_boundary%num_stk_bands > 3) then call MOM_error(FATAL, "Number of Stokes Bands > 3, NUOPC cap not set up for this") diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 0a4e12c647..085dbe818c 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -256,6 +256,14 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- + ! Langmuir enhancement factor + !---- + ice_ocean_boundary%lamult (:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Sw_lamult', & + isc, iec, jsc, jec, ice_ocean_boundary%lamult, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- ! Partitioned Stokes Drift Components !---- diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 7e7716d386..a07b24e8d8 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -367,13 +367,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i use_melt_pot=.false. endif + call get_param(param_file, mdl, "USE_WAVES", OS%use_waves, & + "If true, enables surface wave modules.", default=.false.) + ! 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, use_meltpot=use_melt_pot) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & - OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) + OS%forcing_CSp, OS%restore_salinity, OS%restore_temp, OS%use_waves) if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & @@ -385,8 +388,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call allocate_forcing_type(OS%grid, OS%fluxes, shelf=.true.) endif - call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & - "If true, enables surface wave modules.", default=.false.) if (OS%use_waves) then call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.) call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 7168823fbc..61d900f0b6 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -183,6 +183,7 @@ module MOM_surface_forcing_nuopc !! ice-shelves, expressed as a coefficient !! for divergence damping, as determined !! outside of the ocean model in [m3/s] + real, pointer, dimension(:,:) :: lamult => NULL() !< Langmuir enhancement factor [nondim] real, pointer, dimension(:,:) :: ustk0 => NULL() !< Surface Stokes drift, zonal [m/s] real, pointer, dimension(:,:) :: vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] real, pointer, dimension(:) :: stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m] @@ -688,8 +689,11 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0 if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0 - if ( associated(IOB%ustkb) ) & + if ( associated(IOB%lamult) ) then + call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=0) + elseif ( associated(IOB%ustkb) ) then call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands) + endif ! applied surface pressure from atmosphere and cryosphere if (CS%use_limited_P_SSH) then @@ -849,6 +853,13 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) endif ! endif for wind related fields ! wave to ocean coupling + if ( associated(IOB%lamult)) then + do j=js,je; do i=is,ie + forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + enddo ; enddo + call pass_var(forces%lamult, G%domain, halo=1 ) + endif + if ( associated(IOB%ustkb) ) then forces%stk_wavenumbers(:) = IOB%stk_wavenumbers @@ -1041,7 +1052,7 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & end subroutine forcing_save_restart !> Initialize the surface forcing, including setting parameters and allocating permanent memory. -subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp) +subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, restore_temp, use_waves) type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1054,6 +1065,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, !! restoring will be applied in this model. logical, optional, intent(in) :: restore_temp !< If present and true surface temperature !! restoring will be applied in this model. + logical, optional, intent(in) :: use_waves !< If present and true, use waves and activate + !! the corresponding wave forcing diagnostics ! Local variables real :: utide ! The RMS tidal velocity [Z T-1 ~> m s-1]. @@ -1321,7 +1334,7 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "If true, makes available diagnostics of fluxes from icebergs "//& "as seen by MOM6.", default=.false.) call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, & - use_berg_fluxes=iceberg_flux_diags) + use_berg_fluxes=iceberg_flux_diags, use_waves=use_waves) call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & "If true, allows flux adjustments to specified via the "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 682ad03397..6e2c0040c8 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -246,7 +246,8 @@ module MOM_forcing_type logical :: accumulate_rigidity = .false. !< If true, the rigidity due to various types of !! ice needs to be accumulated, and the rigidity explicitly !! reset to zero at the driver level when appropriate. - + real, pointer, dimension(:,:) :: & + lamult => NULL() !< Langmuir enhancement factor [nondim] real, pointer, dimension(:,:) :: & ustk0 => NULL(), & !< Surface Stokes drift, zonal [m/s] vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] @@ -356,6 +357,9 @@ module MOM_forcing_type ! Iceberg + Ice shelf diagnostic handles integer :: id_ustar_ice_cover = -1 integer :: id_frac_ice_cover = -1 + + ! wave forcing diagnostics handles. + integer :: id_lamult = -1 !>@} integer :: id_clock_forcing = -1 !< CPU clock id @@ -1239,13 +1243,14 @@ end subroutine forcing_SinglePointPrint !> Register members of the forcing type for diagnostics -subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_waves) type(time_type), intent(in) :: Time !< time type type(diag_ctrl), intent(inout) :: diag !< diagnostic control type type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< True if T/S are in use type(forcing_diags), intent(inout) :: handles !< handles for diagnostics logical, optional, intent(in) :: use_berg_fluxes !< If true, allow iceberg flux diagnostics + logical, optional, intent(in) :: use_waves !< If true, allow wave forcing diagnostics ! Clock for forcing diagnostics handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1895,6 +1900,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, handles%id_total_saltFluxAdded = register_scalar_field('ocean_model', 'total_salt_Flux_Added', & Time, diag, long_name='Area integrated surface salt flux due to restoring or flux adjustment', units='kg s-1') + !=============================================================== + ! wave forcing diagnostics + if (present(use_waves)) then + if (use_waves) then + handles%id_lamult = register_diag_field('ocean_model', 'lamult', & + diag%axesT1, Time, long_name='Langmuir enhancement factor received from WW3', units="nondim") + endif + endif end subroutine register_forcing_type_diags @@ -2266,6 +2279,10 @@ subroutine mech_forcing_diags(forces_in, dt, G, time_end, diag, handles) ! endif + ! wave forcing =============================================================== + if (handles%id_lamult > 0) & + call post_data(handles%id_lamult, forces%lamult, diag) + call disable_averaging(diag) if (turns /= 0) then @@ -3034,14 +3051,25 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & !These fields should only be allocated when waves call myAlloc(forces%ustk0,isd,ied,jsd,jed, waves) call myAlloc(forces%vstk0,isd,ied,jsd,jed, waves) - if (present(waves)) then; if (waves) then; if (.not.associated(forces%ustkb)) then - if (.not.(present(num_stk_bands))) call MOM_error(FATAL,"Requested to & + if (present(waves)) then; if (waves) then; + if (.not. present(num_stk_bands)) then + call MOM_error(FATAL,"Requested to & initialize with waves, but no waves are present.") - allocate(forces%stk_wavenumbers(num_stk_bands)) - forces%stk_wavenumbers(:) = 0.0 - allocate(forces%ustkb(isd:ied,jsd:jed,num_stk_bands)) - forces%ustkb(isd:ied,jsd:jed,:) = 0.0 - endif ; endif ; endif + endif + if (num_stk_bands==0) then + if (.not.associated(forces%lamult)) then + allocate(forces%lamult(isd:ied,jsd:jed)) + endif + else !num_stk_bands > 0 + if (.not.associated(forces%ustkb)) then + allocate(forces%stk_wavenumbers(num_stk_bands)) + forces%stk_wavenumbers(:) = 0.0 + allocate(forces%ustkb(isd:ied,jsd:jed,num_stk_bands)) + forces%ustkb(isd:ied,jsd:jed,:) = 0.0 + endif + endif + endif ; endif + if (present(waves)) then; if (waves) then; if (.not.associated(forces%vstkb)) then allocate(forces%vstkb(isd:ied,jsd:jed,num_stk_bands)) From b767a134310c4cc4c124cb393e6db61bdf31a6b3 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 6 Jul 2021 22:32:20 -0600 Subject: [PATCH 06/28] import ifrac and u10_sqr only when CFC module is active --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 1054234e6b..42ad69fb1a 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -184,7 +184,7 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: frunoff_hflx =>NULL() !< heat content of frozen runoff [W/m2] real, pointer, dimension(:,:) :: p =>NULL() !< pressure of overlying ice and atmosphere !< on ocean surface [Pa] - real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< mass of ice [nondim] + real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< fractional ice area [nondim] real, pointer, dimension(:,:) :: u10_sqr =>NULL() !< wind speed squared at 10m [m2/s2] real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] real, pointer, dimension(:,:) :: ice_rigidity =>NULL() !< rigidity of the sea ice, sea-ice and @@ -508,14 +508,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%seaice_melt)) & fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0) - ! sea ice fraction [nondim] - if (associated(IOB%ice_fraction)) & - fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) - - ! 10-m wind speed squared [m2/s2] - if (associated(IOB%u10_sqr)) & - fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) - fluxes%latent(i,j) = 0.0 ! notice minus sign since fprec is positive into the ocean if (associated(IOB%fprec)) then @@ -554,6 +546,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo + if (CS%use_CFC) then + do j=js,je ; do i=is,ie + ! sea ice fraction [nondim] + if (associated(IOB%ice_fraction)) & + fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) + ! 10-m wind speed squared [m2/s2] + if (associated(IOB%u10_sqr)) & + fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) + enddo ; enddo + endif + ! applied surface pressure from atmosphere and cryosphere if (associated(IOB%p)) then if (CS%max_p_surf >= 0.0) then @@ -1474,8 +1477,10 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) write(outunit,100) 'iobt%lrunoff ' , mpp_chksum( iobt%lrunoff ) write(outunit,100) 'iobt%frunoff ' , mpp_chksum( iobt%frunoff ) write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p ) - write(outunit,100) 'iobt%ice_fraction ' , mpp_chksum( iobt%ice_fraction ) - write(outunit,100) 'iobt%u10_sqr ' , mpp_chksum( iobt%u10_sqr ) + if (associated(iobt%ice_fraction)) & + write(outunit,100) 'iobt%ice_fraction ' , mpp_chksum( iobt%ice_fraction ) + if (associated(iobt%u10_sqr)) & + write(outunit,100) 'iobt%u10_sqr ' , mpp_chksum( iobt%u10_sqr ) if (associated(iobt%ustar_berg)) & write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & From 5f756ab7e2312039686fa97e2a642db819c4a2e1 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 18:48:39 -0600 Subject: [PATCH 07/28] Enable the export of ice fraction when wave coupling is on. The export of ice fraction was previously added for the newly added CFC module. --- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 1 + .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 18 +++++++----------- src/core/MOM_forcing_type.F90 | 6 +++++- 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 1af4185a78..d9e5e53fdd 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -395,6 +395,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif if (OS%use_waves) then + call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true.) call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.) call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) if (OS%wave_method == "VR12-MA") then diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index f1429bf867..4b03ccf743 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -545,18 +545,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + & fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j) - enddo ; enddo + ! sea ice fraction [nondim] + if (associated(IOB%ice_fraction) .and. associated(fluxes%ice_fraction)) & + fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) + ! 10-m wind speed squared [m2/s2] + if (associated(IOB%u10_sqr) .and. associated(fluxes%u10_sqr)) & + fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) - if (CS%use_CFC) then - do j=js,je ; do i=is,ie - ! sea ice fraction [nondim] - if (associated(IOB%ice_fraction)) & - fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) - ! 10-m wind speed squared [m2/s2] - if (associated(IOB%u10_sqr)) & - fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) - enddo ; enddo - endif + enddo ; enddo ! applied surface pressure from atmosphere and cryosphere if (associated(IOB%p)) then diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 69114d9339..cd4aa06edd 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -2948,7 +2948,7 @@ end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & - shelf, iceberg, salt, fix_accum_bug, cfc) + shelf, iceberg, salt, fix_accum_bug, cfc, waves) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields logical, optional, intent(in) :: water !< If present and true, allocate water fluxes @@ -2961,6 +2961,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: fix_accum_bug !< If present and true, avoid using a bug in !! accumulation of ustar_gustless logical, optional, intent(in) :: cfc !< If present and true, allocate cfc fluxes + logical, optional, intent(in) :: waves !< If present and true, allocate wave fields ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -3022,6 +3023,9 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) + !These fields should only on allocated when wave coupling is activated. + call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, waves) + if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group From 7d85ab2912ced10af8bcb10e18df3e0c6aec903a Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 19:49:09 -0600 Subject: [PATCH 08/28] set lamult to 1 below ice --- config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 4b03ccf743..228b47dc69 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -879,7 +879,11 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) ! wave to ocean coupling if ( associated(IOB%lamult)) then do j=js,je; do i=is,ie - forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + if (IOB%ice_fraction(i-i0,j-j0) <= 0.05 ) then + forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) + else + forces%lamult(i,j) = 1.0 + endif enddo ; enddo call pass_var(forces%lamult, G%domain, halo=1 ) endif From 997636e3cd42c1ccc2a7e088d46bce61339dc9db Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 8 Jul 2021 20:54:27 -0600 Subject: [PATCH 09/28] move lamult from forces to fluxes --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 24 +++++++++---------- src/core/MOM_forcing_type.F90 | 20 +++++++--------- 2 files changed, 21 insertions(+), 23 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 228b47dc69..4743e0d268 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -554,6 +554,18 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo + ! wave to ocean coupling + if ( associated(IOB%lamult)) then + do j=js,je; do i=is,ie + if (IOB%ice_fraction(i-i0,j-j0) <= 0.05 ) then + fluxes%lamult(i,j) = IOB%lamult(i-i0,j-j0) + else + fluxes%lamult(i,j) = 1.0 + endif + enddo ; enddo + call pass_var(fluxes%lamult, G%domain, halo=1 ) + endif + ! applied surface pressure from atmosphere and cryosphere if (associated(IOB%p)) then if (CS%max_p_surf >= 0.0) then @@ -876,18 +888,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) endif ! endif for wind related fields - ! wave to ocean coupling - if ( associated(IOB%lamult)) then - do j=js,je; do i=is,ie - if (IOB%ice_fraction(i-i0,j-j0) <= 0.05 ) then - forces%lamult(i,j) = IOB%lamult(i-i0,j-j0) - else - forces%lamult(i,j) = 1.0 - endif - enddo ; enddo - call pass_var(forces%lamult, G%domain, halo=1 ) - endif - if ( associated(IOB%ustkb) ) then forces%stk_wavenumbers(:) = IOB%stk_wavenumbers diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index cd4aa06edd..7022c8412c 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -190,6 +190,9 @@ module MOM_forcing_type ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] + real, pointer, dimension(:,:) :: & + lamult => NULL() !< Langmuir enhancement factor [nondim] + ! passive tracer surface fluxes type(coupler_2d_bc_type) :: tr_fluxes !< This structure contains arrays of !! of named fields used for passive tracer fluxes. @@ -253,8 +256,6 @@ module MOM_forcing_type logical :: accumulate_rigidity = .false. !< If true, the rigidity due to various types of !! ice needs to be accumulated, and the rigidity explicitly !! reset to zero at the driver level when appropriate. - real, pointer, dimension(:,:) :: & - lamult => NULL() !< Langmuir enhancement factor [nondim] real, pointer, dimension(:,:) :: & ustk0 => NULL(), & !< Surface Stokes drift, zonal [m/s] vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] @@ -2325,10 +2326,6 @@ subroutine mech_forcing_diags(forces_in, dt, G, time_end, diag, handles) ! endif - ! wave forcing =============================================================== - if (handles%id_lamult > 0) & - call post_data(handles%id_lamult, forces%lamult, diag) - call disable_averaging(diag) if (turns /= 0) then @@ -2934,6 +2931,10 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if ((handles%id_ustar_ice_cover > 0) .and. associated(fluxes%ustar_shelf)) & call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) + ! wave forcing =============================================================== + if (handles%id_lamult > 0) & + call post_data(handles%id_lamult, fluxes%lamult, diag) + ! endif ! query_averaging_enabled call disable_averaging(diag) @@ -3025,6 +3026,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & !These fields should only on allocated when wave coupling is activated. call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, waves) + call myAlloc(fluxes%lamult,isd,ied,jsd,jed, waves) if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group @@ -3126,11 +3128,7 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & call MOM_error(FATAL,"Requested to & initialize with waves, but no waves are present.") endif - if (num_stk_bands==0) then - if (.not.associated(forces%lamult)) then - allocate(forces%lamult(isd:ied,jsd:jed)) - endif - else !num_stk_bands > 0 + if (num_stk_bands > 0) then if (.not.associated(forces%ustkb)) then allocate(forces%stk_wavenumbers(num_stk_bands)) forces%stk_wavenumbers(:) = 0.0 From 762291114748d87ce0045ddc73c5a5b287dbcea1 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 23 Jul 2021 14:30:20 -0600 Subject: [PATCH 10/28] Changes the order of calculation In a previous commit, the expression tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) was change to tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) in 4 lines. This changed the order in which b1(i)*h_tr*tr(i,j,1) is evaluated, even if sfc_src(i,j) is 0.0. Because the CESM configurations use the passive tracer tridiagonal solver for T and S, the modification above changed answers in these configurations. To recover previous answers the expression must be tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) --- src/tracer/MOM_tracer_diabatic.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index cbfe2d7f0a..c1e39598cc 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -152,7 +152,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) @@ -190,7 +190,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & b_denom_1 = h_tr + ea(i,j,1) b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = h_tr * b1(i) - tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) @@ -351,7 +351,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) @@ -389,7 +389,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & b_denom_1 = h_tr + ent(i,j,1) b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = h_tr * b1(i) - tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) + tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + b1(i)*sfc_src(i,j) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) From a4ca3073c0cd8b7d3a106e4f3fbbe3760c3e3cd2 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 29 Jul 2021 15:01:26 +0000 Subject: [PATCH 11/28] remove white space and fix comment --- src/ocean_data_assim/MOM_oda_incupd.F90 | 2 +- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 12d18587bf..334597cda7 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -450,7 +450,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo endif enddo; enddo - + ! remap u to h_obs to get increment if (CS%uv_inc) then call pass_var(h, G%Domain) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 86a83c58a9..fb759b13b2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -111,7 +111,7 @@ module MOM_diabatic_driver !! domain. The exact location and properties of !! those sponges are set by calls to !! initialize_sponge and set_up_sponge_field. - logical :: use_oda_incupd !!< If True, DA incremental update is + logical :: use_oda_incupd !< If True, DA incremental update is !! applied everywhere logical :: use_geothermal !< If true, apply geothermal heating. logical :: use_int_tides !< If true, use the code that advances a separate set From cf46bc90bd04746e72ca3618e3305a2afc2b4e37 Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Thu, 5 Aug 2021 10:20:00 -0600 Subject: [PATCH 12/28] Update MOM_oda_incupd.F90 remove unused index bounds, and fix sum_h2 loop. --- src/ocean_data_assim/MOM_oda_incupd.F90 | 21 ++------------------- 1 file changed, 2 insertions(+), 19 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 334597cda7..d3199dcb74 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -65,22 +65,6 @@ module MOM_oda_incupd type, public :: oda_incupd_CS ; private integer :: nz !< The total number of layers. integer :: nz_data !< The total number of arbritary layers (used by older code). - integer :: isc !< The starting i-index of the computational domain at h. - integer :: iec !< The ending i-index of the computational domain at h. - integer :: jsc !< The starting j-index of the computational domain at h. - integer :: jec !< The ending j-index of the computational domain at h. - integer :: IscB !< The starting I-index of the computational domain at u/v. - integer :: IecB !< The ending I-index of the computational domain at u/v. - integer :: JscB !< The starting J-index of the computational domain at u/v. - integer :: JecB !< The ending J-index of the computational domain at h. - integer :: isd !< The starting i-index of the data domain at h. - integer :: ied !< The ending i-index of the data domain at h. - integer :: jsd !< The starting j-index of the data domain at h. - integer :: jed !< The ending j-index of the data domain at h. - integer :: IsdB !< The starting I-index of the data domain at u/v. - integer :: IedB !< The ending I-index of the data domain at u/v. - integer :: JsdB !< The starting J-index of the data domain at u/v. - integer :: JedB !< The ending J-index of the data domain at h. integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_oda_incupd_field @@ -224,9 +208,6 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res default=.true.) CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%isdB = G%isdB ; CS%iedB = G%iedB; CS%jsdB = G%jsdB ; CS%jedB = G%jedB ! increments on horizontal grid if (.not. CS%incupdDataOngrid) call MOM_error(FATAL,'initialize_oda_incupd: '// & @@ -421,6 +402,8 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) do k=1,nz_data sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) enddo ! get temperature From 141277181b2a58cb870dbb6ff1aa041ef786cade Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 5 Aug 2021 10:56:42 -0600 Subject: [PATCH 13/28] Minor fixes for the WW3 coupling option VR12-MA --- config_src/drivers/nuopc_cap/mom_cap_methods.F90 | 10 ++++++---- config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 | 3 +++ 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 6181c3579e..bb9743bb84 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -274,10 +274,12 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ! Langmuir enhancement factor !---- - ice_ocean_boundary%lamult (:,:) = 0._ESMF_KIND_R8 - call state_getimport(importState, 'Sw_lamult', & - isc, iec, jsc, jec, ice_ocean_boundary%lamult, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + if ( associated(ice_ocean_boundary%lamult) ) then + ice_ocean_boundary%lamult (:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Sw_lamult', & + isc, iec, jsc, jec, ice_ocean_boundary%lamult, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + endif !---- ! Partitioned Stokes Drift Components diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index f520462cb1..04786105ed 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -396,6 +396,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call allocate_forcing_type(OS%grid, OS%fluxes, waves=.true.) call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & "If true, enables surface wave modules.", default=.false.) + if (OS%Use_Waves) then + call get_param(param_file, mdl, "WAVE_METHOD", OS%wave_method, default="EMPTY", do_not_log=.true.) + endif ! MOM_wave_interface_init is called regardless of the value of USE_WAVES because ! it also initializes statistical waves. call MOM_wave_interface_init(OS%Time, OS%grid, OS%GV, OS%US, param_file, OS%Waves, OS%diag) From 2d385a061fe83d632f30d606315a2928169890d4 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 5 Aug 2021 10:57:32 -0600 Subject: [PATCH 14/28] Introduce changes needed for CVMix KPP module for the VR12-MA wave coupling option --- .../vertical/MOM_CVMix_KPP.F90 | 128 +++++++++++++----- .../vertical/MOM_diabatic_driver.F90 | 14 +- src/user/MOM_wave_interface.F90 | 40 +++++- 3 files changed, 139 insertions(+), 43 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 083f5ed000..cafdede19b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -15,7 +15,7 @@ module MOM_CVMix_KPP use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type -use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number +use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number, get_wave_method use MOM_domains, only : pass_var use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE @@ -198,6 +198,9 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) # include "version_variable.h" character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module character(len=20) :: string !< local temporary string + character(len=20) :: langmuir_mixing_opt !< langmuir mixing opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_entrainment_opt !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: wave_method logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function) !! False => compute G'(1) as in LMD94 @@ -456,6 +459,18 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) endif call closeParameterBlock(paramFile) + + ! since wave interface control structure is not initialized yet, temporarily read in wave_method from param file + call get_param(paramFile, mdl, "WAVE_METHOD", wave_method, default="EMPTY", do_not_log=.true.) + select case ( trim(wave_method) ) + case ("VR12-MA") + langmuir_mixing_opt = 'LWF16' + langmuir_entrainment_opt = 'LWF16' + case default + langmuir_mixing_opt = 'NONE' + langmuir_entrainment_opt = 'NONE' + end select + call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) call CVMix_init_kpp( Ri_crit=CS%Ri_crit, & @@ -471,6 +486,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) lenhanced_diff=CS%enhance_diffusion,& lnonzero_surf_nonlocal=Cs_is_one ,& lnoDGat1=lnoDGat1 ,& + langmuir_mixing_str=langmuir_mixing_opt,& + langmuir_entrainment_str=langmuir_entrainment_opt,& CVMix_kpp_params_user=CS%KPP_params ) ! Register diagnostics @@ -600,7 +617,7 @@ end function KPP_init !> KPP vertical diffusivity/viscosity and non-local tracer transport subroutine KPP_calculate(CS, G, GV, US, h, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& - nonLocalTransScalar, waves) + nonLocalTransScalar, waves, lamult) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -621,7 +638,8 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & !! [Z2 T-1 ~> m2 s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport [m s-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local trans. [m s-1] - type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement multiplier ! Local variables integer :: i, j, k ! Loop indices @@ -714,23 +732,44 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & Kviscosity(:) = US%Z2_T_to_m2_s * Kv(i,j,:) endif - call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] - iFaceHeight, & ! (in) Height of interfaces [m] - cellHeight, & ! (in) Height of level centers [m] - Kviscosity(:), & ! (in) Original viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] - CS%OBLdepth(i,j), & ! (in) OBL depth [m] - CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent - nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] - nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] - surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] - surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] - GV%ke, & ! (in) Number of levels to compute coeffs for - GV%ke, & ! (in) Number of levels in array shape - CVMix_kpp_params_user=CS%KPP_params ) + if (get_wave_method(Waves) == "VR12-MA") then + call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] + iFaceHeight, & ! (in) Height of interfaces [m] + cellHeight, & ! (in) Height of level centers [m] + Kviscosity(:), & ! (in) Original viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] + CS%OBLdepth(i,j), & ! (in) OBL depth [m] + CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent + nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] + nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] + surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] + GV%ke, & ! (in) Number of levels to compute coeffs for + GV%ke, & ! (in) Number of levels in array shape + Langmuir_EFactor=lamult(i,j),& ! Langmuir enhancement multiplier + CVMix_kpp_params_user=CS%KPP_params ) + else + call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] + iFaceHeight, & ! (in) Height of interfaces [m] + cellHeight, & ! (in) Height of level centers [m] + Kviscosity(:), & ! (in) Original viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] + CS%OBLdepth(i,j), & ! (in) OBL depth [m] + CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent + nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] + nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] + surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] + GV%ke, & ! (in) Number of levels to compute coeffs for + GV%ke, & ! (in) Number of levels in array shape + CVMix_kpp_params_user=CS%KPP_params ) + endif ! safety check, Kviscosity and Kdiffusivity must be >= 0 do k=1, GV%ke+1 @@ -900,7 +939,7 @@ end subroutine KPP_calculate !> Compute OBL depth -subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFlux, Waves) +subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFlux, Waves, lamult) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -914,8 +953,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Velocity j-component [L T-1 ~> m s-1] type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3] type(wave_parameters_CS), optional, pointer :: Waves !< Wave CS + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult!< Langmuir enhancement factor ! Local variables integer :: i, j, k, km1 ! Loop indices @@ -1147,11 +1187,20 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CVMix_kpp_params_user=CS%KPP_params ) !Compute CVMix VT2 - CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & - zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] - N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] - CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + if (get_wave_method(Waves) == "VR12-MA") then + CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & + zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] + EFactor=lamult(i,j), & ! Langmuir enhancement multiplier + CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + else + CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & + zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] + CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + endif !Modify CVMix VT2 IF (CS%LT_VT2_ENHANCEMENT) then @@ -1195,13 +1244,24 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl enddo ! Calculate Bulk Richardson number from eq (21) of LMD94 - BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & - zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] - delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] - delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] - Vt_sqr_cntr=CS%Vt2(i,j,:), & - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] - N_iface=CS%N(i,j,:)) ! Buoyancy frequency [s-1] + if (get_wave_method(Waves) == "VR12-MA") then + BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & + zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] + delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] + delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] + Vt_sqr_cntr=CS%Vt2(i,j,:), & + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1] + EFactor=lamult(i,j)) ! Langmuir enhancement factor + else + BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & + zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] + delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] + delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] + Vt_sqr_cntr=CS%Vt2(i,j,:), & + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] + N_iface=CS%N(i,j,:)) ! Buoyancy frequency [s-1] + endif surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 7a75802a84..b44d4c84b4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1190,11 +1190,19 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) ! The KPP scheme calculates boundary layer diffusivities and non-local transport. - call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + if ( associated(fluxes%lamult) ) then + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + else + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + endif if (associated(Hml)) then call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index b9b239e5df..aee73e8202 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -39,6 +39,7 @@ module MOM_wave_interface ! serious consideration of the full 3d wave-averaged Navier-Stokes ! CL2 effects. public Waves_end ! public interface to deallocate and free wave related memory. +public get_wave_method ! public interface to obtain the wave method string ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -179,6 +180,14 @@ module MOM_wave_interface integer, parameter :: DATAOVR = 1, COUPLER = 2, INPUT = 3 !>@} +! Strings for the wave method +character*(5), parameter :: NULL_STRING = "EMPTY" +character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" +character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" +character*(5), parameter :: DHH85_STRING = "DHH85" +character*(4), parameter :: LF17_STRING = "LF17" +character*(7), parameter :: VR12MA_STRING = "VR12-MA" + contains !> Initializes parameters related to MOM_wave_interface @@ -196,12 +205,6 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) ! This include declares and sets the variable "version". # include "version_variable.h" character*(13) :: TMPSTRING1, TMPSTRING2 - character*(5), parameter :: NULL_STRING = "EMPTY" - character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" - character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" - character*(5), parameter :: DHH85_STRING = "DHH85" - character*(4), parameter :: LF17_STRING = "LF17" - character*(7), parameter :: VR12MA_STRING = "VR12-MA" character*(12), parameter :: DATAOVR_STRING = "DATAOVERRIDE" character*(7), parameter :: COUPLER_STRING = "COUPLER" character*(5), parameter :: INPUT_STRING = "INPUT" @@ -1055,6 +1058,31 @@ subroutine get_Langmuir_Number( LA, G, GV, US, HBL, ustar, i, j, & end subroutine get_Langmuir_Number +!> function to return the wave method string set in the param file +function get_wave_method(CS) + character(:), allocatable :: get_wave_method + type(wave_parameters_CS), pointer :: CS !< Control structure + + if (associated(CS)) then + select case(CS%WaveMethod) + case (Null_WaveMethod) + get_wave_method = NULL_STRING + case (TESTPROF) + get_wave_method = TESTPROF_STRING + case (SURFBANDS) + get_wave_method = SURFBANDS_STRING + case (DHH85) + get_wave_method = DHH85_STRING + case (LF17) + get_wave_method = LF17_STRING + case (VR12MA) + get_wave_method = VR12MA_STRING + end select + else + get_wave_method = NULL_STRING + endif +end function get_wave_method + !> Get SL averaged Stokes drift from Li/FK 17 method !! !! Original description: From 9fe68ac66b03efc5fd92a9fcbde674dd7eb32a53 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 5 Aug 2021 13:27:33 -0600 Subject: [PATCH 15/28] fix doxygen error: add comments for wave method strings --- src/user/MOM_wave_interface.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index aee73e8202..dc2e51248f 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -181,12 +181,12 @@ module MOM_wave_interface !>@} ! Strings for the wave method -character*(5), parameter :: NULL_STRING = "EMPTY" -character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" -character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" -character*(5), parameter :: DHH85_STRING = "DHH85" -character*(4), parameter :: LF17_STRING = "LF17" -character*(7), parameter :: VR12MA_STRING = "VR12-MA" +character*(5), parameter :: NULL_STRING = "EMPTY" !< null wave method string +character*(12), parameter :: TESTPROF_STRING = "TEST_PROFILE" !< test profile string +character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" !< surface bands string +character*(5), parameter :: DHH85_STRING = "DHH85" !< DHH85 wave method string +character*(4), parameter :: LF17_STRING = "LF17" !< LF17 wave method string +character*(7), parameter :: VR12MA_STRING = "VR12-MA" !< VR12-MA wave method string contains From 80e1e2e88ff847fb06ba40123bf29c057bef45f0 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 5 Aug 2021 13:28:15 -0600 Subject: [PATCH 16/28] fix omp error in KPP module: add lamult to shared clauses --- src/parameterizations/vertical/MOM_CVMix_KPP.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index cafdede19b..cdd205fabf 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -679,7 +679,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & !$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, & !$OMP sigmaRatio) & !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, & - !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves) + !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves, lamult) ! loop over horizontal points on processor do j = G%jsc, G%jec do i = G%isc, G%iec @@ -1027,7 +1027,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, & !$OMP BulkRi_1d, zBottomMinusOffset) & !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, & - !$OMP Temp, Salt, waves, tv, GoRho, u, v) + !$OMP Temp, Salt, waves, tv, GoRho, u, v, lamult) do j = G%jsc, G%jec do i = G%isc, G%iec From f8a8e4c4c721830f08ccc83a662fb8215eefa0cf Mon Sep 17 00:00:00 2001 From: jiandewang Date: Mon, 16 Aug 2021 08:53:19 -0400 Subject: [PATCH 17/28] update to gfdl 20210806 (#74) * remove white space and fix comment * Update MOM_oda_incupd.F90 remove unused index bounds, and fix sum_h2 loop. Co-authored-by: pjpegion Co-authored-by: Marshall Ward --- src/ocean_data_assim/MOM_oda_incupd.F90 | 23 +++---------------- .../vertical/MOM_diabatic_driver.F90 | 2 +- 2 files changed, 4 insertions(+), 21 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 12d18587bf..d3199dcb74 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -65,22 +65,6 @@ module MOM_oda_incupd type, public :: oda_incupd_CS ; private integer :: nz !< The total number of layers. integer :: nz_data !< The total number of arbritary layers (used by older code). - integer :: isc !< The starting i-index of the computational domain at h. - integer :: iec !< The ending i-index of the computational domain at h. - integer :: jsc !< The starting j-index of the computational domain at h. - integer :: jec !< The ending j-index of the computational domain at h. - integer :: IscB !< The starting I-index of the computational domain at u/v. - integer :: IecB !< The ending I-index of the computational domain at u/v. - integer :: JscB !< The starting J-index of the computational domain at u/v. - integer :: JecB !< The ending J-index of the computational domain at h. - integer :: isd !< The starting i-index of the data domain at h. - integer :: ied !< The ending i-index of the data domain at h. - integer :: jsd !< The starting j-index of the data domain at h. - integer :: jed !< The ending j-index of the data domain at h. - integer :: IsdB !< The starting I-index of the data domain at u/v. - integer :: IedB !< The ending I-index of the data domain at u/v. - integer :: JsdB !< The starting J-index of the data domain at u/v. - integer :: JedB !< The ending J-index of the data domain at h. integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_oda_incupd_field @@ -224,9 +208,6 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res default=.true.) CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%isdB = G%isdB ; CS%iedB = G%iedB; CS%jsdB = G%jsdB ; CS%jedB = G%jedB ! increments on horizontal grid if (.not. CS%incupdDataOngrid) call MOM_error(FATAL,'initialize_oda_incupd: '// & @@ -421,6 +402,8 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) do k=1,nz_data sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) enddo ! get temperature @@ -450,7 +433,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo endif enddo; enddo - + ! remap u to h_obs to get increment if (CS%uv_inc) then call pass_var(h, G%Domain) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 86a83c58a9..fb759b13b2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -111,7 +111,7 @@ module MOM_diabatic_driver !! domain. The exact location and properties of !! those sponges are set by calls to !! initialize_sponge and set_up_sponge_field. - logical :: use_oda_incupd !!< If True, DA incremental update is + logical :: use_oda_incupd !< If True, DA incremental update is !! applied everywhere logical :: use_geothermal !< If true, apply geothermal heating. logical :: use_int_tides !< If true, use the code that advances a separate set From 734e537826cc560fb0e0447067519ff687dbb7a1 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 17 Aug 2021 07:23:34 -0600 Subject: [PATCH 18/28] change VR12-MA wave_method string to EFACTOR --- config_src/drivers/nuopc_cap/mom_cap.F90 | 4 +-- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 2 +- src/user/MOM_wave_interface.F90 | 25 ++++++++++--------- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 86219c56a4..67d650aded 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -708,7 +708,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call query_ocean_state(ocean_state, use_waves=use_waves, wave_method=wave_method) if (use_waves) then call query_ocean_state(ocean_state, NumWaveBands=Ice_ocean_boundary%num_stk_bands) - if (wave_method == "VR12-MA") then + if (wave_method == "EFACTOR") then allocate( Ice_ocean_boundary%lamult(isc:iec,jsc:jec) ) Ice_ocean_boundary%lamult = 0.0 else @@ -761,7 +761,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_runoff_heat_flx" , "will provide") !call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_calving_heat_flx" , "will provide") if (use_waves) then - if (wave_method == "VR12-MA") then + if (wave_method == "EFACTOR") then call fld_list_add(fldsToOcn_num, fldsToOcn, "Sw_lamult" , "will provide") else if (Ice_ocean_boundary%num_stk_bands > 3) then diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 04786105ed..39e1046c1c 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -585,7 +585,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid, OS%US) if (OS%use_waves) then - if (OS%wave_method /= "VR12-MA") then + if (OS%wave_method /= "EFACTOR") then call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) endif endif diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index dc2e51248f..2f0d95a62d 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -174,7 +174,7 @@ module MOM_wave_interface ! Switches needed in import_stokes_drift !>@{ Enumeration values for the wave method -integer, parameter :: TESTPROF = 0, SURFBANDS = 1, DHH85 = 2, LF17 = 3, VR12MA = 4, NULL_WaveMethod = -99 +integer, parameter :: TESTPROF = 0, SURFBANDS = 1, DHH85 = 2, LF17 = 3, EFACTOR = 4, NULL_WaveMethod = -99 !>@} !>@{ Enumeration values for the wave data source integer, parameter :: DATAOVR = 1, COUPLER = 2, INPUT = 3 @@ -186,7 +186,7 @@ module MOM_wave_interface character*(13), parameter :: SURFBANDS_STRING = "SURFACE_BANDS" !< surface bands string character*(5), parameter :: DHH85_STRING = "DHH85" !< DHH85 wave method string character*(4), parameter :: LF17_STRING = "LF17" !< LF17 wave method string -character*(7), parameter :: VR12MA_STRING = "VR12-MA" !< VR12-MA wave method string +character*(7), parameter :: EFACTOR_STRING = "EFACTOR" !< EFACTOR (based on vr12-ma) wave method string contains @@ -288,10 +288,11 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) " wave spectrum with prescribed values. \n"// & " LF17 - Infers Stokes drift profile from wind \n"// & " speed following Li and Fox-Kemper 2017.\n"// & - " VR12-MA - Applies an enhancement factor to the KPP\n"//& - " turbulent velocity scale received from \n"// & - " WW3 and is based on the surface layer \n"// & - " and projected Langmuir number. (Li 2016)\n", & + " EFACTOR - Applies an enhancement factor to the KPP\n"//& + " turbulent velocity scale received \n"// & + " directly from WW3 and is based on the \n"// & + " surface layer and projected Langmuir \n"// & + " number (Li 2016)\n", & units='', default=NULL_STRING) select case (TRIM(TMPSTRING1)) case (NULL_STRING)! No Waves @@ -390,8 +391,8 @@ subroutine MOM_wave_interface_init(time, G, GV, US, param_file, CS, diag ) default=.false.) case (LF17_STRING)!Li and Fox-Kemper 17 wind-sea Langmuir number CS%WaveMethod = LF17 - case (VR12MA_STRING)!Li and Fox-Kemper 16 - CS%WaveMethod = VR12MA + case (EFACTOR_STRING)!Li and Fox-Kemper 16 + CS%WaveMethod = EFACTOR case default call MOM_error(FATAL,'Check WAVE_METHOD.') end select @@ -770,8 +771,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo CS%DHH85_is_set = .true. endif - elseif (CS%WaveMethod==VR12MA) then - return ! todo + elseif (CS%WaveMethod==EFACTOR) then + return ! pass else! Keep this else, fallback to 0 Stokes drift do kk= 1,GV%ke do jj = G%jsd,G%jed @@ -1075,8 +1076,8 @@ function get_wave_method(CS) get_wave_method = DHH85_STRING case (LF17) get_wave_method = LF17_STRING - case (VR12MA) - get_wave_method = VR12MA_STRING + case (EFACTOR) + get_wave_method = EFACTOR_STRING end select else get_wave_method = NULL_STRING From d734536c0f4b3e5213bb9088694725421b940107 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 17 Aug 2021 08:04:41 -0600 Subject: [PATCH 19/28] Refactor the way Langmuir entrainment and enhancement factor are computed and applied: (1) When available, pass the enhancement factor received from WW3. Otherwise, let CVmix compute the enhancement factor. (2) Instead of explicitly multiplying diffusivities and viscosities with the enhancement factor, let CVMix handle enhancement internally. --- .../vertical/MOM_CVMix_KPP.F90 | 251 +++++++----------- 1 file changed, 98 insertions(+), 153 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index cdd205fabf..0abc5d1b2a 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -198,8 +198,8 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) # include "version_variable.h" character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module character(len=20) :: string !< local temporary string - character(len=20) :: langmuir_mixing_opt !< langmuir mixing opt to be passed to CVMix, e.g., LWF16 - character(len=20) :: langmuir_entrainment_opt !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_mixing_opt = 'NONE' !< langmuir mixing opt to be passed to CVMix, e.g., LWF16 + character(len=20) :: langmuir_entrainment_opt = 'NONE' !< langmuir entrainment opt to be passed to CVMix, e.g., LWF16 character(len=20) :: wave_method logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function) @@ -418,10 +418,16 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) '\t RW16 = Function of Langmuir number based on RW16', & default='CONSTANT') select case ( trim(string)) - case ("CONSTANT") ; CS%LT_K_METHOD = LT_K_MODE_CONSTANT - case ("VR12") ; CS%LT_K_METHOD = LT_K_MODE_VR12 - case ("RW16") ; CS%LT_K_METHOD = LT_K_MODE_RW16 - case default ; call MOM_error(FATAL,"KPP_init: "//& + case ("CONSTANT") + CS%LT_K_METHOD = LT_K_MODE_CONSTANT + case ("VR12") + CS%LT_K_METHOD = LT_K_MODE_VR12 + langmuir_mixing_opt = 'LWF16' + case ("RW16") + CS%LT_K_METHOD = LT_K_MODE_RW16 + langmuir_mixing_opt = 'RWHGK16' + case default + call MOM_error(FATAL,"KPP_init: "//& "Unrecognized KPP_LT_K_METHOD option: "//trim(string)) end select if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then @@ -444,12 +450,20 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) '\t LF17 = Function of Langmuir number based on LF17', & default='CONSTANT') select case ( trim(string)) - case ("CONSTANT") ; CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT - case ("VR12") ; CS%LT_VT2_METHOD = LT_VT2_MODE_VR12 - case ("RW16") ; CS%LT_VT2_METHOD = LT_VT2_MODE_RW16 - case ("LF17") ; CS%LT_VT2_METHOD = LT_VT2_MODE_LF17 - case default ; call MOM_error(FATAL,"KPP_init: "//& - "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string)) + case ("CONSTANT") + CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT + case ("VR12") + CS%LT_VT2_METHOD = LT_VT2_MODE_VR12 + langmuir_entrainment_opt = 'LWF16' + case ("RW16") + CS%LT_VT2_METHOD = LT_VT2_MODE_RW16 + langmuir_entrainment_opt = 'RWHGK16' + case ("LF17") + CS%LT_VT2_METHOD = LT_VT2_MODE_LF17 + langmuir_entrainment_opt = 'LF17' + case default + call MOM_error(FATAL,"KPP_init: "//& + "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string)) end select if (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then call get_param(paramFile, mdl, "KPP_VT2_ENH_FAC",CS%KPP_VT2_ENH_FAC , & @@ -460,17 +474,6 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) call closeParameterBlock(paramFile) - ! since wave interface control structure is not initialized yet, temporarily read in wave_method from param file - call get_param(paramFile, mdl, "WAVE_METHOD", wave_method, default="EMPTY", do_not_log=.true.) - select case ( trim(wave_method) ) - case ("VR12-MA") - langmuir_mixing_opt = 'LWF16' - langmuir_entrainment_opt = 'LWF16' - case default - langmuir_mixing_opt = 'NONE' - langmuir_entrainment_opt = 'NONE' - end select - call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) call CVMix_init_kpp( Ri_crit=CS%Ri_crit, & @@ -732,62 +735,16 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & Kviscosity(:) = US%Z2_T_to_m2_s * Kv(i,j,:) endif - if (get_wave_method(Waves) == "VR12-MA") then - call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] - iFaceHeight, & ! (in) Height of interfaces [m] - cellHeight, & ! (in) Height of level centers [m] - Kviscosity(:), & ! (in) Original viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] - CS%OBLdepth(i,j), & ! (in) OBL depth [m] - CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent - nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] - nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] - surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] - surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] - GV%ke, & ! (in) Number of levels to compute coeffs for - GV%ke, & ! (in) Number of levels in array shape - Langmuir_EFactor=lamult(i,j),& ! Langmuir enhancement multiplier - CVMix_kpp_params_user=CS%KPP_params ) - else - call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] - iFaceHeight, & ! (in) Height of interfaces [m] - cellHeight, & ! (in) Height of level centers [m] - Kviscosity(:), & ! (in) Original viscosity [m2 s-1] - Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] - Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] - CS%OBLdepth(i,j), & ! (in) OBL depth [m] - CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent - nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] - nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] - surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] - surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] - GV%ke, & ! (in) Number of levels to compute coeffs for - GV%ke, & ! (in) Number of levels in array shape - CVMix_kpp_params_user=CS%KPP_params ) - endif - - ! safety check, Kviscosity and Kdiffusivity must be >= 0 - do k=1, GV%ke+1 - if (Kviscosity(k) < 0. .or. Kdiffusivity(k,1) < 0.) then - call MOM_error(FATAL,"KPP_calculate, after CVMix_coeffs_kpp: "// & - "Negative vertical viscosity or diffusivity has been detected. " // & - "This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //& - "You might consider using the default options for these parameters." ) - endif - enddo - IF (CS%LT_K_ENHANCEMENT) then if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then LangEnhK = CS%KPP_K_ENH_FAC elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then - ! Added minimum value for La_SL, so removed maximum value for LangEnhK. - LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & + if (present(lamult)) then + LangEnhK = lamult(i,j) + else + LangEnhK = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & (5.4*CS%La_SL(i,j))**(-4)) + endif elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then !This maximum value is proposed in Reichl et al., 2016 JPO formula LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j)) @@ -796,26 +753,59 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT") LangEnhK = 1.0 endif + + ! diffusivities don't need to be enhanced below anymore since LangEnhK is applied within CVMix. + ! todo: need to double check if the distinction between the two different options of LT_K_SHAPE may need to be + ! treated specially. do k=1,GV%ke if (CS%LT_K_SHAPE== LT_K_CONSTANT) then if (CS%id_EnhK > 0) CS%EnhK(i,j,:) = LangEnhK - Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK - Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK - Kviscosity(k) = Kviscosity(k) * LangEnhK + !Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK + !Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK + !Kviscosity(k) = Kviscosity(k) * LangEnhK elseif (CS%LT_K_SHAPE == LT_K_SCALED) then sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) SigmaRatio = sigma * (1. - sigma)**2 / 0.148148037 if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio) - Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) - Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) - Kviscosity(k) = Kviscosity(k) * ( 1. + & - ( LangEnhK - 1.)*sigmaRatio) + !Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) + !Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) + !Kviscosity(k) = Kviscosity(k) * ( 1. + & + ! ( LangEnhK - 1.)*sigmaRatio) endif enddo endif + call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1] + iFaceHeight, & ! (in) Height of interfaces [m] + cellHeight, & ! (in) Height of level centers [m] + Kviscosity(:), & ! (in) Original viscosity [m2 s-1] + Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1] + Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1] + CS%OBLdepth(i,j), & ! (in) OBL depth [m] + CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent + nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim] + nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim] + surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1] + surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3] + GV%ke, & ! (in) Number of levels to compute coeffs for + GV%ke, & ! (in) Number of levels in array shape + Langmuir_EFactor=LangEnhK,& ! Langmuir enhancement multiplier + CVMix_kpp_params_user=CS%KPP_params ) + + ! safety check, Kviscosity and Kdiffusivity must be >= 0 + do k=1, GV%ke+1 + if (Kviscosity(k) < 0. .or. Kdiffusivity(k,1) < 0.) then + call MOM_error(FATAL,"KPP_calculate, after CVMix_coeffs_kpp: "// & + "Negative vertical viscosity or diffusivity has been detected. " // & + "This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //& + "You might consider using the default options for these parameters." ) + endif + enddo + ! Over-write CVMix NLT shape function with one of the following choices. ! The CVMix code has yet to update for thse options, so we compute in MOM6. ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with @@ -993,8 +983,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl integer :: kk, ksfc, ktmp ! For Langmuir Calculations - real :: LangEnhW ! Langmuir enhancement for turbulent velocity scale - real, dimension(GV%ke) :: LangEnhVt2 ! Langmuir enhancement for unresolved shear + real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear real, dimension(GV%ke) :: U_H, V_H real :: MLD_GUESS, LA real :: surfHuS, surfHvS, surfUs, surfVs, wavedir, currentdir @@ -1186,87 +1175,43 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1] CVMix_kpp_params_user=CS%KPP_params ) - !Compute CVMix VT2 - if (get_wave_method(Waves) == "VR12-MA") then - CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & - zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] - N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] - EFactor=lamult(i,j), & ! Langmuir enhancement multiplier - CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - else - CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & - zt_cntr=cellHeight(1:GV%ke), & ! Depth of cell center [m] - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] - N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] - CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - endif - - !Modify CVMix VT2 + ! Determine the enhancement factor for unresolved shear IF (CS%LT_VT2_ENHANCEMENT) then IF (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then - do k=1,GV%ke - LangEnhVT2(k) = CS%KPP_VT2_ENH_FAC - enddo + LangEnhVT2 = CS%KPP_VT2_ENH_FAC elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed. - enhvt2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & - (5.4*CS%La_SL(i,j))**(-4)) - do k=1,GV%ke - LangEnhVT2(k) = enhvt2 - enddo - elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_RW16) then - !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed. - enhvt2 = 1. + 2.3*CS%La_SL(i,j)**(-0.5) - do k=1,GV%ke - LangEnhVT2(k) = enhvt2 - enddo - elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_LF17) then - CS%CS=cvmix_get_kpp_real('c_s',CS%KPP_params) - do k=1,GV%ke - WST = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellHeight(k)))**(1./3.) - LangEnhVT2(k) = sqrt((0.15*WST**3. + 0.17*surfFricVel**3.* & - (1.+0.49*CS%La_SL(i,j)**(-2.))) / & - (0.2*ws_1d(k)**3/(CS%cs*CS%surf_layer_ext*CS%vonKarman**4.))) - enddo + if (present(lamult)) then + LangEnhVT2 = lamult(i,j) + else + LangEnhVT2 = sqrt(1.+(1.5*CS%La_SL(i,j))**(-2) + & + (5.4*CS%La_SL(i,j))**(-4)) + endif else - !This shouldn't be reached. - !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2") - LangEnhVT2(:) = 1.0 + ! for other methods (e.g., LT_VT2_MODE_RW16, LT_VT2_MODE_LF17), the enhancement factor is + ! computed internally within CVMix using LaSL, bfsfc, and ustar to be passed to CVMix. + LangEnhVT2 = 1.0 endif else - LangEnhVT2(:) = 1.0 + LangEnhVT2 = 1.0 endif - do k=1,GV%ke - CS%Vt2(i,j,k)=CS%Vt2(i,j,k)*LangEnhVT2(k) - if (CS%id_EnhVt2 > 0) CS%EnhVt2(i,j,k)=LangEnhVT2(k) - enddo + surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! Calculate Bulk Richardson number from eq (21) of LMD94 - if (get_wave_method(Waves) == "VR12-MA") then - BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & - zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] - delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] - delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] - Vt_sqr_cntr=CS%Vt2(i,j,:), & - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] - N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1] - EFactor=lamult(i,j)) ! Langmuir enhancement factor - else - BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & - zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] - delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] - delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] - Vt_sqr_cntr=CS%Vt2(i,j,:), & - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] - N_iface=CS%N(i,j,:)) ! Buoyancy frequency [s-1] - endif + BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( & + zt_cntr = cellHeight(1:GV%ke), & ! Depth of cell center [m] + delta_buoy_cntr=GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1] + delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2] + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency [s-1] + EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] + LaSL = CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] + bfsfc = surfBuoyFlux, & ! surface buoyancy flux [m2 s-3] + uStar = uStar(i,j), & ! surface friction velocity [m s-1] + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters - surfBuoyFlux = buoy_scale * buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit - ! h to Monin-Obukov (default is false, ie. not used) - call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number iFaceHeight, & ! (in) Height of interfaces [m] From 9011801b617ef6c09a2e48f8685532886122d10f Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Mon, 23 Aug 2021 09:31:03 -0600 Subject: [PATCH 20/28] Various dimension rescaling fixes correct scale argument in get_param and chksum calls add dimension scaling terms to some expressions add some chksum calls --- .../mct_cap/mom_surface_forcing_mct.F90 | 3 ++- config_src/drivers/nuopc_cap/mom_cap.F90 | 2 +- .../solo_driver/MESO_surface_forcing.F90 | 3 ++- src/core/MOM_forcing_type.F90 | 2 +- src/parameterizations/lateral/MOM_MEKE.F90 | 18 +++++++++++++----- .../lateral/MOM_lateral_mixing_coeffs.F90 | 11 ++++++----- .../lateral/MOM_thickness_diffuse.F90 | 4 ++-- .../vertical/MOM_CVMix_KPP.F90 | 6 +++--- .../vertical/MOM_CVMix_conv.F90 | 4 ++-- 9 files changed, 32 insertions(+), 21 deletions(-) diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index e2a557daff..7b32270b4c 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -1245,7 +1245,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "If true, use a 2-dimensional gustiness supplied from "//& "an input file", default=.false.) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & - "The background gustiness in the winds.", units="Pa", default=0.0) + "The background gustiness in the winds.", units="Pa", default=0.0, & + scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) if (CS%read_gust_2d) then call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & "The file in which the wind gustiness is found in "//& diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 0090468a05..ef3d9867b9 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -1113,7 +1113,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) k = k + 1 ! Increment position within gindex if (mask(k) /= 0) then mesh_areas(k) = dataPtr_mesh_areas(k) - model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 + model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 mod2med_areacor(k) = model_areas(k) / mesh_areas(k) med2mod_areacor(k) = mesh_areas(k) / model_areas(k) end if diff --git a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 index 679f147797..7c778d9753 100644 --- a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 @@ -243,7 +243,8 @@ subroutine MESO_surface_forcing_init(Time, G, US, param_file, diag, CS) "parameters from vertical units of m to kg m-2.", & units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & - "The background gustiness in the winds.", units="Pa", default=0.0) + "The background gustiness in the winds.", units="Pa", default=0.0, & + scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7d8375a8af..20e9d6e698 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -1117,7 +1117,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%heat_content_lrunoff)) & call hchksum(fluxes%heat_content_lrunoff, mesg//" fluxes%heat_content_lrunoff", G%HI, & - haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_frunoff)) & call hchksum(fluxes%heat_content_frunoff, mesg//" fluxes%heat_content_frunoff", G%HI, & haloshift=hshift, scale=US%QRZ_T_to_W_m2) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 762b2edaea..25946eb631 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -81,7 +81,7 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] - real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [T-1 ~> s-1]. logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. @@ -340,6 +340,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif + if (CS%debug) then + call hchksum(src, "MEKE src", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) + endif + ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -534,6 +538,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif endif ! MEKE_KH>=0 + if (CS%debug) then + call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, scale=US%L_T_to_m_s**2) + endif + ! do j=js,je ; do i=is,ie ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) ! enddo ; enddo @@ -688,7 +696,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 else FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points @@ -824,7 +832,7 @@ subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 enddo ; enddo if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) @@ -957,7 +965,7 @@ subroutine MEKE_lengthScales_0d(CS, US, area, beta, depth, Rd_dx, SN, EKE, & ! Z Leady = 0. endif if (CS%use_min_lscale) then - LmixScale = 1.e7 + LmixScale = 1.e7*US%m_to_L if (CS%aDeform*Ldeform > 0.) LmixScale = min(LmixScale,CS%aDeform*Ldeform) if (CS%aFrict *Lfrict > 0.) LmixScale = min(LmixScale,CS%aFrict *Lfrict) if (CS%aRhines*Lrhines > 0.) LmixScale = min(LmixScale,CS%aRhines*Lrhines) @@ -1069,7 +1077,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) if (CS%MEKE_equilibrium_restoring) then call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & - default=1e6, scale=US%T_to_s) + default=1e6, scale=US%s_to_T) CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index fb70f5d679..dabeb27159 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -506,10 +506,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency - !! at u-points [T-2 ~> s-2] + !! at u-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency - !! at v-points [T-2 ~> s-2] + !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), pointer :: CS !< Variable mixing coefficients type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -654,9 +654,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C endif if (CS%debug) then - call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) + call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, & + scale=US%Z_to_L, haloshift=1) call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, & - scale=US%s_to_T**2, scalar_pair=.true.) + scale=US%L_to_Z**2 * US%s_to_T**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & scale=US%s_to_T, scalar_pair=.true.) endif @@ -1484,7 +1485,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB)) ; CS%Depth_fn_v(:,:) = 0.0 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, & "The depth above which KHTH is scaled away.",& - units="m", default=1000.) + units="m", scale=US%m_to_Z, default=1000.) call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, & "The exponent used in the depth dependent scaling function for KHTH.",& units="nondim", default=3.0) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index da62ffc6b7..4d038ed31e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -433,7 +433,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m) if (use_stored_slopes) then call uvchksum("VarMix%slope_[xy]", VarMix%slope_x, VarMix%slope_y, & - G%HI, haloshift=0) + G%HI, haloshift=0, scale=US%Z_to_L) endif if (associated(tv%eqn_of_state)) then call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1) @@ -505,7 +505,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(GV%m_to_H,GV%Z_to_H*G%bathyT(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 083f5ed000..dec8fad571 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -644,7 +644,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%debug) then call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m) call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T) - call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0) + call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) endif @@ -967,8 +967,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%debug) then call hchksum(Salt, "KPP in: S",G%HI,haloshift=0) call hchksum(Temp, "KPP in: T",G%HI,haloshift=0) - call hchksum(u, "KPP in: u",G%HI,haloshift=0) - call hchksum(v, "KPP in: v",G%HI,haloshift=0) + call hchksum(u, "KPP in: u",G%HI,haloshift=0,scale=US%L_T_to_m_s) + call hchksum(v, "KPP in: v",G%HI,haloshift=0,scale=US%L_T_to_m_s) endif call cpu_clock_begin(id_clock_KPP_compute_BLD) diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 8b007a0b11..a615c9f40b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -288,8 +288,8 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) ! if (CS%id_kv_conv > 0) & ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kd)) call hchksum(Kv, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) + if (present(Kd)) call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) endif ! send diagnostics to post_data From e87adb88a710323c3d3e8a069f2f2635450a6aa7 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Mon, 23 Aug 2021 10:27:40 -0600 Subject: [PATCH 21/28] dimension rescaling fixes for MOM_CFC_cap tracers --- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 2 +- src/core/MOM_forcing_type.F90 | 14 +++-- src/tracer/MOM_CFC_cap.F90 | 52 ++++++++++--------- 3 files changed, 37 insertions(+), 31 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index c83e32711c..4eb31d2434 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -576,7 +576,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! CFCs if (CS%use_CFC) then - call CFC_cap_fluxes(fluxes, sfc_state, G, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) endif if (associated(IOB%salt_flux)) then diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 20e9d6e698..77eda959ad 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -185,8 +185,8 @@ module MOM_forcing_type ! CFC-related arrays needed in the MOM_CFC_cap module real, pointer, dimension(:,:) :: & - cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [mol m-2 s-1]. - cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [mol m-2 s-1]. + cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> mol m-2 s-1]. + cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> mol m-2 s-1]. ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] @@ -1101,9 +1101,13 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) if (associated(fluxes%p_surf)) & call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift , scale=US%RL2_T2_to_Pa) if (associated(fluxes%u10_sqr)) & - call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%Z_to_m**2*US%s_to_T**2) + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%L_to_m**2*US%s_to_T**2) if (associated(fluxes%ice_fraction)) & call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) + if (associated(fluxes%cfc11_flux)) & + call hchksum(fluxes%cfc11_flux, mesg//" fluxes%cfc11_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + if (associated(fluxes%cfc12_flux)) & + call hchksum(fluxes%cfc12_flux, mesg//" fluxes%cfc12_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & @@ -1314,14 +1318,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, if (use_cfcs) then handles%id_cfc11 = register_diag_field('ocean_model', 'cfc11_flux', diag%axesT1, Time, & 'Gas exchange flux of CFC11 into the ocean ', 'mol m-2 s-1', & - conversion= US%m_to_Z**2*US%s_to_T,& + conversion= US%Z_to_m*US%s_to_T,& cmor_field_name='fgcfc11', & cmor_long_name='Surface Downward CFC11 Flux', & cmor_standard_name='surface_downward_cfc11_flux') handles%id_cfc12 = register_diag_field('ocean_model', 'cfc12_flux', diag%axesT1, Time, & 'Gas exchange flux of CFC12 into the ocean ', 'mol m-2 s-1', & - conversion= US%m_to_Z**2*US%s_to_T,& + conversion= US%Z_to_m*US%s_to_T,& cmor_field_name='fgcfc12', & cmor_long_name='Surface Downward CFC12 Flux', & cmor_standard_name='surface_downward_cfc12_flux') diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 58ba419c62..583c906a0b 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -302,20 +302,18 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C ! Local variables real, pointer, dimension(:,:,:) :: CFC11 => NULL(), CFC12 => NULL() real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] - real :: scale_factor ! convert from [Conc. m s-1] to [Conc. kg m-2 T-1] + real :: scale_factor ! convert from cfc1[12]_flux to units of sfc_flux in tracer_vertdiff integer :: i, j, k, m, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - ! CFC_flux **unscaled** units are [mol m-2 s-1] which is the same as [CU kg m-2 s-1], - ! where CU = mol kg-1. However, CFC_flux has been scaled already. - ! CFC_flux **scaled** units are [CU L T-1 R]. We want [CU kg m-2 T-1] because - ! these units are what tracer_vertdiff needs. - ! Therefore, we need to convert [L R] to [kg m-2] using scale_factor. - scale_factor = US%L_to_Z * US%RZ_to_kg_m2 ! this will give [CU kg m-2 T-1], which is what verdiff wants - if (.not.associated(CS)) return + ! set factor to convert from cfc1[12]_flux units to tracer_vertdiff argument sfc_flux units + ! cfc1[12]_flux units are CU Z T-1 kg m-3 + ! tracer_vertdiff argument sfc_flux units are CU kg m-2 T-1 + scale_factor = US%Z_to_m + CFC11 => CS%CFC11 ; CFC12 => CS%CFC12 ! Use a tridiagonal solver to determine the concentrations after the @@ -340,8 +338,8 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C endif ! If needed, write out any desired diagnostics from tracer sources & sinks here. - if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, CFC11*GV%Rho0, CS%diag) - if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, CFC12*GV%Rho0, CS%diag) + if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC11, CS%diag) + if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC12, CS%diag) end subroutine CFC_cap_column_physics @@ -421,14 +419,15 @@ end subroutine CFC_cap_surface_state !> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM !! concentration, and calculating the solubility, Schmidt number, and gas exchange. -subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc12_atm) +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type type(surface), intent(in ) :: sfc_state !< A structure containing fields !! that describe the surface state of the ocean. type(forcing), intent(inout) :: fluxes !< A structure containing pointers !! to thermodynamic and tracer forcing fields. Unused fields !! have NULL ptrs. - real, intent(in ) :: Rho0 !< The mean ocean density [kg m-3] + real, intent(in ) :: Rho0 !< The mean ocean density [R ~> kg m-3] type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the !! CFC's concentration in the atmosphere. integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. @@ -440,7 +439,7 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc CFC12_Csurf, & ! The CFC-12 surface concentrations times the Schmidt number term [mol kg-1]. CFC11_alpha, & ! The CFC-11 solubility [mol kg-1 atm-1]. CFC12_alpha, & ! The CFC-12 solubility [mol kg-1 atm-1]. - kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [m s-1]. + kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1]. cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) ! [mol kg-1]. cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] @@ -451,7 +450,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12. real :: sc_no_term ! A term related to the Schmidt number. - real :: kw_coeff ! A coefficient used to scale the piston velocity [L T-1 ~> m s-1] + real :: kw_coeff ! A coefficient used to scale the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + real :: Rho0_kg_m3 ! Rho0 in non-scaled units [kg m-3] real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. integer :: i, j, m, is, ie, js, je @@ -479,6 +479,15 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc "has not been implemented yet.") endif + !--------------------------------------------------------------------- + ! Gas exchange/piston velocity parameter + !--------------------------------------------------------------------- + ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 + ! = 6.97e-7 m/s s^2/m^2 [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + kw_coeff = (US%m_to_Z*US%s_to_T*US%L_to_m**2) * 6.97e-7 + + Rho0_kg_m3 = Rho0 * US%R_to_kg_m3 + do j=js,je ; do i=is,ie ! ta in hectoKelvin ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) @@ -498,21 +507,14 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc CFC12_alpha(i,j) = alpha_12 * sc_no_term CFC12_Csurf(i,j) = sfc_state%sfc_CFC12(i,j) * sc_no_term - !--------------------------------------------------------------------- - ! Gas exchange/piston velocity parameter - !--------------------------------------------------------------------- - ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 - ! 6.97e-07 is used to convert from cm/hr to m/s, kw = m/s [L/T] - kw_coeff = 6.97e-07 * G%US%m_to_L * G%US%T_to_s - kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) ! air concentrations and cfcs BC's fluxes ! CFC flux units: mol kg-1 s-1 kg m-2 - cair(i,j) = pa_to_atm * CFC11_alpha(i,j) * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc11_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC11_Csurf(i,j)) * Rho0 - cair(i,j) = pa_to_atm * CFC12_alpha(i,j) * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc12_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC12_Csurf(i,j)) * Rho0 + cair(i,j) = US%R_to_kg_m3*US%L_T_to_m_s**2*pa_to_atm * CFC11_alpha(i,j) * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC11_Csurf(i,j)) * Rho0_kg_m3 + cair(i,j) = US%R_to_kg_m3*US%L_T_to_m_s**2*pa_to_atm * CFC12_alpha(i,j) * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC12_Csurf(i,j)) * Rho0_kg_m3 enddo ; enddo end subroutine CFC_cap_fluxes From bc198197c4cd0f4e78bf0b6194112902be02c35e Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Mon, 23 Aug 2021 11:06:18 -0600 Subject: [PATCH 22/28] let register_tracer handle flux_units for MOM_CFC_cap tracers --- src/tracer/MOM_CFC_cap.F90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 583c906a0b..ffa05d9217 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -87,7 +87,6 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers. - character(len=48) :: flux_units ! The units for tracer fluxes. character(len=48) :: dummy ! Dummy variable to store params that need to be logged here. logical :: register_CFC_cap integer :: isd, ied, jsd, jed, nz, m @@ -142,8 +141,6 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) CS%CFC11_name = "CFC_11" ; CS%CFC12_name = "CFC_12" CS%CFC11_desc = var_desc(CS%CFC11_name,"mol kg-1","Moles Per Unit Mass of CFC-11 in sea water", caller=mdl) CS%CFC12_desc = var_desc(CS%CFC12_name,"mol kg-1","Moles Per Unit Mass of CFC-12 in sea water", caller=mdl) - if (GV%Boussinesq) then ; flux_units = "mol s-1" - else ; flux_units = "mol m-3 kg s-1" ; endif allocate(CS%CFC11(isd:ied,jsd:jed,nz)) ; CS%CFC11(:,:,:) = 0.0 allocate(CS%CFC12(isd:ied,jsd:jed,nz)) ; CS%CFC12(:,:,:) = 0.0 @@ -154,13 +151,11 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! Register CFC11 for horizontal advection, diffusion, and restarts. call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & tr_desc=CS%CFC11_desc, registry_diags=.true., & - flux_units=flux_units, & restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) ! Do the same for CFC12 tr_ptr => CS%CFC12 call register_tracer(tr_ptr, Tr_Reg, param_file, HI, GV, & tr_desc=CS%CFC12_desc, registry_diags=.true., & - flux_units=flux_units, & restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) CS%tr_Reg => tr_Reg From fa96ae332cd90ec693751d7418e2e11b3949d596 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Mon, 23 Aug 2021 11:09:24 -0600 Subject: [PATCH 23/28] correct prepending inputdir to CFC_BC_file --- .../drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 | 4 +++- src/tracer/MOM_CFC_cap.F90 | 8 +++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 4eb31d2434..5d1070954b 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -1396,7 +1396,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "internal BC generation (TODO).", default=" ", do_not_log=.true.) if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then ! Add the directory if CFC_BC_file is not already a complete path. - CS%CFC_BC_file = trim(slasher(CS%inputdir))//trim(CS%CFC_BC_file) + CS%CFC_BC_file = trim(CS%inputdir) // trim(CS%CFC_BC_file) + endif + if (len_trim(CS%CFC_BC_file) > 0) then call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, & "The name of the variable representing CFC-11 in "//& "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.) diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index ffa05d9217..2a7516b4d2 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -87,7 +87,7 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers. - character(len=48) :: dummy ! Dummy variable to store params that need to be logged here. + character(len=200) :: dummy ! Dummy variable to store params that need to be logged here. logical :: register_CFC_cap integer :: isd, ied, jsd, jed, nz, m @@ -128,6 +128,12 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) "found (units must be parts per trillion), or an empty string for "//& "internal BC generation (TODO).", default=" ") if ((len_trim(dummy) > 0) .and. (scan(dummy,'/') == 0)) then + ! Add the directory if dummy is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + dummy = trim(slasher(inputdir))//trim(dummy) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", dummy) + endif + if (len_trim(dummy) > 0) then call get_param(param_file, mdl, "CFC11_VARIABLE", dummy, & "The name of the variable representing CFC-11 in "//& "CFC_BC_FILE.", default="CFC_11") From f6524adeeb5e4f0a2738ae4b104c7941aa0b5453 Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Tue, 24 Aug 2021 06:56:14 -0600 Subject: [PATCH 24/28] Add run-time parameters CFC11_IC_VAL and CFC12_IC_VAL. --- src/tracer/MOM_CFC_cap.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 2a7516b4d2..959348b803 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -120,6 +120,12 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) "if they are not found in the restart files. Otherwise "//& "it is a fatal error if tracers are not found in the "//& "restart files of a restarted run.", default=.false.) + call get_param(param_file, mdl, "CFC11_IC_VAL", CS%CFC11_IC_val, & + "Value that CFC_11 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) + call get_param(param_file, mdl, "CFC12_IC_VAL", CS%CFC12_IC_val, & + "Value that CFC_12 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) ! the following params are not used in this module. Instead, they are used in ! the cap but are logged here to keep all the CFC cap params together. From 197cd187042766e059d15ecf7726c0920866cfcd Mon Sep 17 00:00:00 2001 From: Keith Lindsay Date: Tue, 24 Aug 2021 14:32:33 -0600 Subject: [PATCH 25/28] MOM_CFC_cap cleanup remove unused and unneeded local variables correct unit-related comments rework CFC_cap_fluxes to better align with Orr et al., GMD, 2017 --- src/tracer/MOM_CFC_cap.F90 | 114 +++++++++++++++---------------------- 1 file changed, 46 insertions(+), 68 deletions(-) diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 959348b803..e1770b0d52 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -84,9 +84,6 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! This include declares and sets the variable "version". #include "version_variable.h" real, dimension(:,:,:), pointer :: tr_ptr => NULL() - real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients - real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and - real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers. character(len=200) :: dummy ! Dummy variable to store params that need to be logged here. logical :: register_CFC_cap integer :: isd, ied, jsd, jed, nz, m @@ -442,11 +439,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & - CFC11_Csurf, & ! The CFC-11 surface concentrations times the Schmidt number term [mol kg-1]. - CFC12_Csurf, & ! The CFC-12 surface concentrations times the Schmidt number term [mol kg-1]. - CFC11_alpha, & ! The CFC-11 solubility [mol kg-1 atm-1]. - CFC12_alpha, & ! The CFC-12 solubility [mol kg-1 atm-1]. kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1]. + kw, & ! gas transfer velocity [Z T-1 ~> m s-1]. cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) ! [mol kg-1]. cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] @@ -456,10 +450,10 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12. - real :: sc_no_term ! A term related to the Schmidt number. - real :: kw_coeff ! A coefficient used to scale the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + real :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] real :: Rho0_kg_m3 ! Rho0 in non-scaled units [kg m-3] real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. + real :: press_to_atm ! converts from model pressure units to atm integer :: i, j, m, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -493,7 +487,9 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id ! = 6.97e-7 m/s s^2/m^2 [Z T-1 T2 L-2 = Z T L-2 ~> s / m] kw_coeff = (US%m_to_Z*US%s_to_T*US%L_to_m**2) * 6.97e-7 + ! set unit conversion factors Rho0_kg_m3 = Rho0 * US%R_to_kg_m3 + press_to_atm = US%R_to_kg_m3*US%L_T_to_m_s**2 * pa_to_atm do j=js,je ; do i=is,ie ! ta in hectoKelvin @@ -506,77 +502,59 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id ! Calculate Schmidt numbers using coefficients given by ! Wanninkhof (2014); doi:10.4319/lom.2014.12.351. call comp_CFC_schmidt(sfc_state%SST(i,j), sc_11, sc_12) - sc_no_term = sqrt(660.0 / sc_11) - CFC11_alpha(i,j) = alpha_11 * sc_no_term - CFC11_Csurf(i,j) = sfc_state%sfc_CFC11(i,j) * sc_no_term - sc_no_term = sqrt(660.0 / sc_12) - CFC12_alpha(i,j) = alpha_12 * sc_no_term - CFC12_Csurf(i,j) = sfc_state%sfc_CFC12(i,j) * sc_no_term - - kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) + kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) ! air concentrations and cfcs BC's fluxes - ! CFC flux units: mol kg-1 s-1 kg m-2 - cair(i,j) = US%R_to_kg_m3*US%L_T_to_m_s**2*pa_to_atm * CFC11_alpha(i,j) * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc11_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC11_Csurf(i,j)) * Rho0_kg_m3 - cair(i,j) = US%R_to_kg_m3*US%L_T_to_m_s**2*pa_to_atm * CFC12_alpha(i,j) * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) - fluxes%cfc12_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC12_Csurf(i,j)) * Rho0_kg_m3 + ! CFC flux units: CU Z T-1 kg m-3 = mol kg-1 Z T-1 kg m-3 ~> mol m-2 s-1 + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_11) + cair(i,j) = press_to_atm * alpha_11 * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC11(i,j)) * Rho0_kg_m3 + + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_12) + cair(i,j) = press_to_atm * alpha_12 * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC12(i,j)) * Rho0_kg_m3 enddo ; enddo end subroutine CFC_cap_fluxes !> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32. subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask) - real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1] - real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1] - real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin] - real, intent(in ) :: sal !< Surface salinity [PSU]. - real, intent(in ) :: mask !< ocean mask + real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1] + real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1] + real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin] + real, intent(in ) :: sal !< Surface salinity [PSU]. + real, intent(in ) :: mask !< ocean mask ! Local variables - real :: d11_dflt(4), d12_dflt(4) ! values of the various coefficients - real :: e11_dflt(3), e12_dflt(3) ! in the expressions for the solubility - real :: d1_11, d1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [nondim] - real :: d2_11, d2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-1] - real :: d3_11, d3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [log(hectoKelvin)-1] - real :: d4_11, d4_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-2] - real :: e1_11, e1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1] - real :: e2_11, e2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1 hectoKelvin-1] - real :: e3_11, e3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-2 hectoKelvin-2] - real :: factor ! factor to use in the solubility conversion - - !----------------------------------------------------------------------- - ! Solubility coefficients for alpha in mol/(kg atm) for CFC11 (_11) and CFC12 (_12) + + ! Coefficients for calculating CFC11 solubilities + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + + real, parameter :: d1_11 = -232.0411 ! [nondim] + real, parameter :: d2_11 = 322.5546 ! [hectoKelvin-1] + real, parameter :: d3_11 = 120.4956 ! [log(hectoKelvin)-1] + real, parameter :: d4_11 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_11 = -0.146531 ! [PSU-1] + real, parameter :: e2_11 = 0.093621 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_11 = -0.0160693 ! [PSU-2 hectoKelvin-2] + + ! Coefficients for calculating CFC12 solubilities ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. - !----------------------------------------------------------------------- - d11_dflt(:) = (/ -232.0411, 322.5546, 120.4956, -1.39165 /) - e11_dflt(:) = (/ -0.146531, 0.093621, -0.0160693 /) - d12_dflt(:) = (/ -220.2120, 301.8695, 114.8533, -1.39165 /) - e12_dflt(:) = (/ -0.147718, 0.093175, -0.0157340 /) - - d1_11 = d11_dflt(1) - d2_11 = d11_dflt(2) - d3_11 = d11_dflt(3) - d4_11 = d11_dflt(4) - - e1_11 = e11_dflt(1) - e2_11 = e11_dflt(2) - e3_11 = e11_dflt(3) - - d1_12 = d12_dflt(1) - d2_12 = d12_dflt(2) - d3_12 = d12_dflt(3) - d4_12 = d12_dflt(4) - - e1_12 = e12_dflt(1) - e2_12 = e12_dflt(2) - e3_12 = e12_dflt(3) - - ! Calculate solubilities using Warner and Weiss (1985) DSR, vol 32. - ! The following is from Eq. 9 in Warner and Weiss (1985) - ! The factor 1.0e+03 is for the conversion from mol/(l * atm) to mol/(m3 * atm) - ! The factor 1.e-09 converts from mol/(l * atm) to mol/(m3 * pptv) + + real, parameter :: d1_12 = -220.2120 ! [nondim] + real, parameter :: d2_12 = 301.8695 ! [hectoKelvin-1] + real, parameter :: d3_12 = 114.8533 ! [log(hectoKelvin)-1] + real, parameter :: d4_12 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_12 = -0.147718 ! [PSU-1] + real, parameter :: e2_12 = 0.093175 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_12 = -0.0157340 ! [PSU-2 hectoKelvin-2] + + real :: factor ! introduce units to result [mol kg-1 atm-1] + + ! Eq. 9 from Warner and Weiss (1985) DSR, vol 32. factor = 1.0 alpha_11 = exp(d1_11 + d2_11/ta + d3_11*log(ta) + d4_11*ta**2 +& sal * ((e3_11 * ta + e2_11) * ta + e1_11)) * & From 88b36422fb1d61585b671cea5c5d7c64528fd0ee Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 26 Aug 2021 18:37:10 -0600 Subject: [PATCH 26/28] Restore Vt2 diagnostic calculation when requested and retain the functionality of prescribed constant enhancement. --- .../vertical/MOM_CVMix_KPP.F90 | 28 ++++++++++--------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 0abc5d1b2a..a26f251a2b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -420,6 +420,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) select case ( trim(string)) case ("CONSTANT") CS%LT_K_METHOD = LT_K_MODE_CONSTANT + langmuir_mixing_opt = 'LWF16' case ("VR12") CS%LT_K_METHOD = LT_K_MODE_VR12 langmuir_mixing_opt = 'LWF16' @@ -452,6 +453,7 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) select case ( trim(string)) case ("CONSTANT") CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT + langmuir_entrainment_opt = 'LWF16' case ("VR12") CS%LT_VT2_METHOD = LT_VT2_MODE_VR12 langmuir_entrainment_opt = 'LWF16' @@ -856,18 +858,6 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & Kviscosity(:) = 0.0 endif - - ! compute unresolved squared velocity for diagnostics - if (CS%id_Vt2 > 0) then -!BGR Now computing VT2 above so can modify for LT -! therefore, don't repeat this operation here -! CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & -! cellHeight(1:GV%ke), & ! Depth of cell center [m] -! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] -! N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] -! CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters - endif - ! Copy 1d data into 3d diagnostic arrays !/ grabbing obldepth_0d for next time step. CS%OBLdepthprev(i,j)=CS%OBLdepth(i,j) @@ -915,7 +905,6 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%id_OBLdepth_original > 0) call post_data(CS%id_OBLdepth_original,CS%OBLdepth_original,CS%diag) if (CS%id_sigma > 0) call post_data(CS%id_sigma, CS%sigma, CS%diag) if (CS%id_Ws > 0) call post_data(CS%id_Ws, CS%Ws, CS%diag) - if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) if (CS%id_uStar > 0) call post_data(CS%id_uStar, uStar, CS%diag) if (CS%id_buoyFlux > 0) call post_data(CS%id_buoyFlux, buoyFlux, CS%diag) if (CS%id_Kt_KPP > 0) call post_data(CS%id_Kt_KPP, CS%Kt_KPP, CS%diag) @@ -1236,6 +1225,18 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + ! compute unresolved squared velocity for diagnostics + if (CS%id_Vt2 > 0) then + CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( & + cellHeight(1:GV%ke), & ! Depth of cell center [m] + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1] + N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1] + EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim] + LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim] + bfsfc=surfBuoyFlux, & ! surface buoyancy flux [m2 s-3] + uStar=uStar(i,j), & ! surface friction velocity [m s-1] + CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters + endif ! recompute wscale for diagnostics, now that we in fact know boundary layer depth !BGR consider if LTEnhancement is wanted for diagnostics @@ -1278,6 +1279,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%id_EnhK > 0) call post_data(CS%id_EnhK, CS%EnhK, CS%diag) if (CS%id_EnhVt2 > 0) call post_data(CS%id_EnhVt2, CS%EnhVt2, CS%diag) if (CS%id_La_SL > 0) call post_data(CS%id_La_SL, CS%La_SL, CS%diag) + if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag) ! BLD smoothing: if (CS%n_smooth > 0) call KPP_smooth_BLD(CS,G,GV,h) From 11c550f55a2da9b60577d40dc6f7504e580a71aa Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 1 Sep 2021 18:32:36 -0600 Subject: [PATCH 27/28] add references to KPP enhancements methods --- src/parameterizations/vertical/MOM_CVMix_KPP.F90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index fc8525ab4e..4dcaa70bc2 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -411,11 +411,14 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) "Unrecognized KPP_LT_K_SHAPE option: "//trim(string)) end select call get_param(paramFile, mdl, "KPP_LT_K_METHOD", string , & - 'Method to enhance mixing coefficient in KPP. '// & + 'Method to enhance mixing coefficient in KPP. '// & 'Valid options are: \n'// & '\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'// & '\t VR12 = Function of Langmuir number based on VR12\n'// & - '\t RW16 = Function of Langmuir number based on RW16', & + '\t (Van Roekel et al. 2012)\n'// & + '\t (Li et al. 2016, OM) \n'// & + '\t RW16 = Function of Langmuir number based on RW16\n'// & + '\t (Reichl et al., 2016, JPO)', & default='CONSTANT') select case ( trim(string)) case ("CONSTANT") @@ -443,12 +446,16 @@ logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves) 'in Bulk Richardson Number.', units="", Default=.false.) if (CS%LT_Vt2_Enhancement) then call get_param(paramFile, mdl, "KPP_LT_VT2_METHOD",string , & - 'Method to enhance Vt2 in KPP. '// & + 'Method to enhance Vt2 in KPP. '// & 'Valid options are: \n'// & '\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'// & '\t VR12 = Function of Langmuir number based on VR12\n'// & + '\t (Van Roekel et al., 2012) \n'// & + '\t (Li et al. 2016, OM) \n'// & '\t RW16 = Function of Langmuir number based on RW16\n'// & - '\t LF17 = Function of Langmuir number based on LF17', & + '\t (Reichl et al., 2016, JPO) \n'// & + '\t LF17 = Function of Langmuir number based on LF17\n'// & + '\t (Li and Fox-Kemper, 2017, JPO)', & default='CONSTANT') select case ( trim(string)) case ("CONSTANT") From 5f8c446d91d6942a78fc31322a1c690ec7e28c79 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 1 Sep 2021 18:34:00 -0600 Subject: [PATCH 28/28] update all calls of KPP_compute_BLD and KPP_calculate to support the passage of lamult enhancement factor --- .../vertical/MOM_diabatic_driver.F90 | 32 ++++++++++++++----- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b44d4c84b4..01cf7f52e4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -624,11 +624,19 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) ! The KPP scheme calculates boundary layer diffusivities and non-local transport. - call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + if ( associated(fluxes%lamult) ) then + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + else + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + endif if (associated(Hml)) then call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US) @@ -1771,11 +1779,19 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e enddo ; enddo ; enddo endif - call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & - fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) + if ( associated(fluxes%lamult) ) then + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves, lamult=fluxes%lamult) + + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves, lamult=fluxes%lamult) + else + call KPP_compute_BLD(CS%KPP_CSp, G, GV, US, h, tv%T, tv%S, u, v, tv, & + fluxes%ustar, CS%KPP_buoy_flux, Waves=Waves) - call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & - Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + call KPP_calculate(CS%KPP_CSp, G, GV, US, h, fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, & + Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar, Waves=Waves) + endif if (associated(Hml)) then call KPP_get_BLD(CS%KPP_CSp, Hml(:,:), G, US)