From 83e50eb0dfe8454a0431a079bf716c1bac9a9b37 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 13 Jun 2016 20:26:35 -0400 Subject: [PATCH] Split out DOME OBC positions - Previously, s/r DOME_set_Open_Bdry_Conds() set both the OBC masks and the OBC data (inflow conditions). - DOME_set_OBC_positions() now sets the OBC masks. - DOME_set_OBC_data() replaces DOME_set_Open_Bdry_Conds() and only sets inflow data. - This is one step in a larger re-factor of OBC code. --- .../MOM_state_initialization.F90 | 5 +- src/user/DOME_initialization.F90 | 138 ++++++++---------- 2 files changed, 67 insertions(+), 76 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f9568cd2c..fe310f51a3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -42,7 +42,7 @@ module MOM_state_initialization use user_initialization, only : user_init_temperature_salinity use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness -use DOME_initialization, only : DOME_set_Open_Bdry_Conds +use DOME_initialization, only : DOME_set_OBC_positions, DOME_set_OBC_data use DOME_initialization, only : DOME_initialize_sponges use ISOMIP_initialization, only : ISOMIP_initialize_thickness use ISOMIP_initialization, only : ISOMIP_initialize_sponges @@ -440,7 +440,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & " \t USER - call a user modified routine.", default="file", & fail_if_missing=.true.) if (trim(config) == "DOME") then - call DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) + call DOME_set_OBC_positions(G, PF, OBC) + call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) elseif (trim(config) == "USER") then call user_set_Open_Bdry_Conds(OBC, tv, G, PF, tracer_Reg) else diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index b94fb69a18..455cb559c5 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -37,7 +37,8 @@ module DOME_initialization public DOME_initialize_topography public DOME_initialize_thickness public DOME_initialize_sponges -public DOME_set_Open_Bdry_Conds +public DOME_set_OBC_positions +public DOME_set_OBC_data contains @@ -230,31 +231,69 @@ subroutine DOME_initialize_sponges(G, GV, tv, PF, CSp) endif end subroutine DOME_initialize_sponges -! ----------------------------------------------------------------------------- -! ----------------------------------------------------------------------------- +!> Set the positions of the open boundary needed for the DOME experiment. +subroutine DOME_set_OBC_positions(G, param_file, OBC) + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(param_file_type), intent(in) :: param_file !< Parameter file handle. + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure. + ! Local variables + character(len=40) :: mod = "DOME_set_OBC_positions" ! This subroutine's name. + integer :: i, j + logical :: any_OBC ! Set to true if any points in this subdomain use OBCs + + if (.not.associated(OBC)) allocate(OBC) + + call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & + "If true, open boundary conditions may be set at some \n"//& + "u-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & + "If true, open boundary conditions may be set at some \n"//& + "v-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + + any_OBC = .false. + if (OBC%apply_OBC_u) then + ! Set where u points are determined by OBCs. + !allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. + call MOM_error(FATAL,"DOME_initialization, DOME_set_OBC_positions: "//& + "APPLY_OBC_U=True is not coded for the DOME experiment") + endif + if (OBC%apply_OBC_v) then + ! Set where v points are determined by OBCs. + allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. + do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + if ((G%geoLonCv(i,J) > 1000.0) .and. (G%geoLonCv(i,J) < 1100.0) .and. & + (abs(G%geoLatCv(i,J) - G%gridLatB(G%JegB)) < 0.1)) then + OBC%OBC_mask_v(i,J) = .true. + any_OBC = .true. + endif + enddo ; enddo + endif + if (.not.any_OBC) then + ! If this PE does not have any OBC points then we do not need the mask + OBC%apply_OBC_v = .false. + deallocate(OBC%OBC_mask_v) + endif +end subroutine DOME_set_OBC_positions + !> This subroutine sets the properties of flow at open boundary conditions. !! This particular example is for the DOME inflow describe in Legg et al. 2006. -subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) +subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies !! whether, where, and what open boundary !! conditions are used. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any !! available thermodynamic fields, including potential - !! temperature and salinity or mixed layer density. Absent - !! fields have NULL ptrs. + !! temperature and salinity or mixed layer density. Absent + !! fields have NULL ptrs. 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(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. - logical :: any_OBC ! Set to true if any points in this subdomain use - ! open boundary conditions. - logical, pointer, dimension(:,:) :: & - OBC_mask_u => NULL(), & ! These arrays are true at zonal or meridional - OBC_mask_v => NULL() ! velocity points that have prescribed open boundary - ! conditions. real, pointer, dimension(:,:,:) :: & OBC_T_u => NULL(), & ! These arrays should be allocated and set to OBC_T_v => NULL(), & ! specify the values of T and S that should come @@ -276,7 +315,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! thickness D_edge, in the same units as lat. real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile. - character(len=40) :: mod = "DOME_set_Open_Bdry_Conds" ! This subroutine's name. + character(len=40) :: mod = "DOME_set_OBC_data" ! This subroutine's name. integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB @@ -289,56 +328,10 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) Ri_trans = 1.0/3.0 ! The shear Richardson number in the transition region ! region of the specified shear profile. - call get_param(param_file, mod, "APPLY_OBC_U", apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V", apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - - if (apply_OBC_u) then - ! Determine where u points are applied. - allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. - any_OBC = .false. - do j=jsd,jed ; do I=IsdB,IedB - ! if (SOME_TEST_FOR_U_OPEN_BCS) then - ! OBC_mask_u(I,j) = .true. ; any_OBC = .true. - ! endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no u points at which open boundary conditions are - ! to be applied. - apply_OBC_u = .false. - deallocate(OBC_mask_u) - endif - endif - if (apply_OBC_v) then - ! Determine where v points are applied. - allocate(OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC_mask_v(:,:) = .false. - any_OBC = .false. - do J=JsdB,JedB ; do i=isd,ied - if ((G%geoLonCv(i,J) > 1000.0) .and. (G%geoLonCv(i,J) < 1100.0) .and. & - (abs(G%geoLatCv(i,J) - G%gridLatB(G%JegB)) < 0.1)) then - OBC_mask_v(i,J) = .true. ; any_OBC = .true. - endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no v points at which open boundary conditions are - ! to be applied. - apply_OBC_v = .false. - deallocate(OBC_mask_v) - endif - endif - - if (.not.(apply_OBC_u .or. apply_OBC_v)) return + if (.not.associated(OBC)) return + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v)) return - if (.not.associated(OBC)) allocate(OBC) - - if (apply_OBC_u) then - OBC%apply_OBC_u = .true. - OBC%OBC_mask_u => OBC_mask_u + if (OBC%apply_OBC_u) then allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE @@ -346,9 +339,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE enddo ; enddo endif - if (apply_OBC_v) then - OBC%apply_OBC_v = .true. - OBC%OBC_mask_v => OBC_mask_v + if (OBC%apply_OBC_v) then allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE @@ -357,7 +348,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) enddo ; enddo endif - if (apply_OBC_v) then + if (OBC%apply_OBC_v) then g_prime_tot = (G%g_Earth/GV%Rho0)*2.0 Def_Rad = sqrt(D_edge*g_prime_tot) / (1.0e-4*1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*Def_Rad) * GV%m_to_H @@ -380,7 +371,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) if (k == nz) tr_k = tr_k + tr_0 * (2.0/(Ri_trans*(2.0+Ri_trans))) * & log((2.0+Ri_trans)/(2.0-Ri_trans)) do J=JsdB,JedB ; do i=isd,ied - if (OBC_mask_v(i,J)) then + if (OBC%OBC_mask_v(i,J)) then ! This needs to be unneccesarily complicated without symmetric memory. lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J) ! if (isd > IsdB) lon_im1 = G%geoLonBu(I-1,J) @@ -394,9 +385,9 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) enddo endif - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - if (OBC_mask_u(I,j)) then + if (OBC%OBC_mask_u(I,j)) then ! An appropriate expression for the zonal inflow velocities and ! transports should go here. OBC%uh(I,j,k) = 0.0 * GV%m_to_H ; OBC%u(I,j,k) = 0.0 @@ -408,7 +399,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! The inflow values of temperature and salinity also need to be set here if ! these variables are used. The following code is just a naive example. - if (apply_OBC_u .or. apply_OBC_v) then + if (OBC%apply_OBC_u .or. OBC%apply_OBC_v) then if (associated(tv%S)) then ! In this example, all S inflows have values of 35 psu. call add_tracer_OBC_values("S", tr_Reg, OBC_inflow=35.0) @@ -428,13 +419,13 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) do k=1,nz ; T0(k) = T0(k) + (GV%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo enddo - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB OBC_T_u(I,j,k) = T0(k) enddo ; enddo ; enddo endif - if (apply_OBC_v) then + if (OBC%apply_OBC_v) then allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied OBC_T_v(i,J,k) = T0(k) @@ -445,8 +436,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) endif endif -end subroutine DOME_set_Open_Bdry_Conds -! ----------------------------------------------------------------------------- +end subroutine DOME_set_OBC_data !> \class DOME_initialization !!