This subroutine provides an interface for a user to use to modify the main code to alter the diffusivities as needed. The specific example implemented here augments the diffusivity for a specified range of latitude and coordinate potential density.
48 type(ocean_grid_type),
intent(in) :: G
49 type(verticalGrid_type),
intent(in) :: GV
50 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
51 type(thermo_var_ptrs),
intent(in) :: tv
54 type(user_change_diff_CS),
pointer :: CS
55 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(inout) :: Kd_lay
57 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(inout) :: Kd_int
59 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: T_f
61 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
optional,
intent(in) :: S_f
63 real,
dimension(:,:,:),
optional,
pointer :: Kd_int_add
67 real :: Rcv(SZI_(G),SZK_(G))
68 real :: p_ref(SZI_(G))
73 logical :: store_Kd_add
74 integer :: i, j, k, is, ie, js, je, nz
75 integer :: isd, ied, jsd, jed
79 character(len=200) :: mesg
81 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
82 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
84 if (.not.
associated(cs))
call mom_error(fatal,
"user_set_diffusivity: "//&
85 "Module must be initialized before it is used.")
87 use_eos =
associated(tv%eqn_of_state)
88 if (.not.use_eos)
return
89 store_kd_add = .false.
90 if (
present(kd_int_add)) store_kd_add =
associated(kd_int_add)
92 if (.not.range_ok(cs%lat_range))
then
93 write(mesg,
'(4(1pe15.6))') cs%lat_range(1:4)
94 call mom_error(fatal,
"user_set_diffusivity: bad latitude range: \n "//&
97 if (.not.range_ok(cs%rho_range))
then
98 write(mesg,
'(4(1pe15.6))') cs%rho_range(1:4)
99 call mom_error(fatal,
"user_set_diffusivity: bad density range: \n "//&
103 if (store_kd_add) kd_int_add(:,:,:) = 0.0
105 do i=is,ie ; p_ref(i) = tv%P_Ref ;
enddo
107 if (
present(t_f) .and.
present(s_f))
then
109 call calculate_density(t_f(:,j,k),s_f(:,j,k),p_ref,rcv(:,k),&
110 is,ie-is+1,tv%eqn_of_state)
114 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p_ref,rcv(:,k),&
115 is,ie-is+1,tv%eqn_of_state)
119 if (
present(kd_lay))
then
120 do k=1,nz ;
do i=is,ie
121 if (cs%use_abs_lat)
then
122 lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
124 lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
126 rho_fn = val_weights(rcv(i,k), cs%rho_range)
127 if (rho_fn * lat_fn > 0.0) &
128 kd_lay(i,j,k) = kd_lay(i,j,k) + cs%Kd_add * rho_fn * lat_fn
131 if (
present(kd_int))
then
132 do k=2,nz ;
do i=is,ie
133 if (cs%use_abs_lat)
then
134 lat_fn = val_weights(abs(g%geoLatT(i,j)), cs%lat_range)
136 lat_fn = val_weights(g%geoLatT(i,j), cs%lat_range)
139 rho_fn = val_weights( 0.5*(rcv(i,k-1) + rcv(i,k)), cs%rho_range)
140 if (rho_fn * lat_fn > 0.0)
then
141 kd_int(i,j,k) = kd_int(i,j,k) + cs%Kd_add * rho_fn * lat_fn
142 if (store_kd_add) kd_int_add(i,j,k) = cs%Kd_add * rho_fn * lat_fn