diff --git a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 index 3c8f084b4a..7075fb7c10 100644 --- a/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing_gfdl.F90 @@ -189,12 +189,6 @@ module MOM_surface_forcing_gfdl !! ice-shelves, expressed as a coefficient !! for divergence damping, as determined !! outside of the ocean model [m3 s-1] - real, pointer, dimension(:,:) :: ustk0 => NULL() !< - real, pointer, dimension(:,:) :: vstk0 => NULL() !< - real, pointer, dimension(:) :: stk_wavenumbers => NULL() !< - real, pointer, dimension(:,:,:) :: ustkb => NULL() !< - real, pointer, dimension(:,:,:) :: vstkb => NULL() !< - integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of named fields !! used for passive tracer fluxes. @@ -202,7 +196,6 @@ module MOM_surface_forcing_gfdl !! This flag may be set by the flux-exchange code, based on what !! the sea-ice model is providing. Otherwise, the value from !! the surface_forcing_CS is used. - integer :: num_stk_bands !< Number of Stokes drift bands passed through the coupler end type ice_ocean_boundary_type integer :: id_clock_forcing !< A CPU time clock @@ -684,7 +677,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ real :: mass_eff ! effective mass of sea ice for rigidity [R Z ~> kg m-2] real :: wt1, wt2 ! Relative weights of previous and current values of ustar, ND. - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0, istk + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd @@ -725,9 +718,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ (associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) & call allocate_mech_forcing(G, forces, iceberg=.true.) - if ( associated(IOB%ustk0) ) & - call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands) - if (associated(IOB%ice_rigidity)) then rigidity_at_h(:,:) = 0.0 call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) @@ -787,19 +777,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_ forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j) enddo ; enddo endif - forces%stk_wavenumbers(:) = IOB%stk_wavenumbers - do j=js,je; do i=is,ie - forces%ustk0(i,j) = IOB%ustk0(i-I0,j-J0) ! How to be careful here that the domains are right? - forces%vstk0(i,j) = IOB%vstk0(i-I0,j-J0) - enddo ; enddo - call pass_vector(forces%ustk0,forces%vstk0, G%domain ) - do istk = 1,IOB%num_stk_bands - do j=js,je; do i=is,ie - forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk) - forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk) - enddo; enddo - call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain ) - enddo ! Find the net mass source in the input forcing without other adjustments. if (CS%approx_net_mass_src .and. associated(forces%net_mass_src)) then diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 223c179bfb..b091af9bb0 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -566,7 +566,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda ! For now, the waves are only updated on the thermodynamics steps, because that is where ! the wave intensities are actually used to drive mixing. At some point, the wave updates ! might also need to become a part of the ocean dynamics, according to B. Reichl. - call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces) + call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves) endif if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model. diff --git a/config_src/mct_driver/mom_surface_forcing_mct.F90 b/config_src/mct_driver/mom_surface_forcing_mct.F90 index 82105e040e..ef0527dd1c 100644 --- a/config_src/mct_driver/mom_surface_forcing_mct.F90 +++ b/config_src/mct_driver/mom_surface_forcing_mct.F90 @@ -486,17 +486,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! latent heat flux (W/m^2) fluxes%latent(i,j) = 0.0 - ! contribution from frozen ppt + ! contribution from frozen ppt (notice minus sign since fprec is positive into the ocean) if (associated(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif - ! contribution from frozen runoff + ! contribution from frozen runoff (notice minus sign since rofi_flux is positive into the ocean) if (associated(fluxes%frunoff)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif ! contribution from evaporation if (associated(IOB%q_flux)) then diff --git a/config_src/nuopc_driver/mom_cap.F90 b/config_src/nuopc_driver/mom_cap.F90 index e078c0d74f..3bd81b5fcd 100644 --- a/config_src/nuopc_driver/mom_cap.F90 +++ b/config_src/nuopc_driver/mom_cap.F90 @@ -36,6 +36,7 @@ module MOM_cap_mod use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh use MOM_cap_time, only: AlarmInit use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, state_diagnose +use MOM_cap_methods, only: ChkErr #ifdef CESMCOUPLED use shr_file_mod, only: shr_file_setLogUnit, shr_file_getLogUnit #endif @@ -2060,18 +2061,6 @@ subroutine shr_file_getLogUnit(nunit) end subroutine shr_file_getLogUnit #endif -logical function chkerr(rc, line, file) - integer, intent(in) :: rc - integer, intent(in) :: line - character(len=*), intent(in) :: file - integer :: lrc - chkerr = .false. - lrc = rc - if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then - chkerr = .true. - endif -end function chkerr - !> !! @page nuopc_cap NUOPC Cap !! @author Fei Liu (fei.liu@gmail.com) diff --git a/config_src/nuopc_driver/mom_cap_methods.F90 b/config_src/nuopc_driver/mom_cap_methods.F90 index 1d51c1e6dd..78014f1c63 100644 --- a/config_src/nuopc_driver/mom_cap_methods.F90 +++ b/config_src/nuopc_driver/mom_cap_methods.F90 @@ -30,6 +30,7 @@ module MOM_cap_methods public :: mom_import public :: mom_export public :: state_diagnose +public :: ChkErr private :: State_getImport private :: State_setExport @@ -251,9 +252,9 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- - ! Partitioned Stokes Drift Components + ! Partitioned Stokes Drift Components !---- - if ( associated(ice_ocean_boundary%ustkb) ) then + if ( associated(ice_ocean_boundary%ustkb) ) then allocate(stkx1(isc:iec,jsc:jec)) allocate(stky1(isc:iec,jsc:jec)) allocate(stkx2(isc:iec,jsc:jec)) @@ -765,15 +766,18 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid end subroutine State_SetExport +!> This subroutine writes the minimum, maximum and sum of each field +!! contained within an ESMF state. subroutine state_diagnose(State, string, rc) ! ---------------------------------------------- ! Diagnose status of State ! ---------------------------------------------- - type(ESMF_State), intent(in) :: state - character(len=*), intent(in) :: string - integer , intent(out) :: rc + type(ESMF_State), intent(in) :: state !< An ESMF State + character(len=*), intent(in) :: string !< A string indicating whether the State is an + !! import or export State + integer , intent(out) :: rc !< Return code ! local variables integer :: i,j,n @@ -787,19 +791,19 @@ subroutine state_diagnose(State, string, rc) ! ---------------------------------------------- call ESMF_StateGet(state, itemCount=fieldCount, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate(lfieldnamelist(fieldCount)) call ESMF_StateGet(state, itemNameList=lfieldnamelist, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return do n = 1, fieldCount call ESMF_StateGet(state, itemName=lfieldnamelist(n), field=lfield, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call field_getfldptr(lfield, fldptr1=dataPtr1d, fldptr2=dataPtr2d, rank=lrank, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return if (lrank == 0) then ! no local data @@ -829,23 +833,16 @@ subroutine state_diagnose(State, string, rc) end subroutine state_diagnose -!=============================================================================== - +!> Obtain a pointer to a rank 1 or rank 2 ESMF field subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) - ! ---------------------------------------------- - ! for a field, determine rank and return fldptr1 or fldptr2 - ! abort is true by default and will abort if fldptr is not yet allocated in field - ! rank returns 0, 1, or 2. 0 means fldptr not allocated and abort=false - ! ---------------------------------------------- - ! input/output variables - type(ESMF_Field) , intent(in) :: field - real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:) - real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:) - integer , intent(out) , optional :: rank - logical , intent(in) , optional :: abort - integer , intent(out) , optional :: rc + type(ESMF_Field) , intent(in) :: field !< An ESMF field + real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:) !< A pointer to a rank 1 ESMF field + real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:) !< A pointer to a rank 2 ESMF field + integer , intent(out) , optional :: rank !< Field rank + logical , intent(in) , optional :: abort !< Abort code + integer , intent(out) , optional :: rc !< Return code ! local variables type(ESMF_GeomType_Flag) :: geomtype @@ -872,7 +869,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) lrank = -99 call ESMF_FieldGet(field, status=status, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return if (status /= ESMF_FIELDSTATUS_COMPLETE) then lrank = 0 @@ -886,20 +883,20 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) else call ESMF_FieldGet(field, geomtype=geomtype, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return if (geomtype == ESMF_GEOMTYPE_GRID) then call ESMF_FieldGet(field, rank=lrank, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return elseif (geomtype == ESMF_GEOMTYPE_MESH) then call ESMF_FieldGet(field, rank=lrank, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_FieldGet(field, mesh=lmesh, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_MeshGet(lmesh, numOwnedNodes=nnodes, numOwnedElements=nelements, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return if (nnodes == 0 .and. nelements == 0) lrank = 0 - else + else call ESMF_LogWrite(trim(subname)//": ERROR geomtype not supported ", & ESMF_LOGMSG_INFO) rc = ESMF_FAILURE @@ -917,7 +914,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) return endif call ESMF_FieldGet(field, farrayPtr=fldptr1, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return elseif (lrank == 2) then if (.not.present(fldptr2)) then call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", & @@ -926,7 +923,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) return endif call ESMF_FieldGet(field, farrayPtr=fldptr2, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + if (ChkErr(rc,__LINE__,u_FILE_u)) return else call ESMF_LogWrite(trim(subname)//": ERROR in rank ", & ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u) @@ -942,16 +939,17 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc) end subroutine field_getfldptr -logical function chkerr(rc, line, file) - integer, intent(in) :: rc - integer, intent(in) :: line - character(len=*), intent(in) :: file +!> Returns true if ESMF_LogFoundError() determines that rc is an error code. Otherwise false. +logical function ChkErr(rc, line, file) + integer, intent(in) :: rc !< return code to check + integer, intent(in) :: line !< Integer source line number + character(len=*), intent(in) :: file !< User-provided source file name integer :: lrc - chkerr = .false. + ChkErr = .false. lrc = rc if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then - chkerr = .true. + ChkErr = .true. endif -end function chkerr +end function ChkErr end module MOM_cap_methods diff --git a/config_src/nuopc_driver/mom_cap_time.F90 b/config_src/nuopc_driver/mom_cap_time.F90 index aa1a6c7072..7f210bda71 100644 --- a/config_src/nuopc_driver/mom_cap_time.F90 +++ b/config_src/nuopc_driver/mom_cap_time.F90 @@ -16,6 +16,7 @@ module MOM_cap_time use ESMF , only : ESMF_RC_ARG_BAD use ESMF , only : operator(<), operator(/=), operator(+), operator(-), operator(*) , operator(>=) use ESMF , only : operator(<=), operator(>), operator(==) +use MOM_cap_methods , only : ChkErr implicit none; private @@ -336,16 +337,4 @@ subroutine date2ymd (date, year, month, day) end subroutine date2ymd -logical function chkerr(rc, line, file) - integer, intent(in) :: rc - integer, intent(in) :: line - character(len=*), intent(in) :: file - integer :: lrc - chkerr = .false. - lrc = rc - if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then - chkerr = .true. - endif -end function chkerr - end module MOM_cap_time diff --git a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 index ada0a7c302..7168823fbc 100644 --- a/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 +++ b/config_src/nuopc_driver/mom_surface_forcing_nuopc.F90 @@ -183,11 +183,15 @@ 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(:,:) :: ustk0 => NULL() !< - real, pointer, dimension(:,:) :: vstk0 => NULL() !< - real, pointer, dimension(:) :: stk_wavenumbers => NULL() !< - real, pointer, dimension(:,:,:) :: ustkb => NULL() !< - real, pointer, dimension(:,:,:) :: vstkb => NULL() !< + 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] + real, pointer, dimension(:,:,:) :: ustkb => NULL() !< Stokes Drift spectrum, zonal [m/s] + !! Horizontal - u points + !! 3rd dimension - wavenumber + real, pointer, dimension(:,:,:) :: vstkb => NULL() !< Stokes Drift spectrum, meridional [m/s] + !! Horizontal - v points + !! 3rd dimension - wavenumber integer :: num_stk_bands !< Number of Stokes drift bands passed through the coupler integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of @@ -493,15 +497,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(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 - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion endif + ! notice minus sign since frunoff is positive into the ocean if (associated(IOB%frunoff)) then - fluxes%latent(i,j) = fluxes%latent(i,j) + & + fluxes%latent(i,j) = fluxes%latent(i,j) - & IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * & + fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * & IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion endif if (associated(IOB%q_flux)) then @@ -858,7 +864,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS) enddo; enddo call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain ) enddo - endif + endif ! sea ice related dynamic fields if (associated(IOB%ice_rigidity)) then diff --git a/config_src/nuopc_driver/time_utils.F90 b/config_src/nuopc_driver/time_utils.F90 index db056f4bf5..81efcd2765 100644 --- a/config_src/nuopc_driver/time_utils.F90 +++ b/config_src/nuopc_driver/time_utils.F90 @@ -14,6 +14,7 @@ module time_utils_mod use ESMF, only: ESMF_Time, ESMF_TimeGet, ESMF_LogFoundError use ESMF, only: ESMF_LOGERR_PASSTHRU,ESMF_TimeInterval use ESMF, only: ESMF_TimeIntervalGet, ESMF_TimeSet, ESMF_SUCCESS +use MOM_cap_methods, only: ChkErr implicit none; private @@ -160,16 +161,4 @@ function string_to_date(string, rc) end function string_to_date -logical function chkerr(rc, line, file) - integer, intent(in) :: rc - integer, intent(in) :: line - character(len=*), intent(in) :: file - integer :: lrc - chkerr = .false. - lrc = rc - if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then - chkerr = .true. - endif -end function chkerr - end module time_utils_mod diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 6720414b2b..f0cc8f553c 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -250,13 +250,17 @@ module MOM_forcing_type !! reset to zero at the driver level when appropriate. real, pointer, dimension(:,:) :: & - ustk0 => NULL(), & - vstk0 => NULL() + ustk0 => NULL(), & !< Surface Stokes drift, zonal [m/s] + vstk0 => NULL() !< Surface Stokes drift, meridional [m/s] real, pointer, dimension(:) :: & - stk_wavenumbers => NULL() + stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m] real, pointer, dimension(:,:,:) :: & - ustkb => NULL(), & - vstkb => NULL() + ustkb => NULL(), & !< Stokes Drift spectrum, zonal [m/s] + !! Horizontal - u points + !! 3rd dimension - wavenumber + vstkb => NULL() !< Stokes Drift spectrum, meridional [m/s] + !! Horizontal - v points + !! 3rd dimension - wavenumber logical :: initialized = .false. !< This indicates whether the appropriate arrays have been initialized. end type mech_forcing @@ -3003,8 +3007,8 @@ subroutine allocate_mech_forcing_by_group(G, forces, stress, ustar, shelf, & logical, optional, intent(in) :: shelf !< If present and true, allocate forces for ice-shelf logical, optional, intent(in) :: press !< If present and true, allocate p_surf and related fields logical, optional, intent(in) :: iceberg !< If present and true, allocate forces for icebergs - logical, optional, intent(in) :: waves !< If present and true, allocate wave fields - integer, optional, intent(in) :: num_stk_bands !< Number of Stokes bands to allocate + logical, optional, intent(in) :: waves !< If present and true, allocate wave fields + integer, optional, intent(in) :: num_stk_bands !< Number of Stokes bands to allocate ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -3034,19 +3038,19 @@ 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 (present(waves)) then; if (waves) then; if (.not.associated(forces%ustkb)) then if (.not.(present(num_stk_bands))) 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)) + 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 ; endif ; endif - if (present(waves)) then; if (waves) then; if (.not.associated(forces%vstkb)) then + if (present(waves)) then; if (waves) then; if (.not.associated(forces%vstkb)) then allocate(forces%vstkb(isd:ied,jsd:jed,num_stk_bands)) forces%vstkb(isd:ied,jsd:jed,:) = 0.0 - endif; endif; endif + endif ; endif ; endif end subroutine allocate_mech_forcing_by_group diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 2b7a5646cc..01d8e4163b 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -579,13 +579,13 @@ subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_ type(lbd_CS), pointer :: CS !< Lateral diffusion control structure !! the boundary layer ! Local variables - real, dimension(:), allocatable :: dz_top !< The LBD z grid to be created [L ~ m] - real, dimension(:), allocatable :: phi_L_z !< Tracer values in the ztop grid (left) [conc] - real, dimension(:), allocatable :: phi_R_z !< Tracer values in the ztop grid (right) [conc] - real, dimension(:), allocatable :: F_layer_z !< Diffusive flux at U- or V-point in the ztop grid [H L2 conc ~> m3 conc] + real, dimension(:), allocatable :: dz_top !< The LBD z grid to be created [L ~ m] + real, dimension(:), allocatable :: phi_L_z !< Tracer values in the ztop grid (left) [conc] + real, dimension(:), allocatable :: phi_R_z !< Tracer values in the ztop grid (right) [conc] + real, dimension(:), allocatable :: F_layer_z !< Diffusive flux at U/V-point in the ztop grid [H L2 conc ~> m3 conc] real, dimension(ke) :: h_vel !< Thicknesses at u- and v-points in the native grid - !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] - real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] + !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2] + real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1] !! This is just to remind developers that khtr_avg should be !! computed once khtr is 3D. real :: htot !< Total column thickness [H ~> m or kg m-2] diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 4ba1b779e3..0780ca5e3d 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -466,6 +466,12 @@ subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces) if (DataSource==DATAOVR) then call Surface_Bands_by_data_override(day_center, G, GV, US, CS) elseif (DataSource==Coupler) then + if (.not.present(FORCES)) then + call MOM_error(FATAL,"The option SURFBAND = COUPLER can not be used with "//& + "this driver. If you are using a coupled driver with a wave model then "//& + "check the arguments in the subroutine call to Update_Surface_Waves, "//& + "otherwise select another option for SURFBAND_SOURCE.") + endif if (size(CS%WaveNum_Cen).ne.size(forces%stk_wavenumbers)) then call MOM_error(FATAL, "Number of wavenumber bands in WW3 does not match that in MOM6. "//& "Make sure that STK_BAND_COUPLER in MOM6 input is equal to the number of bands in "//& @@ -475,13 +481,13 @@ subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces) do b=1,CS%NumBands CS%WaveNum_Cen(b) = forces%stk_wavenumbers(b) !Interpolate from a grid to c grid - do II=G%iscB,G%iecB - do jj=G%jsc,G%jec + do jj=G%jsc,G%jec + do II=G%iscB,G%iecB CS%STKx0(II,jj,b) = 0.5*(forces%UStkb(ii,jj,b)+forces%UStkb(ii+1,jj,b)) enddo enddo - do ii=G%isc,G%iec - do JJ=G%jscB, G%jecB + do JJ=G%jscB, G%jecB + do ii=G%isc,G%iec CS%STKY0(ii,JJ,b) = 0.5*(forces%VStkb(ii,jj,b)+forces%VStkb(ii,jj+1,b)) enddo enddo @@ -489,13 +495,13 @@ subroutine Update_Surface_Waves(G, GV, US, Day, dt, CS, forces) enddo elseif (DataSource==Input) then do b=1,CS%NumBands - do II=G%isdB,G%iedB - do jj=G%jsd,G%jed + do jj=G%jsd,G%jed + do II=G%isdB,G%iedB CS%STKx0(II,jj,b) = CS%PrescribedSurfStkX(b) enddo enddo - do ii=G%isd,G%ied - do JJ=G%jsdB, G%jedB + do JJ=G%jsdB, G%jedB + do ii=G%isd,G%ied CS%STKY0(ii,JJ,b) = CS%PrescribedSurfStkY(b) enddo enddo @@ -531,8 +537,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) ! Computing mid-point value from surface value and decay wavelength if (WaveMethod==TESTPROF) then DecayScale = 4.*PI / TP_WVL !4pi - do II = G%isdB,G%iedB - do jj = G%jsd,G%jed + do jj = G%jsd,G%jed + do II = G%isdB,G%iedB IIm1 = max(1,II-1) Bottom = 0.0 MidPoint = 0.0 @@ -544,8 +550,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo enddo enddo - do ii = G%isd,G%ied - do JJ = G%jsdB,G%jedB + do JJ = G%jsdB,G%jedB + do ii = G%isd,G%ied JJm1 = max(1,JJ-1) Bottom = 0.0 MidPoint = 0.0 @@ -566,8 +572,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) CS%Us0_x(:,:) = 0.0 CS%Us0_y(:,:) = 0.0 ! Computing X direction Stokes drift - do II = G%isdB,G%iedB - do jj = G%jsd,G%jed + do jj = G%jsd,G%jed + do II = G%isdB,G%iedB ! 1. First compute the surface Stokes drift ! by integrating over the partitionas. do b = 1,CS%NumBands @@ -624,8 +630,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo enddo ! Computing Y direction Stokes drift - do ii = G%isd,G%ied - do JJ = G%jsdB,G%jedB + do JJ = G%jsdB,G%jedB + do ii = G%isd,G%ied ! Compute the surface values. do b = 1,CS%NumBands if (PartitionMode==0) then @@ -682,8 +688,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo elseif (WaveMethod==DHH85) then if (.not.(StaticWaves .and. DHH85_is_set)) then - do II = G%isdB,G%iedB - do jj = G%jsd,G%jed + do jj = G%jsd,G%jed + do II = G%isdB,G%iedB bottom = 0.0 do kk = 1,G%ke Top = Bottom @@ -700,8 +706,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) enddo enddo enddo - do ii = G%isd,G%ied - do JJ = G%jsdB,G%jedB + do JJ = G%jsdB,G%jedB + do ii = G%isd,G%ied Bottom = 0.0 do kk=1, G%ke Top = Bottom @@ -726,13 +732,13 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) endif else! Keep this else, fallback to 0 Stokes drift do kk= 1,G%ke - do II = G%isdB,G%iedB - do jj = G%jsd,G%jed + do jj = G%jsd,G%jed + do II = G%isdB,G%iedB CS%Us_x(II,jj,kk) = 0. enddo enddo - do ii = G%isd,G%ied - do JJ = G%jsdB,G%jedB + do JJ = G%jsdB,G%jedB + do ii = G%isd,G%ied CS%Us_y(ii,JJ,kk) = 0. enddo enddo @@ -742,8 +748,8 @@ subroutine Update_Stokes_Drift(G, GV, US, CS, h, ustar) ! Turbulent Langmuir number is computed here and available to use anywhere. ! SL Langmuir number requires mixing layer depth, and therefore is computed ! in the routine it is needed by (e.g. KPP or ePBL). - do ii = G%isc,G%iec - do jj = G%jsc, G%jec + do jj = G%jsc, G%jec + do ii = G%isc,G%iec Top = h(ii,jj,1)*GV%H_to_Z call get_Langmuir_Number( La, G, GV, US, Top, ustar(ii,jj), ii, jj, & H(ii,jj,:),Override_MA=.false.,WAVES=CS)